From 8dac52cd04c29de4093f6f884ea5fe372f81396a Mon Sep 17 00:00:00 2001 From: Klaus Mosthaf <klmos@env.dtu.dk> Date: Fri, 6 Dec 2013 10:39:47 +0000 Subject: [PATCH] multidomainproblem receives the mdGrid now, this allows to use references instead of shared pointers; a copy of the multidomain gridview is now stored in the multidomainproblem to get rid of the "reference to temporary" warning cleanup in multidomainnewtoncontroller, which is now derived from the standard newtoncontroller. The MultiDomain TypeTag is now also derived from NewtonMethod. Reviewed by Markus git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@12140 2fb0f335-1f38-0410-981e-8018bf24f1b0 --- .../common/multidomainassembler.hh | 7 + dumux/multidomain/common/multidomainmodel.hh | 11 + .../common/multidomainnewtoncontroller.hh | 371 ++---------------- .../multidomain/common/multidomainproblem.hh | 128 +++--- .../common/multidomainproperties.hh | 3 +- 5 files changed, 118 insertions(+), 402 deletions(-) diff --git a/dumux/multidomain/common/multidomainassembler.hh b/dumux/multidomain/common/multidomainassembler.hh index f13300fb75..92762907c9 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 db7328b052..c99c7861da 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 2049f0597d..ab34bdf0aa 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 7b48f754dd..9b6064b491 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 5774aceed5..8ffc4b7ccf 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 -- GitLab