diff --git a/dumux/multidomain/common/multidomainassembler.hh b/dumux/multidomain/common/multidomainassembler.hh index f13300fb75be5ef227c869c45ecf088d68a8f659..92762907c90aba3b5890dfc7cd75c57645f6c986 100644 --- a/dumux/multidomain/common/multidomainassembler.hh +++ b/dumux/multidomain/common/multidomainassembler.hh @@ -198,6 +198,13 @@ public: // printvector(std::cout, residual_, "residual", "row", 200, 1, 3); } + /*! + * \brief docme + */ + void reassembleAll() + { + } + /* * \brief docme * diff --git a/dumux/multidomain/common/multidomainmodel.hh b/dumux/multidomain/common/multidomainmodel.hh index db7328b052f2ac3959bcec1f3fda5bbaf567a2d2..c99c7861dabc7948458c88395548b6a19d794cd4 100644 --- a/dumux/multidomain/common/multidomainmodel.hh +++ b/dumux/multidomain/common/multidomainmodel.hh @@ -262,6 +262,17 @@ public: } + /*! + * \brief Check the plausibility of the current solution + * + * This has to be done by the actual model, it knows + * best, what (ranges of) variables to check. + * This is primarily a hook + * which the actual model can overload. + */ + void checkPlausibility() const + { } + /*! * \brief Called by the update() method before it tries to * apply the newton method. This is primary a hook diff --git a/dumux/multidomain/common/multidomainnewtoncontroller.hh b/dumux/multidomain/common/multidomainnewtoncontroller.hh index 2049f0597dd098be0a14aa6cf83d65844eafb9ea..ab34bdf0aa6958175c3ca7177def9c6703abc092 100644 --- a/dumux/multidomain/common/multidomainnewtoncontroller.hh +++ b/dumux/multidomain/common/multidomainnewtoncontroller.hh @@ -23,23 +23,7 @@ #ifndef DUMUX_MULTIDOMAIN_NEWTON_CONTROLLER_HH #define DUMUX_MULTIDOMAIN_NEWTON_CONTROLLER_HH -#include <dumux/common/exceptions.hh> -#include <dumux/linear/linearsolverproperties.hh> -#include <dumux/linear/boxlinearsolver.hh> - -#include <dune/istl/overlappingschwarz.hh> -#include <dune/istl/schwarz.hh> -#include <dune/istl/preconditioners.hh> -#include <dune/istl/solvers.hh> -#include "dune/istl/owneroverlapcopy.hh" - -#include <dune/istl/io.hh> -#include <dune/common/mpihelper.hh> -#include <iostream> -#include <boost/format.hpp> - -#include <dune/pdelab/backend/istlsolverbackend.hh> - +#include <dumux/nonlinear/newtoncontroller.hh> #include "multidomainconvergencewriter.hh" /*! @@ -56,49 +40,11 @@ class MultiDomainNewtonController; namespace Properties { - -//! Specifies the implementation of the Newton controller -NEW_PROP_TAG(NewtonController); - -//! Specifies the type of the actual Newton method -NEW_PROP_TAG(NewtonMethod); - -//! specifies whether the convergence rate and the global residual -//! gets written out to disk for every newton iteration (default is false) -NEW_PROP_TAG(NewtonWriteConvergence); - -/*! - * \brief Specifies whether the update should be done using the line search - * method instead of the plain Newton method. - * - * Whether this property has any effect depends on wether the line - * search method is implemented for the actual model's Newton - * controller's update() method. By default line search is not used. - */ -NEW_PROP_TAG(NewtonUseLineSearch); - -//! the value for the relative error below which convergence is declared -NEW_PROP_TAG(NewtonRelTolerance); - -/*! - * \brief The number of iterations at which the Newton method - * should aim at. - * - * This is used to control the time step size. The heuristic used - * is to scale the last time step size by the deviation of the - * number of iterations used from the target steps. - */ -NEW_PROP_TAG(NewtonTargetSteps); - -//! Number of maximum iterations for the Newton method. -NEW_PROP_TAG(NewtonMaxSteps); - -// set default values for Newton +// set default values for Newton for multidomain problems // they can be overwritten in the parameter file SET_INT_PROP(MultiDomain, NewtonTargetSteps, 8); SET_INT_PROP(MultiDomain, NewtonMaxSteps, 15); SET_SCALAR_PROP(MultiDomain, NewtonRelTolerance, 1e-5); -SET_BOOL_PROP(MultiDomain, NewtonWriteConvergence, false); } @@ -111,8 +57,9 @@ SET_BOOL_PROP(MultiDomain, NewtonWriteConvergence, false); * methods. */ template <class TypeTag> -class MultiDomainNewtonController +class MultiDomainNewtonController : public NewtonController<TypeTag> { + typedef NewtonController<TypeTag> ParentType; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, NewtonController) Implementation; @@ -135,148 +82,26 @@ class MultiDomainNewtonController * \brief docme * \param problem docme */ - public: MultiDomainNewtonController(const Problem &problem) - : endIterMsgStream_(std::ostringstream::out) + : ParentType(problem) + , endIterMsgStream_(std::ostringstream::out) , linearSolver_(problem) , convergenceWriter_(asImp_()) { - verbose_ = true; - numSteps_ = 0; - - // maximum tolerated deflection between two iterations - setRelTolerance(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Newton, RelTolerance)); - setTargetSteps(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Newton, TargetSteps)); - setMaxSteps(GET_PARAM_FROM_GROUP(TypeTag, Scalar, Newton, MaxSteps)); - // Writes out, where the relative tolerance is defined -// std::cout << "ParameterNewtonRelTol= " << PROP_DIAGNOSTIC(TypeTag, NewtonRelTolerance) << std::endl; + std::cout << "ParameterNewtonRelTol= " + << PROP_DIAGNOSTIC(TypeTag, NewtonRelTolerance) + << ", " + << GET_PROP_VALUE(TypeTag, NewtonRelTolerance) + << std::endl; }; /*! * \brief Destructor */ ~MultiDomainNewtonController() - { - }; - - /*! - * \brief Specifies if the newton method ought to be chatty. - * - * \param val docme - * - */ - void setVerbose(bool val) - { verbose_ = val; } - - /*! - * \brief Returns true if the newton method ought to be chatty. - */ - bool verbose() const - { return verbose_ && gridView_().comm().rank() == 0; } - - /*! - * \brief Set the maximum acceptable difference for convergence of - * any primary variable between two iterations. - * - * \param tolerance The maximum relative error between two Newton - * iterations at which the scheme is considered - * finished - */ - void setRelTolerance(Scalar tolerance) - { tolerance_ = tolerance; } - - /*! - * \brief Set the maximum acceptable residual norm reduction. - * - * \param tolerance The maximum reduction of the residual norm - * at which the scheme is considered finished - */ -// void setAbsTolerance(Scalar tolerance) -// { absoluteTolerance_ = tolerance; } - - /*! - * \brief Set the number of iterations at which the Newton method - * should aim at. - * - * This is used to control the time step size. The heuristic used - * is to scale the last time step size by the deviation of the - * number of iterations used from the target steps. - * - * \param targetSteps Number of iterations which are considered "optimal" - */ - void setTargetSteps(int targetSteps) - { targetSteps_ = targetSteps; } - - /*! - * \brief Set the number of iterations after which the Newton - * method gives up. - * \param maxSteps docme - */ - void setMaxSteps(int maxSteps) - { maxSteps_ = maxSteps; } - - /*! - * \brief Returns true if another iteration should be done. - * - * \param uCurrentIter The solution of the current newton iteration - */ - bool newtonProceed(SolutionVector &uCurrentIter) - { - if (numSteps_ < 2) - return true; // we always do at least two iterations - else - if (numSteps_ >= maxSteps_) - { - // we have exceeded the allowed number of steps. if the - // relative error was reduced by a factor of at least 5, - // we proceed anyway - return error_*4.0 < lastError_ && !asImp_().newtonConverged(); - } - else - if (asImp_().newtonConverged()) - return false; // we are below the desired tolerance - - return true; // do another round - } - - /*! - * \brief Returns true if the error of the solution is below the - * tolerance. - */ - bool newtonConverged() const - { return error_ <= tolerance_; } - - /*! - * \brief Called before the newton method is applied to an - * non-linear system of equations. - * - * \param method docme - * \param uCurrentIter docme - * - */ - void newtonBegin(NewtonMethod &method, const SolutionVector &uCurrentIter) - { - method_ = &method; - numSteps_ = 0; - - if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence)) - convergenceWriter_.beginTimestep(); - } - - /*! - * \brief Indidicates the beginning of a newton iteration. - */ - void newtonBeginStep() - { lastError_ = error_; } - - /*! - * \brief Returns the number of steps done since newtonBegin() was - * called. - */ - int newtonNumSteps() - { return numSteps_; } + { }; /*! @@ -292,7 +117,7 @@ public: { // calculate the relative error as the maximum relative // deflection in any degree of freedom. - error_ = 0; + this->error_ = 0; SolutionVector uNewI = uLastIter; uNewI -= deltaU; @@ -302,12 +127,11 @@ public: Scalar vertexError = std::abs(deltaU[i][j]); vertexError /= std::max<Scalar>(1.0, std::abs(uLastIter[i][j] + uNewI[i][j])/2); - error_ = std::max(error_, vertexError); + this->error_ = std::max(this->error_, vertexError); } } } - /*! * \brief Solve the linear system of equations \f$ \mathbf{A}x - b * = 0\f$. @@ -332,14 +156,8 @@ public: try { int converged = linearSolver_.solve(A, x, b); - // make sure all processes converged -/* -* \brief docme -* \param convergedSend docme -* \param converged docme -*/ - #if HAVE_MPI + // make sure all processes converged int convergedSend = 1; MPI_Allreduce(/*sendBuf=*/&convergedSend, /*recvBuf=*/&converged, @@ -354,14 +172,8 @@ public: } } catch (const Dune::MatrixBlockError &e) { - // make sure all processes converged -/* -* \brief docme -* \param convergedSend docme -* \param converged docme -*/ - #if HAVE_MPI + // make sure all processes converged int convergedSend = 0; int converged; @@ -381,15 +193,8 @@ public: throw p; } catch (const Dune::Exception &e) { - // make sure all processes converged - -/* -* \brief docme -* \param convergedSend docme -* \param converged docme -*/ - #if HAVE_MPI + // make sure all processes converged int convergedSend = 0; int converged; @@ -449,36 +254,22 @@ public: { typedef Dumux::SplitAndMerge<TypeTag> Common; - Common::splitSolVector(this->model().curSol(), - this->model().subModel1().curSol(), - this->model().subModel2().curSol()); + Common::splitSolVector(this->model_().curSol(), + this->model_().subModel1().curSol(), + this->model_().subModel2().curSol()); - ++numSteps_; - - if (verbose()) - std::cout << "\rNewton iteration " << numSteps_ << " done: " - << "error=" << error_ << endIterMsg().str() << "\n"; - endIterMsgStream_.str(""); + ParentType::newtonEndStep(uCurrentIter, uLastIter); } - /*! - * \brief Indicates that we're done solving the non-linear system of equations. - */ - void newtonEnd() - { - if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence)) - convergenceWriter_.endTimestep(); - } - - /*! - * \brief Called if the newton method broke down. - * - * This method is called _after_ newtonEnd() - */ - void newtonFail() - { - numSteps_ = std::max(maxSteps_, targetSteps_*2); - } +// /*! +// * \brief Called if the newton method broke down. +// * +// * This method is called _after_ newtonEnd() +// */ +// void newtonFail() +// { +// numSteps_ = std::max(maxSteps_, targetSteps_*2); +// } /*! * \brief Called when the newton method was sucessful. @@ -489,70 +280,6 @@ public: { } - /*! - * \brief Suggest a new time stepsize based on the old time step size. - * - * The default behaviour is to suggest the old time step size - * scaled by the ratio between the target iterations and the - * iterations required to actually solve the last time step. - * - * \param oldTimeStep docme - * - */ - Scalar suggestTimeStepSize(Scalar oldTimeStep) const - { - // be agressive reducing the timestep size but - // conservative when increasing it. the rationale is - // that we want to avoid failing in the next newton - // iteration which would require another linearization - // of the problem. - if (numSteps_ > targetSteps_) { - Scalar percent = ((Scalar) numSteps_ - targetSteps_)/targetSteps_; - return oldTimeStep/(1.0 + percent); - } - else { - /*Scalar percent = (Scalar(1))/targetSteps_; - return oldTimeStep*(1 + percent); - */ - Scalar percent = ((Scalar) targetSteps_ - numSteps_)/targetSteps_; - return oldTimeStep*(1.0 + percent/1.2); - } - } - - /*! - * \brief Returns a reference to the current newton method - * which is controlled by this controller. - */ - NewtonMethod &method() - { return *method_; } - - /*! - * \brief Returns a reference to the current newton method - * which is controlled by this controller. - */ - const NewtonMethod &method() const - { return *method_; } - - /*! - * \brief Returns a reference to the current numeric model. - */ - Model &model() - { return method_->model(); } - - /*! - * \brief Returns a const reference to the current numeric model. - */ - const Model &model() const - { return method_->model(); } - - std::ostringstream &endIterMsg() - { return endIterMsgStream_; } - - const GridView1 &gridView1() const - { return problem_().gridView1(); } - const GridView2 &gridView2() const - { return problem_().gridView2(); } - /*! * \brief the convergence writer produces the output * @@ -566,58 +293,32 @@ public: const SolutionVector &deltaU) { if (GET_PARAM_FROM_GROUP(TypeTag, bool, Newton, WriteConvergence)) { - convergenceWriter_.beginIteration(gridView1_(), gridView2_()); + convergenceWriter_.beginIteration(sdGridView1_(), sdGridView2_()); convergenceWriter_.writeFields(uLastIter, deltaU); convergenceWriter_.endIteration(); }; }; - /*! - * \brief Returns a copy of the the grid view. - */ - const GridView gridView_() const - { return problem_().gridView(); } - /*! * \brief the subdomain gridviews */ - const GridView1 gridView1_() const - { return problem_().subProblem1().gridView(); } - const GridView2 gridView2_() const - { return problem_().subProblem2().gridView(); } - - /*! - * \brief the coupled problem - * - * \returns the actual implementation for the controller we do - * \it this way in order to allow "poor man's virtual methods", - * \i.e. methods of subclasses which can be called by the base - *\class. - */ - Problem &problem_() - { return method_->problem(); } - const Problem &problem_() const - { return method_->problem(); } + const GridView1 sdGridView1_() const + { return this->problem_().sdGridView1(); } + const GridView2 sdGridView2_() const + { return this->problem_().sdGridView2(); } +private: Implementation &asImp_() { return *static_cast<Implementation*>(this); } const Implementation &asImp_() const { return *static_cast<const Implementation*>(this); } - /* - * \brief docme - */ - bool verbose_; std::ostringstream endIterMsgStream_; NewtonMethod *method_; - Scalar tolerance_; - Scalar error_; - Scalar lastError_; - // optimal number of iterations we want to achieve int targetSteps_; // maximum number of iterations we do before giving up diff --git a/dumux/multidomain/common/multidomainproblem.hh b/dumux/multidomain/common/multidomainproblem.hh index 7b48f754ddbbec4e2c41aa2f506b33a86f2c0b82..9b6064b491a0cde14c881e8c6ba7dca53993ad29 100644 --- a/dumux/multidomain/common/multidomainproblem.hh +++ b/dumux/multidomain/common/multidomainproblem.hh @@ -76,8 +76,8 @@ private: typedef typename GET_PROP_TYPE(SubTypeTag1, Problem) SubProblem1; typedef typename GET_PROP_TYPE(SubTypeTag2, Problem) SubProblem2; - typedef typename GET_PROP_TYPE(SubTypeTag1, GridView) GridView1; - typedef typename GET_PROP_TYPE(SubTypeTag2, GridView) GridView2; + typedef typename GET_PROP_TYPE(SubTypeTag1, GridView) SubDomainGridView1; + typedef typename GET_PROP_TYPE(SubTypeTag2, GridView) SubDomainGridView2; typedef typename GET_PROP_TYPE(TypeTag, Grid) HostGrid; typedef typename GET_PROP_TYPE(TypeTag, MultiDomainGrid) MDGrid; @@ -97,26 +97,21 @@ public: * \param timeManager The TimeManager which is used by the simulation * */ - MultiDomainProblem(HostGrid& hostGrid, + MultiDomainProblem(MDGrid &mdGrid, TimeManager &timeManager) - : timeManager_(timeManager), - newtonMethod_(asImp_()), - newtonCtl_(asImp_()) - { - // subdivide grid in two subdomains - subID1_ = 0; - subID2_ = 1; - - mdGrid_ = Dune::make_shared<MDGrid> (hostGrid); - - sdGrid1_ = &(mdGrid_->subDomain(subID1_)); - sdGrid2_ = &(mdGrid_->subDomain(subID2_)); - - subProblem1_ = Dune::make_shared<SubProblem1>(timeManager, sdGrid1_->leafView()); - subProblem2_ = Dune::make_shared<SubProblem2>(timeManager, sdGrid2_->leafView()); - - mdVertexMapper_ = Dune::make_shared<VertexMapper>(mdGrid_->leafView()); - }; + : timeManager_(timeManager) + , newtonMethod_(asImp_()) + , newtonCtl_(asImp_()) + , mdGrid_(mdGrid) + , mdGridView_(mdGrid.leafView()) + , mdVertexMapper_(mdGrid_.leafView()) + , subID1_(0) + , subID2_(1) + , sdGrid1_(mdGrid.subDomain(subID1_)) + , sdGrid2_(mdGrid.subDomain(subID2_)) + , subProblem1_(timeManager, sdGrid1_.leafView()) + , subProblem2_(timeManager, sdGrid2_.leafView()) + { }; /*! * \brief Called by the Dumux::TimeManager in order to @@ -440,27 +435,25 @@ public: * \brief Returns a reference to subproblem1 */ SubProblem1& subProblem1() - { return *subProblem1_; } + { return subProblem1_; } /*! * \brief Returns a const reference to subproblem1 */ - const SubProblem1& subProblem1() const - { return *subProblem1_; } +// const SubProblem1& subProblem1() const +// { return subProblem1_; } /*! * \brief Returns a reference to subproblem2 */ SubProblem2& subProblem2() - { return *subProblem2_; } + { return subProblem2_; } /*! * \brief Returns a const reference to subproblem2 */ - const SubProblem2& subProblem2() const - { return *subProblem2_; } - /////////////////////////////// - +// const SubProblem2& subProblem2() const +// { return subProblem2_; } /*! * \brief Returns a reference to the localresidual1 @@ -478,69 +471,79 @@ public: * \brief Returns a reference to the multidomain grid */ MDGrid& mdGrid() - { return *mdGrid_; } + { return mdGrid_; } /*! * \brief Returns a const reference to the multidomain grid */ const MDGrid& mdGrid() const - { return *mdGrid_; } + { return mdGrid_; } /*! * \brief Returns the multidomain gridview */ - const MDGridView mdGridView() const - { return mdGrid().leafView(); } + const MDGridView& mdGridView() const + { return mdGridView_; } /*! * \brief Returns the multidomain gridview */ - const MDGridView gridView() const - { return mdGrid().leafView(); } + const MDGridView& gridView() const + { return mdGridView_; } /*! - * \brief Returns a reference to the subdomain1 grid + * \brief Provides a vertex mapper for the multidomain + * */ - SDGrid& sdGrid1() - { return *sdGrid1_; } + VertexMapper& mdVertexMapper() + { return mdVertexMapper_; } /*! * \brief Returns a const reference to the subdomain1 grid */ const SDGrid& sdGrid1() const - { return *sdGrid1_; } - - /*! - * \brief Returns a reference to the subdomain2 grid - */ - SDGrid& sdGrid2() - { return *sdGrid2_; } + { return sdGrid1_; } /*! * \brief Returns a const reference to the subdomain2 grid */ const SDGrid& sdGrid2() const - { return *sdGrid2_; } + { return sdGrid2_; } /*! * \brief Returns the gridview of subdomain1 */ - const GridView1 gridView1() const + DUNE_DEPRECATED_MSG("use sdGridView1 instead") + const SubDomainGridView1 gridView1() const { return mdGrid().subDomain(subID1_).leafView(); } /*! * \brief Returns the gridview of subdomain2 */ - const GridView2 gridView2() const + DUNE_DEPRECATED_MSG("use sdGridView2 instead") + const SubDomainGridView2 gridView2() const { return mdGrid().subDomain(subID2_).leafView(); } + /*! + * \brief Returns the gridview of subdomain1 + */ + const SubDomainGridView1 sdGridView1() const + { return sdGrid1_.leafView(); } +// mdGrid().subDomain(subID1_).leafView(); } + + /*! + * \brief Returns the gridview of subdomain2 + */ + const SubDomainGridView2 sdGridView2() const + { return sdGrid2_.leafView(); } + /*! * \brief Returns a pointer to the subdomain1 element * \param mdElement1 docme */ SDElementPointer sdElementPointer1(const MDElement& mdElement1) - { return mdGrid().subDomain(subID1_).subDomainEntityPointer(mdElement1); } + { return sdGrid1_.subDomainEntityPointer(mdElement1); } /*! * \brief Returns a pointer to the subdomain2 element @@ -548,14 +551,7 @@ public: * \param mdElement2 docme */ SDElementPointer sdElementPointer2(const MDElement& mdElement2) - { return mdGrid().subDomain(subID2_).subDomainEntityPointer(mdElement2); } - - /*! - * \brief Provides a vertex mapper for the multidomain - * - */ - VertexMapper& mdVertexMapper() - { return *mdVertexMapper_; } + { return sdGrid2_.subDomainEntityPointer(mdElement2); } protected: @@ -582,21 +578,21 @@ private: TimeManager &timeManager_; NewtonMethod newtonMethod_; NewtonController newtonCtl_; + Model model_; + MDGrid &mdGrid_; + const MDGridView mdGridView_; + VertexMapper mdVertexMapper_; + typename MDGrid::SubDomainType subID1_; typename MDGrid::SubDomainType subID2_; - Dune::shared_ptr<MDGrid> mdGrid_; - Dune::shared_ptr<VertexMapper> mdVertexMapper_; - - SDGrid *sdGrid1_; - SDGrid *sdGrid2_; - - Dune::shared_ptr<SubProblem1> subProblem1_; - Dune::shared_ptr<SubProblem2> subProblem2_; - + const SDGrid& sdGrid1_; + const SDGrid& sdGrid2_; + SubProblem1 subProblem1_; + SubProblem2 subProblem2_; }; // definition of the static class member simname_, diff --git a/dumux/multidomain/common/multidomainproperties.hh b/dumux/multidomain/common/multidomainproperties.hh index 5774aceed54d2288ce3d51f0eccfd78ed1f13e60..8ffc4b7ccfb07396fd1de444d27e481636fc4af0 100644 --- a/dumux/multidomain/common/multidomainproperties.hh +++ b/dumux/multidomain/common/multidomainproperties.hh @@ -24,6 +24,7 @@ #define DUMUX_MULTIDOMAIN_PROPERTIES_HH #include <dumux/implicit/common/implicitproperties.hh> +#include <dumux/nonlinear/newtonmethod.hh> #include <dumux/linear/linearsolverproperties.hh> /*! @@ -47,7 +48,7 @@ namespace Properties ////////////////////////////////////////////////////////////////// //! The type tag for problems which utilize the coupling approach -NEW_TYPE_TAG(MultiDomain, INHERITS_FROM(LinearSolverTypeTag, ImplicitBase)); +NEW_TYPE_TAG(MultiDomain, INHERITS_FROM(NewtonMethod, LinearSolverTypeTag, ImplicitBase)); ////////////////////////////////////////////////////////////////// // Property tags