diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 49100b41b90a1b8a76dc4418ea4702b1f0be1147..7f057a90507924dbef6cdae956f4d1798f1af337 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -50,6 +50,8 @@ #include #include "newtonconvergencewriter.hh" +#include "primaryvariableswitchadapter.hh" +#include "newtonvariablesbackend.hh" namespace Dumux { namespace Detail { @@ -64,13 +66,6 @@ struct supportsPartialReassembly {} }; -//! helper aliases to extract a primary variable switch from the VolumeVariables (if defined, yields int otherwise) -template -using DetectPVSwitch = typename Assembler::GridVariables::VolumeVariables::PrimaryVariableSwitch; - -template -using GetPVSwitch = Dune::Std::detected_or; - } // end namespace Detail /*! @@ -91,14 +86,18 @@ class NewtonSolver : public PDESolver using ParentType = PDESolver; using Scalar = typename Assembler::Scalar; using JacobianMatrix = typename Assembler::JacobianMatrix; - using SolutionVector = typename Assembler::ResidualType; - using ConvergenceWriter = ConvergenceWriterInterface; using TimeLoop = TimeLoopBase; - using PrimaryVariableSwitch = typename Detail::GetPVSwitch::type; - using HasPriVarsSwitch = typename Detail::GetPVSwitch::value_t; // std::true_type or std::false_type - static constexpr bool hasPriVarsSwitch() { return HasPriVarsSwitch::value; }; + // enable models with primary variable switch + // TODO: Use more general Assembler::Variables once it is available + using PrimaryVariableSwitchAdapter = Dumux::PrimaryVariableSwitchAdapter; + // the variable backend + using Variables = typename Assembler::Variables; + using Backend = NewtonVariablesBackend; + + using SolutionVector = typename Backend::DofVector; + using ConvergenceWriter = ConvergenceWriterInterface; public: using Communication = Comm; @@ -114,6 +113,7 @@ public: , endIterMsgStream_(std::ostringstream::out) , comm_(comm) , paramGroup_(paramGroup) + , priVarSwitchAdapter_(std::make_unique(paramGroup)) { initParams_(paramGroup); @@ -127,12 +127,6 @@ public: // initialize the partial reassembler if (enablePartialReassembly_) partialReassembler_ = std::make_unique(this->assembler()); - - if (hasPriVarsSwitch()) - { - const int priVarSwitchVerbosity = getParamFromGroup(paramGroup, "PrimaryVariableSwitch.Verbosity", 1); - priVarSwitch_ = std::make_unique(priVarSwitchVerbosity); - } } //! the communicator for parallel runs @@ -263,7 +257,7 @@ public: virtual void newtonBegin(SolutionVector& u) { numSteps_ = 0; - initPriVarSwitch_(u, HasPriVarsSwitch{}); + priVarSwitchAdapter_->initialize(u, this->assembler().gridVariables()); // write the initial residual if a convergence writer was set if (convergenceWriter_) @@ -274,11 +268,11 @@ public: convergenceWriter_->write(u, delta, this->assembler().residual()); } - if (enablePartialReassembly_) - { - partialReassembler_->resetColors(); - resizeDistanceFromLastLinearization_(u, distanceFromLastLinearization_); - } + // if (enablePartialReassembly_) + // { + // partialReassembler_->resetColors(); + // resizeDistanceFromLastLinearization_(u, distanceFromLastLinearization_); + // } } /*! @@ -354,14 +348,7 @@ public: try { if (numSteps_ == 0) - { - Scalar norm2 = b.two_norm2(); - if (comm_.size() > 1) - norm2 = comm_.sum(norm2); - - using std::sqrt; - initialResidual_ = sqrt(norm2); - } + initialResidualNorm_ = this->linearSolver().norm(b); // solve by calling the appropriate implementation depending on whether the linear solver // is capable of handling MultiType matrices or not @@ -418,40 +405,40 @@ public: if (enableShiftCriterion_ || enablePartialReassembly_) newtonUpdateShift_(uLastIter, deltaU); - if (enablePartialReassembly_) { - // Determine the threshold 'eps' that is used for the partial reassembly. - // Every entity where the primary variables exhibit a relative shift - // summed up since the last linearization above 'eps' will be colored - // red yielding a reassembly. - // The user can provide three parameters to influence the threshold: - // 'minEps' by 'Newton.ReassemblyMinThreshold' (1e-1*shiftTolerance_ default) - // 'maxEps' by 'Newton.ReassemblyMaxThreshold' (1e2*shiftTolerance_ default) - // 'omega' by 'Newton.ReassemblyShiftWeight' (1e-3 default) - // The threshold is calculated from the currently achieved maximum - // relative shift according to the formula - // eps = max( minEps, min(maxEps, omega*shift) ). - // Increasing/decreasing 'minEps' leads to less/more reassembly if - // 'omega*shift' is small, i.e., for the last Newton iterations. - // Increasing/decreasing 'maxEps' leads to less/more reassembly if - // 'omega*shift' is large, i.e., for the first Newton iterations. - // Increasing/decreasing 'omega' leads to more/less first and last - // iterations in this sense. - using std::max; - using std::min; - auto reassemblyThreshold = max(reassemblyMinThreshold_, - min(reassemblyMaxThreshold_, - shift_*reassemblyShiftWeight_)); - - updateDistanceFromLastLinearization_(uLastIter, deltaU); - partialReassembler_->computeColors(this->assembler(), - distanceFromLastLinearization_, - reassemblyThreshold); - - // set the discrepancy of the red entities to zero - for (unsigned int i = 0; i < distanceFromLastLinearization_.size(); i++) - if (partialReassembler_->dofColor(i) == EntityColor::red) - distanceFromLastLinearization_[i] = 0; - } + // if (enablePartialReassembly_) { + // // Determine the threshold 'eps' that is used for the partial reassembly. + // // Every entity where the primary variables exhibit a relative shift + // // summed up since the last linearization above 'eps' will be colored + // // red yielding a reassembly. + // // The user can provide three parameters to influence the threshold: + // // 'minEps' by 'Newton.ReassemblyMinThreshold' (1e-1*shiftTolerance_ default) + // // 'maxEps' by 'Newton.ReassemblyMaxThreshold' (1e2*shiftTolerance_ default) + // // 'omega' by 'Newton.ReassemblyShiftWeight' (1e-3 default) + // // The threshold is calculated from the currently achieved maximum + // // relative shift according to the formula + // // eps = max( minEps, min(maxEps, omega*shift) ). + // // Increasing/decreasing 'minEps' leads to less/more reassembly if + // // 'omega*shift' is small, i.e., for the last Newton iterations. + // // Increasing/decreasing 'maxEps' leads to less/more reassembly if + // // 'omega*shift' is large, i.e., for the first Newton iterations. + // // Increasing/decreasing 'omega' leads to more/less first and last + // // iterations in this sense. + // using std::max; + // using std::min; + // auto reassemblyThreshold = max(reassemblyMinThreshold_, + // min(reassemblyMaxThreshold_, + // shift_*reassemblyShiftWeight_)); + // + // updateDistanceFromLastLinearization_(uLastIter, deltaU); + // partialReassembler_->computeColors(this->assembler(), + // distanceFromLastLinearization_, + // reassemblyThreshold); + // + // // set the discrepancy of the red entities to zero + // for (unsigned int i = 0; i < distanceFromLastLinearization_.size(); i++) + // if (partialReassembler_->dofColor(i) == EntityColor::red) + // distanceFromLastLinearization_[i] = 0; + // } if (useLineSearch_) lineSearchUpdate_(uCurrentIter, uLastIter, deltaU); @@ -486,7 +473,7 @@ public: virtual void newtonEndStep(SolutionVector &uCurrentIter, const SolutionVector &uLastIter) { - invokePriVarSwitch_(uCurrentIter, HasPriVarsSwitch{}); + priVarSwitchAdapter_->invoke(uCurrentIter, this->assembler().gridVariables()); ++numSteps_; @@ -535,7 +522,7 @@ public: { // in case the model has a priVar switch and some some primary variables // actually switched their state in the last iteration, enforce another iteration - if (priVarsSwitchedInLastIteration_) + if (priVarSwitchAdapter_->switched()) return false; if (enableShiftCriterion_ && !enableResidualCriterion_) @@ -706,66 +693,16 @@ protected: void computeResidualReduction_(const SolutionVector &uCurrentIter) { - residualNorm_ = this->assembler().residualNorm(uCurrentIter); + auto residual = Backend::makeZeroDofVector(Backend::size(uCurrentIter)); + this->assembler().assembleResidual(residual, uCurrentIter); + residualNorm_ = this->linearSolver().norm(residual); reduction_ = residualNorm_; - reduction_ /= initialResidual_; + reduction_ /= initialResidualNorm_; } bool enableResidualCriterion() const { return enableResidualCriterion_; } - /*! - * \brief Initialize the privar switch, noop if there is no priVarSwitch - */ - void initPriVarSwitch_(SolutionVector&, std::false_type) {} - - /*! - * \brief Initialize the privar switch - */ - void initPriVarSwitch_(SolutionVector& sol, std::true_type) - { - priVarSwitch_->reset(sol.size()); - priVarsSwitchedInLastIteration_ = false; - - const auto& problem = this->assembler().problem(); - const auto& gridGeometry = this->assembler().gridGeometry(); - auto& gridVariables = this->assembler().gridVariables(); - priVarSwitch_->updateBoundary(problem, gridGeometry, gridVariables, sol); - } - - /*! - * \brief Switch primary variables if necessary, noop if there is no priVarSwitch - */ - void invokePriVarSwitch_(SolutionVector&, std::false_type) {} - - /*! - * \brief Switch primary variables if necessary - */ - void invokePriVarSwitch_(SolutionVector& uCurrentIter, std::true_type) - { - // update the variable switch (returns true if the pri vars at at least one dof were switched) - // for disabled grid variable caching - const auto& gridGeometry = this->assembler().gridGeometry(); - const auto& problem = this->assembler().problem(); - auto& gridVariables = this->assembler().gridVariables(); - - // invoke the primary variable switch - priVarsSwitchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, gridVariables, - problem, gridGeometry); - - if (priVarsSwitchedInLastIteration_) - { - for (const auto& element : elements(gridGeometry.gridView())) - { - // if the volume variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, gridVariables, uCurrentIter); - - // if the flux variables are cached globally, we need to update those where the primary variables have been switched - priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, gridVariables, uCurrentIter); - } - } - } - //! optimal number of iterations we want to achieve int targetSteps_; //! minimum number of iterations we do @@ -779,7 +716,7 @@ protected: Scalar reduction_; Scalar residualNorm_; Scalar lastReduction_; - Scalar initialResidual_; + Scalar initialResidualNorm_; // shift criterion variables Scalar shift_; @@ -803,8 +740,8 @@ private: newtonBegin(uCurrentIter); // the given solution is the initial guess - SolutionVector uLastIter(uCurrentIter); - SolutionVector deltaU(uCurrentIter); + auto uLastIter = Backend::makeDofVector(uCurrentIter); + auto deltaU = Backend::makeDofVector(uCurrentIter); // setup timers Dune::Timer assembleTimer(false); @@ -951,43 +888,41 @@ private: virtual void newtonUpdateShift_(const SolutionVector &uLastIter, const SolutionVector &deltaU) { - shift_ = 0; - newtonUpdateShiftImpl_(uLastIter, deltaU); - + shift_ = Backend::maxRelativeShift(uLastIter, deltaU); if (comm_.size() > 1) shift_ = comm_.max(shift_); } - template - void newtonUpdateShiftImpl_(const SolVec &uLastIter, - const SolVec &deltaU) - { - for (int i = 0; i < int(uLastIter.size()); ++i) - { - auto uNewI = uLastIter[i]; - uNewI -= deltaU[i]; - - Scalar shiftAtDof = relativeShiftAtDof_(uLastIter[i], uNewI); - using std::max; - shift_ = max(shift_, shiftAtDof); - } - } - - template - void newtonUpdateShiftImpl_(const Dune::MultiTypeBlockVector &uLastIter, - const Dune::MultiTypeBlockVector &deltaU) - { - // There seems to be a bug in g++5 which which prevents compilation when - // passing the call to the implementation directly to Dune::Hybrid::forEach. - // We therefore store this call in a lambda and pass it to the for loop afterwards. - auto doUpdate = [&](const auto subVectorIdx) - { - this->newtonUpdateShiftImpl_(uLastIter[subVectorIdx], deltaU[subVectorIdx]); - }; - - using namespace Dune::Hybrid; - forEach(integralRange(Dune::Hybrid::size(uLastIter)), doUpdate); - } + // template + // void newtonUpdateShiftImpl_(const SolVec &uLastIter, + // const SolVec &deltaU) + // { + // for (int i = 0; i < int(uLastIter.size()); ++i) + // { + // auto uNewI = uLastIter[i]; + // uNewI -= deltaU[i]; + // + // Scalar shiftAtDof = relativeShiftAtDof_(uLastIter[i], uNewI); + // using std::max; + // shift_ = max(shift_, shiftAtDof); + // } + // } + + // template + // void newtonUpdateShiftImpl_(const Dune::MultiTypeBlockVector &uLastIter, + // const Dune::MultiTypeBlockVector &deltaU) + // { + // // There seems to be a bug in g++5 which which prevents compilation when + // // passing the call to the implementation directly to Dune::Hybrid::forEach. + // // We therefore store this call in a lambda and pass it to the for loop afterwards. + // auto doUpdate = [&](const auto subVectorIdx) + // { + // this->newtonUpdateShiftImpl_(uLastIter[subVectorIdx], deltaU[subVectorIdx]); + // }; + // + // using namespace Dune::Hybrid; + // forEach(integralRange(Dune::Hybrid::size(uLastIter)), doUpdate); + // } virtual void lineSearchUpdate_(SolutionVector &uCurrentIter, const SolutionVector &uLastIter, @@ -1046,25 +981,30 @@ private: SolutionVector& x, SolutionVector& b) { - //! Copy into a standard block vector. - //! This is necessary for all model _not_ using a FieldVector as - //! primary variables vector in combination with UMFPack or SuperLU as their interfaces are hard coded - //! to this field vector type in Dune ISTL - //! Could be avoided for vectors that already have the right type using SFINAE - //! but it shouldn't impact performance too much - constexpr auto blockSize = std::decay_t::dimension; - using BlockType = Dune::FieldVector; - Dune::BlockVector xTmp; xTmp.resize(b.size()); - Dune::BlockVector bTmp(xTmp); - for (unsigned int i = 0; i < b.size(); ++i) - for (unsigned int j = 0; j < blockSize; ++j) - bTmp[i][j] = b[i][j]; + // //! Copy into a standard block vector. + // //! This is necessary for all model _not_ using a FieldVector as + // //! primary variables vector in combination with UMFPack or SuperLU as their interfaces are hard coded + // //! to this field vector type in Dune ISTL + // //! Could be avoided for vectors that already have the right type using SFINAE + // //! but it shouldn't impact performance too much + // constexpr auto blockSize = std::decay_t::dimension; + // using BlockType = Dune::FieldVector; + // Dune::BlockVector xTmp; xTmp.resize(b.size()); + // Dune::BlockVector bTmp(xTmp); + // for (unsigned int i = 0; i < b.size(); ++i) + // for (unsigned int j = 0; j < blockSize; ++j) + // bTmp[i][j] = b[i][j]; + + auto xTmp = Backend::makeZeroDofVectorForSolver(Backend::size(b)); + auto bTmp = Backend::makeDofVectorForSolver(b); const int converged = ls.solve(A, xTmp, bTmp); - for (unsigned int i = 0; i < x.size(); ++i) - for (unsigned int j = 0; j < blockSize; ++j) - x[i][j] = xTmp[i][j]; + x = Backend::reconstructDofVectorFromSolver(xTmp); + + // for (unsigned int i = 0; i < x.size(); ++i) + // for (unsigned int j = 0; j < blockSize; ++j) + // x[i][j] = xTmp[i][j]; return converged; } @@ -1181,63 +1121,63 @@ private: reportParams(); } - template - void updateDistanceFromLastLinearization_(const Sol& u, const Sol& uDelta) - { - for (size_t i = 0; i < u.size(); ++i) { - const auto& currentPriVars(u[i]); - auto nextPriVars(currentPriVars); - nextPriVars -= uDelta[i]; - - // add the current relative shift for this degree of freedom - auto shift = relativeShiftAtDof_(currentPriVars, nextPriVars); - distanceFromLastLinearization_[i] += shift; - } - } - - template - void updateDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector& uLastIter, - const Dune::MultiTypeBlockVector& deltaU) - { - DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector"); - } - - template - void resizeDistanceFromLastLinearization_(const Sol& u, std::vector& dist) - { - dist.assign(u.size(), 0.0); - } - - template - void resizeDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector& u, - std::vector& dist) - { - DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector"); - } - - /*! - * \brief Returns the maximum relative shift between two vectors of - * primary variables. - * - * \param priVars1 The first vector of primary variables - * \param priVars2 The second vector of primary variables - */ - template - Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1, - const PrimaryVariables &priVars2) const - { - Scalar result = 0.0; - using std::abs; - using std::max; - // iterate over all primary variables - for (int j = 0; j < PrimaryVariables::dimension; ++j) { - Scalar eqErr = abs(priVars1[j] - priVars2[j]); - eqErr /= max(1.0,abs(priVars1[j] + priVars2[j])/2); - - result = max(result, eqErr); - } - return result; - } + // template + // void updateDistanceFromLastLinearization_(const Sol& u, const Sol& uDelta) + // { + // for (size_t i = 0; i < u.size(); ++i) { + // const auto& currentPriVars(u[i]); + // auto nextPriVars(currentPriVars); + // nextPriVars -= uDelta[i]; + // + // // add the current relative shift for this degree of freedom + // auto shift = relativeShiftAtDof_(currentPriVars, nextPriVars); + // distanceFromLastLinearization_[i] += shift; + // } + // } + // + // template + // void updateDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector& uLastIter, + // const Dune::MultiTypeBlockVector& deltaU) + // { + // DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector"); + // } + + // template + // void resizeDistanceFromLastLinearization_(const Sol& u, std::vector& dist) + // { + // dist.assign(u.size(), 0.0); + // } + // + // template + // void resizeDistanceFromLastLinearization_(const Dune::MultiTypeBlockVector& u, + // std::vector& dist) + // { + // DUNE_THROW(Dune::NotImplemented, "Reassembly for MultiTypeBlockVector"); + // } + + // /*! + // * \brief Returns the maximum relative shift between two vectors of + // * primary variables. + // * + // * \param priVars1 The first vector of primary variables + // * \param priVars2 The second vector of primary variables + // */ + // template + // Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1, + // const PrimaryVariables &priVars2) const + // { + // Scalar result = 0.0; + // using std::abs; + // using std::max; + // // iterate over all primary variables + // for (int j = 0; j < PrimaryVariables::dimension; ++j) { + // Scalar eqErr = abs(priVars1[j] - priVars2[j]); + // eqErr /= max(1.0,abs(priVars1[j] + priVars2[j])/2); + // + // result = max(result, eqErr); + // } + // return result; + // } //! The communication object Communication comm_; @@ -1280,9 +1220,7 @@ private: std::size_t numLinearSolverBreakdowns_ = 0; //! total number of linear solves that failed //! the class handling the primary variable switch - std::unique_ptr priVarSwitch_; - //! if we switched primary variables in the last iteration - bool priVarsSwitchedInLastIteration_ = false; + std::unique_ptr priVarSwitchAdapter_; //! convergence writer std::shared_ptr convergenceWriter_ = nullptr; diff --git a/dumux/nonlinear/newtonvariablesbackend.hh b/dumux/nonlinear/newtonvariablesbackend.hh new file mode 100644 index 0000000000000000000000000000000000000000..a379eb8ec127846b72ef08ebfb90778ef3fb860a --- /dev/null +++ b/dumux/nonlinear/newtonvariablesbackend.hh @@ -0,0 +1,86 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/**************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Nonlinear + * \brief Backends for the Newton for different variable types + */ +#ifndef DUMUX_NEWTON_VARIABLES_BACKEND_HH +#define DUMUX_NEWTON_VARIABLES_BACKEND_HH + +#include +#include + +namespace Dumux { + +/*! + * \file + * \ingroup Nonlinear + * \brief Class providing operations with variables needed in the Newton solver + */ +template +class NewtonVariablesBackend; + +/*! + * \file + * \ingroup Nonlinear + * \brief Class providing Newton operations for scalar/number types + */ +template +class NewtonVariablesBackend::value, Scalar>> +{ +public: + using Variables = Scalar; //!< the type of the variables object + using DofVector = Scalar; //!< the type of the dofs parametrizing the variables object + + static std::size_t size(const DofVector& d) + { return 1; } + + static DofVector makeDofVector(const DofVector& d) + { return d; } + + static DofVector makeZeroDofVector(std::size_t size) + { return 0.0; } + + static Scalar makeDofVectorForSolver(const DofVector& d) + { return d; } + + static Scalar makeZeroDofVectorForSolver(std::size_t size) + { return 0.0; } + + static Scalar reconstructDofVectorFromSolver(const Scalar& d) + { return d; } + + static Scalar maxRelativeShift(const DofVector& previous, const DofVector& delta) + { + const auto current = previous - delta; + using std::abs; using std::max; + auto error = abs(previous - current); + error /= max(1.0, abs(previous + current)/2.0); + return error; + } + + // //! operations on variables + // static DofVector& getDofVector(Variables& v) + // { return v; } +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/nonlinear/primaryvariableswitchadapter.hh b/dumux/nonlinear/primaryvariableswitchadapter.hh new file mode 100644 index 0000000000000000000000000000000000000000..693ecc78a70554bf99383d9deebb1f3be5058874 --- /dev/null +++ b/dumux/nonlinear/primaryvariableswitchadapter.hh @@ -0,0 +1,137 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/**************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Nonlinear + * \brief An adapter for the Newton to manage models with primary variable switch + */ +#ifndef DUMUX_NONLINEAR_PRIMARY_VARIABLE_SWITCH_ADAPTER_HH +#define DUMUX_NONLINEAR_PRIMARY_VARIABLE_SWITCH_ADAPTER_HH + +#include +#include + +namespace Dumux { +namespace Impl { + +//! helper aliases to extract a primary variable switch from the VolumeVariables (if defined, yields int otherwise) +template +using DetectPVSwitch = typename Variables::VolumeVariables::PrimaryVariableSwitch; + +template +using GetPVSwitch = Dune::Std::detected_or; + +template +using PrimaryVariableSwitch = typename GetPVSwitch::type; + +template +constexpr bool hasPriVarsSwitch() { return typename GetPVSwitch::value_t(); }; + +} // end namespace Impl + +/*! + * \ingroup Nonlinear + * \brief An adapter for the Newton to manage models with primary variable switch + */ +template ()> +class PrimaryVariableSwitchAdapter +{ + using PrimaryVariableSwitch = typename Impl::PrimaryVariableSwitch; + +public: + PrimaryVariableSwitchAdapter(const std::string& paramGroup = "") + { + const int priVarSwitchVerbosity = getParamFromGroup(paramGroup, "PrimaryVariableSwitch.Verbosity", 1); + priVarSwitch_ = std::make_unique(priVarSwitchVerbosity); + } + + /*! + * \brief Initialize the privar switch + */ + template + void initialize(SolutionVector& sol, Variables& vars) + { + priVarSwitch_->reset(sol.size()); + priVarsSwitchedInLastIteration_ = false; + const auto& problem = vars.curGridVolVars().problem(); + const auto& gridGeometry = problem.gridGeometry(); + priVarSwitch_->updateBoundary(problem, gridGeometry, vars, sol); + } + + /*! + * \brief Switch primary variables if necessary + */ + template + void invoke(SolutionVector& uCurrentIter, Variables& vars) + { + // update the variable switch (returns true if the pri vars at at least one dof were switched) + // for disabled grid variable caching + const auto& problem = vars.curGridVolVars().problem(); + const auto& gridGeometry = problem.gridGeometry(); + + // invoke the primary variable switch + priVarsSwitchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, vars, problem, gridGeometry); + if (priVarsSwitchedInLastIteration_) + { + for (const auto& element : elements(gridGeometry.gridView())) + { + // if the volume variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedVolVars(problem, element, gridGeometry, vars, uCurrentIter); + + // if the flux variables are cached globally, we need to update those where the primary variables have been switched + priVarSwitch_->updateSwitchedFluxVarsCache(problem, element, gridGeometry, vars, uCurrentIter); + } + } + } + + /*! + * \brief Whether the primary variables have been switched in the last call to invoke + */ + bool switched() const + { return priVarsSwitchedInLastIteration_; } + +private: + //! the class handling the primary variable switch + std::unique_ptr priVarSwitch_; + //! if we switched primary variables in the last iteration + bool priVarsSwitchedInLastIteration_ = false; +}; + +/*! + * \ingroup Nonlinear + * \brief An empty adapter for the Newton for models without primary variable switch + */ +template +class PrimaryVariableSwitchAdapter +{ +public: + PrimaryVariableSwitchAdapter(const std::string& paramGroup = "") {} + + template + void initialize(SolutionVector&, Variables&) {} + + template + void invoke(SolutionVector&, Variables&) {} + + bool switched() const { return false; } +}; + +} // end namespace Dumux + +#endif diff --git a/test/nonlinear/newton/test_newton.cc b/test/nonlinear/newton/test_newton.cc index 5f37e2c3656bbdaa8f0574888023e0823139cd6b..24fdc56e90f51b11146265baabb4ca3b624c41b1 100644 --- a/test/nonlinear/newton/test_newton.cc +++ b/test/nonlinear/newton/test_newton.cc @@ -27,40 +27,11 @@ namespace Dumux { -class MockScalarVariables -{ -public: - static constexpr int dimension = 1; - MockScalarVariables() : var_(0.0) {} - explicit MockScalarVariables(double v) : var_(v) {} - MockScalarVariables& operator-=(const MockScalarVariables& other) { var_ -= other.var_; return *this; } - double& operator[] (int i) { return var_; } - const double& operator[] (int i) const { return var_; } -private: - double var_; -}; - -class MockScalarResidual -{ -public: - MockScalarResidual() : res_(0.0) {} - explicit MockScalarResidual(double r) : res_(r) {} - MockScalarResidual& operator=(double r) { res_[0] = r; return *this; } - MockScalarResidual& operator*=(double a) { res_[0] *= a; return *this; } - MockScalarResidual& operator+=(const MockScalarResidual& other) { res_[0] += other.res_[0]; return *this; } - MockScalarResidual& operator-=(const MockScalarResidual& other) { res_[0] -= other.res_[0]; return *this; } - constexpr std::size_t size() const { return 1; } - MockScalarVariables& operator[] (int i) { return res_; } - const MockScalarVariables& operator[] (int i) const { return res_; } - double two_norm2() const { return res_[0]*res_[0]; } -private: - MockScalarVariables res_; -}; - class MockScalarAssembler { public: - using ResidualType = MockScalarResidual; + using Variables = double; + using ResidualType = double; using Scalar = double; using JacobianMatrix = double; @@ -70,34 +41,36 @@ public: ResidualType prevSol() { return ResidualType(0.0); } - void resetTimeStep(const ResidualType& sol) {} + void resetTimeStep(const ResidualType& x) {} - void assembleResidual(const ResidualType& sol) + void assembleResidual(const ResidualType& x) { - res_[0][0] = sol[0][0]*sol[0][0] - 5.0; + res_ = x*x - 5.0; } - void assembleJacobianAndResidual (const ResidualType& sol) + void assembleResidual(ResidualType& r, const ResidualType& x) const { - assembleResidual(sol); - jac_ = 2.0*sol[0][0]; + r = x*x - 5.0; + } + + void assembleJacobianAndResidual (const ResidualType& x) + { + assembleResidual(x); + jac_ = 2.0*x; } JacobianMatrix& jacobian() { return jac_; } ResidualType& residual() { return res_; } - double residualNorm(const ResidualType& sol) - { - assembleResidual(sol); - return res_[0][0]; - } - void updateGridVariables(const ResidualType& sol) {} + Variables& gridVariables() { return vars_; } + private: JacobianMatrix jac_; ResidualType res_; + Variables vars_; }; class MockScalarLinearSolver @@ -105,12 +78,11 @@ class MockScalarLinearSolver public: void setResidualReduction(double residualReduction) {} - template - bool solve (const double& A, Vector& x, const Vector& b) const - { - x[0][0] = b[0][0]/A; - return true; - } + bool solve (const double& A, double& x, const double& b) const + { x = b/A; return true; } + + double norm (const double& r) + { return std::abs(r); } }; } // end namespace Dumux @@ -137,16 +109,16 @@ int main(int argc, char* argv[]) try auto solver = std::make_shared(assembler, linearSolver); double initialGuess = 0.1; - MockScalarResidual x(initialGuess); + double x(initialGuess); std::cout << "Solving: x^2 - 5 = 0" << std::endl; solver->solve(x); - std::cout << "Solution: " << std::setprecision(15) << x[0][0] + std::cout << "Solution: " << std::setprecision(15) << x << ", exact: " << std::sqrt(5.0) - << ", error: " << std::abs(x[0][0]-std::sqrt(5.0))/std::sqrt(5.0)*100 << "%" << std::endl; + << ", error: " << std::abs(x-std::sqrt(5.0))/std::sqrt(5.0)*100 << "%" << std::endl; - if (Dune::FloatCmp::ne(x[0][0], std::sqrt(5.0), 1e-13)) - DUNE_THROW(Dune::Exception, "Didn't find correct root: " << std::setprecision(15) << x[0][0] << ", exact: " << std::sqrt(5.0)); + if (Dune::FloatCmp::ne(x, std::sqrt(5.0), 1e-13)) + DUNE_THROW(Dune::Exception, "Didn't find correct root: " << std::setprecision(15) << x << ", exact: " << std::sqrt(5.0)); return 0;