From d99d5076b415a1d2bd517d7330ba15a11fdec9ea Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Mon, 9 Apr 2018 07:26:20 +0200
Subject: [PATCH] Delete obsolete NewtonController and NewtonMethod headers
---
dumux/nonlinear/CMakeLists.txt | 3 +-
dumux/nonlinear/newtoncontroller.hh | 651 ------------------
dumux/nonlinear/newtonmethod.hh | 222 ------
dumux/nonlinear/staggerednewtoncontroller.hh | 299 --------
dumux/porousmediumflow/2pnc/CMakeLists.txt | 1 -
.../privarswitchnewtoncontroller.hh | 220 ------
.../nonequilibrium/CMakeLists.txt | 2 +-
.../nonequilibrium/newtoncontroller.hh | 66 --
.../porousmediumflow/richards/CMakeLists.txt | 3 +-
.../richards/newtoncontroller.hh | 122 ----
dumux/porousmediumflow/richardsnc/model.hh | 1 -
11 files changed, 4 insertions(+), 1586 deletions(-)
delete mode 100644 dumux/nonlinear/newtoncontroller.hh
delete mode 100644 dumux/nonlinear/newtonmethod.hh
delete mode 100644 dumux/nonlinear/staggerednewtoncontroller.hh
delete mode 100644 dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh
delete mode 100644 dumux/porousmediumflow/nonequilibrium/newtoncontroller.hh
delete mode 100644 dumux/porousmediumflow/richards/newtoncontroller.hh
diff --git a/dumux/nonlinear/CMakeLists.txt b/dumux/nonlinear/CMakeLists.txt
index 25e858ce85..0bc9109ce1 100644
--- a/dumux/nonlinear/CMakeLists.txt
+++ b/dumux/nonlinear/CMakeLists.txt
@@ -1,7 +1,6 @@
#install headers
install(FILES
-newtoncontroller.hh
newtonconvergencewriter.hh
-newtonmethod.hh
+newtonsolver.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/nonlinear)
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
deleted file mode 100644
index a64002a69d..0000000000
--- a/dumux/nonlinear/newtoncontroller.hh
+++ /dev/null
@@ -1,651 +0,0 @@
-// -*- 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 2 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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Nonlinear
- * \brief Reference implementation of a controller class for the Newton solver.
- *
- * Usually this controller should be sufficient.
- */
-#ifndef DUMUX_NEWTON_CONTROLLER_HH
-#define DUMUX_NEWTON_CONTROLLER_HH
-
-#include <cmath>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/parallel/mpicollectivecommunication.hh>
-#include <dune/common/parallel/mpihelper.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dumux/common/exceptions.hh>
-#include <dumux/common/timeloop.hh>
-#include <dune/common/deprecated.hh>
-
-#warning "This file is deprecated. Use NewtonSolver instead."
-
-namespace Dumux {
-
-/*!
- * \ingroup Nonlinear
- * \brief An implementation of a Newton controller
- * \tparam Scalar the scalar type
- * \tparam Comm the communication object used to communicate with all processes
- * \note If you want to specialize only some methods but are happy with the
- * defaults of the reference controller, derive your controller from
- * this class and simply overload the required methods.
- */
-template <class Scalar,
- class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
-class DUNE_DEPRECATED_MSG("Use NewtonSolver instead.") NewtonController
-{
-
-public:
- //! the communication type used to communicate with all processes
- using Communication = Comm;
-
- /*!
- * \brief Constructor for stationary problems
- */
- NewtonController(const Communication& comm = Dune::MPIHelper::getCollectiveCommunication(),
- const std::string& paramGroup = "")
- : comm_(comm)
- , endIterMsgStream_(std::ostringstream::out)
- , paramGroup_(paramGroup)
- {
- initParams_(paramGroup);
- }
-
- /*!
- * \brief Constructor for instationary problems
- */
- NewtonController(std::shared_ptr<TimeLoop<Scalar>> timeLoop,
- const Communication& comm = Dune::MPIHelper::getCollectiveCommunication(),
- const std::string& paramGroup = "")
- : comm_(comm)
- , timeLoop_(timeLoop)
- , endIterMsgStream_(std::ostringstream::out)
- , paramGroup_(paramGroup)
- {
- initParams_(paramGroup);
- }
-
- //! the communicator for parallel runs
- const Communication& comm() const
- { return comm_; }
-
- /*!
- * \brief Set the maximum acceptable difference of any primary variable
- * between two iterations for declaring convergence.
- *
- * \param tolerance The maximum relative shift between two Newton
- * iterations at which the scheme is considered finished
- */
- void setMaxRelativeShift(Scalar tolerance)
- { shiftTolerance_ = tolerance; }
-
- /*!
- * \brief Set the maximum acceptable absolute residual for declaring convergence.
- *
- * \param tolerance The maximum absolute residual at which
- * the scheme is considered finished
- */
- void setMaxAbsoluteResidual(Scalar tolerance)
- { residualTolerance_ = 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 setResidualReduction(Scalar tolerance)
- { reductionTolerance_ = 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 Number of iterations after we give up
- */
- void setMaxSteps(int maxSteps)
- { maxSteps_ = maxSteps; }
-
- /*!
- * \brief Returns true if another iteration should be done.
- *
- * \param uCurrentIter The solution of the current Newton iteration
- * \param converged if the Newton method's convergence criterion was met in this step
- */
- template<class SolutionVector>
- bool newtonProceed(const SolutionVector &uCurrentIter, bool converged)
- {
- if (numSteps_ < 2)
- return true; // we always do at least two iterations
- else if (converged) {
- return false; // we are below the desired tolerance
- }
- else if (numSteps_ >= maxSteps_) {
- // We have exceeded the allowed number of steps. If the
- // maximum relative shift was reduced by a factor of at least 4,
- // we proceed even if we are above the maximum number of steps.
- if (enableShiftCriterion_)
- return shift_*4.0 < lastShift_;
- else
- return reduction_*4.0 < lastReduction_;
- }
-
- return true;
- }
-
- /*!
- * \brief Returns true if the error of the solution is below the
- * tolerance.
- */
- bool newtonConverged() const
- {
- if (enableShiftCriterion_ && !enableResidualCriterion_)
- {
- return shift_ <= shiftTolerance_;
- }
- else if (!enableShiftCriterion_ && enableResidualCriterion_)
- {
- if(enableAbsoluteResidualCriterion_)
- return residualNorm_ <= residualTolerance_;
- else
- return reduction_ <= reductionTolerance_;
- }
- else if (satisfyResidualAndShiftCriterion_)
- {
- if(enableAbsoluteResidualCriterion_)
- return shift_ <= shiftTolerance_
- && residualNorm_ <= residualTolerance_;
- else
- return shift_ <= shiftTolerance_
- && reduction_ <= reductionTolerance_;
- }
- else
- {
- return shift_ <= shiftTolerance_
- || reduction_ <= reductionTolerance_
- || residualNorm_ <= residualTolerance_;
- }
-
- return false;
- }
-
- /*!
- * \brief Called before the Newton method is applied to an
- * non-linear system of equations.
- *
- * \param u The initial solution
- */
- template<class SolutionVector>
- void newtonBegin(const SolutionVector& u)
- {
- numSteps_ = 0;
- }
-
- /*!
- * \brief Indicates the beginning of a Newton iteration.
- */
- template<class SolutionVector>
- void newtonBeginStep(const SolutionVector& u)
- {
- lastShift_ = shift_;
- if (numSteps_ == 0)
- {
- lastReduction_ = 1.0;
- }
- else
- {
- lastReduction_ = reduction_;
- }
- }
-
- /*!
- * \brief Returns the number of steps done since newtonBegin() was
- * called.
- */
- int newtonNumSteps()
- { return numSteps_; }
-
- /*!
- * \brief Update the maximum relative shift of the solution compared to
- * the previous iteration.
- *
- * \param uLastIter The current iterative solution
- * \param deltaU The difference between the current and the next solution
- */
- template<class SolutionVector>
- void newtonUpdateShift(const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- shift_ = 0;
-
- for (int i = 0; i < int(uLastIter.size()); ++i) {
- typename SolutionVector::block_type uNewI = uLastIter[i];
- uNewI -= deltaU[i];
-
- Scalar shiftAtDof = relativeShiftAtDof_(uLastIter[i], uNewI);
- using std::max;
- shift_ = max(shift_, shiftAtDof);
- }
-
- if (comm().size() > 1)
- shift_ = comm().max(shift_);
- }
-
- /*!
- * \brief Assemble the linear system of equations \f$\mathbf{A}x - b = 0\f$.
- *
- * \param assembler The jacobian assembler
- * \param uCurrentIter The current iteration's solution vector
- */
- template<class JacobianAssembler, class SolutionVector>
- void assembleLinearSystem(JacobianAssembler& assembler,
- const SolutionVector& uCurrentIter)
- {
- assembler.assembleJacobianAndResidual(uCurrentIter);
- }
-
- /*!
- * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
- *
- * \throws Dumux::NumericalProblem if the linear solver didn't
- * converge.
- *
- * \param ls The linear solver to be used
- * \param A The matrix of the linear system of equations
- * \param x The vector which solves the linear system
- * \param b The right hand side of the linear system
- */
- template<class LinearSolver, class JacobianMatrix, class SolutionVector>
- void solveLinearSystem(LinearSolver& ls,
- JacobianMatrix& A,
- SolutionVector& x,
- SolutionVector& b)
- {
- try {
- if (numSteps_ == 0)
- {
- Scalar norm2 = b.two_norm2();
- if (comm().size() > 1)
- norm2 = comm().sum(norm2);
-
- using std::sqrt;
- initialResidual_ = sqrt(norm2);
- }
-
- //! Copy into a standard block vector.
- //! This is necessary for all model _not_ using a FieldVector<Scalar, blockSize> 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 = JacobianMatrix::block_type::rows;
- using BlockType = Dune::FieldVector<Scalar, blockSize>;
- Dune::BlockVector<BlockType> xTmp; xTmp.resize(b.size());
- Dune::BlockVector<BlockType> 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];
-
- 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];
-
- // make sure all processes converged
- int convergedRemote = converged;
- if (comm().size() > 1)
- convergedRemote = comm().min(converged);
-
- if (!converged) {
- DUNE_THROW(NumericalProblem,
- "Linear solver did not converge");
- }
- else if (!convergedRemote) {
- DUNE_THROW(NumericalProblem,
- "Linear solver did not converge on a remote process");
- }
- }
- catch (const Dune::Exception &e) {
- // make sure all processes converged
- int converged = 0;
- if (comm().size() > 1)
- converged = comm().min(converged);
-
- NumericalProblem p;
- p.message(e.what());
- throw p;
- }
- }
-
- /*!
- * \brief Update the current solution with a delta vector.
- *
- * The error estimates required for the newtonConverged() and
- * newtonProceed() methods should be updated inside this method.
- *
- * Different update strategies, such as line search and chopped
- * updates can be implemented. The default behavior is just to
- * subtract deltaU from uLastIter, i.e.
- * \f[ u^{k+1} = u^k - \Delta u^k \f]
- *
- * \param assembler The assembler (needed for global residual evaluation)
- * \param uCurrentIter The solution vector after the current iteration
- * \param uLastIter The solution vector after the last iteration
- * \param deltaU The delta as calculated from solving the linear
- * system of equations. This parameter also stores
- * the updated solution.
- */
- template<class JacobianAssembler, class SolutionVector>
- void newtonUpdate(JacobianAssembler& assembler,
- SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- if (enableShiftCriterion_)
- newtonUpdateShift(uLastIter, deltaU);
-
- if (useLineSearch_)
- {
- lineSearchUpdate_(assembler, uCurrentIter, uLastIter, deltaU);
- }
- else {
- for (unsigned int i = 0; i < uLastIter.size(); ++i) {
- uCurrentIter[i] = uLastIter[i];
- uCurrentIter[i] -= deltaU[i];
- }
-
- if (enableResidualCriterion_)
- {
- residualNorm_ = assembler.residualNorm(uCurrentIter);
- reduction_ = residualNorm_;
- reduction_ /= initialResidual_;
- }
- else
- {
- // If we get here, the convergence criterion does not require
- // additional residual evalutions. Thus, the grid variables have
- // not yet been updated to the new uCurrentIter.
- assembler.gridVariables().update(uCurrentIter);
- }
- }
- }
-
- /*!
- * \brief Indicates that one Newton iteration was finished.
- *
- * \param assembler The jacobian assembler
- * \param uCurrentIter The solution after the current Newton iteration
- * \param uLastIter The solution at the beginning of the current Newton iteration
- */
- template<class JacobianAssembler, class SolutionVector>
- void newtonEndStep(JacobianAssembler& assembler,
- SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter)
- {
- ++numSteps_;
-
- if (verbose())
- {
- std::cout << "\rNewton iteration " << numSteps_ << " done";
- if (enableShiftCriterion_)
- std::cout << ", maximum relative shift = " << shift_;
- if (enableResidualCriterion_ && enableAbsoluteResidualCriterion_)
- std::cout << ", residual = " << residualNorm_;
- else if (enableResidualCriterion_)
- std::cout << ", residual reduction = " << reduction_;
- std::cout << endIterMsg().str() << "\n";
- }
- endIterMsgStream_.str("");
-
- // When the Newton iterations are done: ask the model to check whether it makes sense
- // TODO: how do we realize this? -> do this here in the newton controller
- // model_().checkPlausibility();
- }
-
- /*!
- * \brief Called if the Newton method ended
- * (not known yet if we failed or succeeded)
- */
- void newtonEnd() {}
-
- /*!
- * \brief Called if the Newton method ended succcessfully
- * This method is called _after_ newtonEnd()
- */
- void newtonSucceed() {}
-
- /*!
- * \brief Called if the Newton method broke down.
- * This method is called _after_ newtonEnd()
- */
- template<class Assembler, class SolutionVector>
- void newtonFail(Assembler& assembler, SolutionVector& u)
- {
- if (!assembler.isStationaryProblem())
- {
- // set solution to previous solution
- u = assembler.prevSol();
-
- // reset the grid variables to the previous solution
- assembler.gridVariables().resetTimeStep(u);
-
- if (verbose())
- {
- std::cout << "Newton solver did not converge with dt = "
- << timeLoop_->timeStepSize() << " seconds. Retrying with time step of "
- << timeLoop_->timeStepSize()/2 << " seconds\n";
- }
-
- // try again with dt = dt/2
- timeLoop_->setTimeStepSize(timeLoop_->timeStepSize()/2);
- }
- else
- DUNE_THROW(Dune::MathError, "Newton solver did not converge");
- }
-
- /*!
- * \brief Suggest a new time-step size based on the old time-step
- * size.
- *
- * The default behavior 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.
- */
- Scalar suggestTimeStepSize(Scalar oldTimeStep) const
- {
- // be aggressive reducing the time-step 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);
- }
-
- Scalar percent = Scalar(targetSteps_ - numSteps_)/targetSteps_;
- return oldTimeStep*(1.0 + percent/1.2);
- }
-
- std::ostringstream &endIterMsg()
- { return endIterMsgStream_; }
-
- /*!
- * \brief Specifies if the Newton method ought to be chatty.
- */
- void setVerbose(bool val)
- { verbose_ = val; }
-
- /*!
- * \brief Returns true if the Newton method ought to be chatty.
- */
- bool verbose() const
- { return verbose_ && comm().rank() == 0; }
-
- /*!
- * \brief Returns the parameter group
- */
- const std::string& paramGroup() const
- { return paramGroup_; }
-
-protected:
-
- //! initialize the parameters by reading from the parameter tree
- void initParams_(const std::string& group = "")
- {
- useLineSearch_ = getParamFromGroup<bool>(group, "Newton.UseLineSearch");
- enableAbsoluteResidualCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableAbsoluteResidualCriterion");
- enableShiftCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableShiftCriterion");
- enableResidualCriterion_ = getParamFromGroup<bool>(group, "Newton.EnableResidualCriterion") || enableAbsoluteResidualCriterion_;
- satisfyResidualAndShiftCriterion_ = getParamFromGroup<bool>(group, "Newton.SatisfyResidualAndShiftCriterion");
- if (!enableShiftCriterion_ && !enableResidualCriterion_)
- {
- DUNE_THROW(Dune::NotImplemented,
- "at least one of NewtonEnableShiftCriterion or "
- << "NewtonEnableResidualCriterion has to be set to true");
- }
-
- setMaxRelativeShift(getParamFromGroup<Scalar>(group, "Newton.MaxRelativeShift"));
- setMaxAbsoluteResidual(getParamFromGroup<Scalar>(group, "Newton.MaxAbsoluteResidual"));
- setResidualReduction(getParamFromGroup<Scalar>(group, "Newton.ResidualReduction"));
- setTargetSteps(getParamFromGroup<int>(group, "Newton.TargetSteps"));
- setMaxSteps(getParamFromGroup<int>(group, "Newton.MaxSteps"));
-
- verbose_ = true;
- numSteps_ = 0;
- }
-
- template<class JacobianAssembler, class SolutionVector>
- void lineSearchUpdate_(JacobianAssembler& assembler,
- SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- Scalar lambda = 1.0;
- SolutionVector tmp(uLastIter);
-
- while (true)
- {
- uCurrentIter = deltaU;
- uCurrentIter *= -lambda;
- uCurrentIter += uLastIter;
-
- residualNorm_ = assembler.residualNorm(uCurrentIter);
- reduction_ = residualNorm_;
- reduction_ /= initialResidual_;
-
- if (reduction_ < lastReduction_ || lambda <= 0.125) {
- endIterMsg() << ", residual reduction " << lastReduction_ << "->" << reduction_ << "@lambda=" << lambda;
- return;
- }
-
- // try with a smaller update
- lambda /= 2.0;
- }
- }
-
- /*!
- * \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<class PrimaryVariables>
- Scalar relativeShiftAtDof_(const PrimaryVariables &priVars1,
- const PrimaryVariables &priVars2)
- {
- 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<Scalar>(1.0,abs(priVars1[j] + priVars2[j])/2);
-
- result = max(result, eqErr);
- }
- return result;
- }
-
- //! The communication object
- Communication comm_;
-
- //! The time loop for stationary simulations
- std::shared_ptr<TimeLoop<Scalar>> timeLoop_;
-
- //! message stream to be displayed at the end of iterations
- std::ostringstream endIterMsgStream_;
-
- //! switches on/off verbosity
- bool verbose_;
-
- // shift criterion variables
- Scalar shift_;
- Scalar lastShift_;
- Scalar shiftTolerance_;
-
- // residual criterion variables
- Scalar reduction_;
- Scalar residualNorm_;
- Scalar lastReduction_;
- Scalar initialResidual_;
- Scalar reductionTolerance_;
- Scalar residualTolerance_;
-
- //! optimal number of iterations we want to achieve
- int targetSteps_;
- //! maximum number of iterations we do before giving up
- int maxSteps_;
- //! actual number of steps done so far
- int numSteps_;
-
- // further parameters
- bool enablePartialReassemble_;
- bool useLineSearch_;
- bool enableAbsoluteResidualCriterion_;
- bool enableShiftCriterion_;
- bool enableResidualCriterion_;
- bool satisfyResidualAndShiftCriterion_;
-
- //! the parameter group for getting parameters from the parameter tree
- std::string paramGroup_;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh
deleted file mode 100644
index 7264b3cedf..0000000000
--- a/dumux/nonlinear/newtonmethod.hh
+++ /dev/null
@@ -1,222 +0,0 @@
-// -*- 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 2 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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Nonlinear
- * \brief The algorithmic part of the multi dimensional newton method.
- *
- * In order to use the method you need a Newtoncontroller
- */
-#ifndef DUMUX_NEWTONMETHOD_HH
-#define DUMUX_NEWTONMETHOD_HH
-
-#include <memory>
-#include <iostream>
-
-#include <dune/common/timer.hh>
-#include <dune/istl/istlexception.hh>
-
-#include <dumux/common/exceptions.hh>
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-#include <dune/common/deprecated.hh>
-
-#warning "This file is deprecated. Use NewtonSolver instead."
-
-namespace Dumux {
-
-/*!
- * \ingroup Newton
- * \ingroup Nonlinear
- * \brief The algorithmic part of the multi dimensional newton method.
- *
- * \tparam NewtonController the controller class is the driver implementation
- * \tparam JacobianAssembler an assembler for the Jacobian and the residual
- * \tparam LinearSolver the linear solver used to solve one iteration
- */
-template <class NewtonController, class JacobianAssembler, class LinearSolver>
-class DUNE_DEPRECATED_MSG("Use NewtonSolver instead.") NewtonMethod
-{
- //! provide an interface as a form of type erasure
- //! this is the minimal requirements a convergence write passed to a newton method has to fulfill
- struct ConvergenceWriterInferface
- {
- template<class SolutionVector>
- void write(const SolutionVector &uLastIter, const SolutionVector &deltaU, const SolutionVector &residual) {}
- };
-
-public:
- NewtonMethod(std::shared_ptr<NewtonController> controller,
- std::shared_ptr<JacobianAssembler> assembler,
- std::shared_ptr<LinearSolver> linearSolver,
- const std::string& modelParamGroup = "")
- : controller_(controller)
- , assembler_(assembler)
- , linearSolver_(linearSolver)
- {
- // set the linear system (matrix & residual) in the assembler
- assembler_->setLinearSystem();
-
- // set a different default for the linear solver residual reduction
- // within the Newton the linear solver doesn't need to solve too exact
- using Scalar = typename LinearSolver::Scalar;
- linearSolver_->setResidualReduction(getParamFromGroup<Scalar>(modelParamGroup, "LinearSolver.ResidualReduction", 1e-6));
- }
-
- /*!
- * \brief Run the newton method to solve a non-linear system.
- * The controller is responsible for all the strategic decisions.
- */
- template<class SolutionVector, class ConvergenceWriter = ConvergenceWriterInferface>
- bool solve(SolutionVector& uCurrentIter, const std::unique_ptr<ConvergenceWriter>& convWriter = nullptr)
- {
- try
- {
- // the given solution is the initial guess
- SolutionVector uLastIter(uCurrentIter);
- SolutionVector deltaU(uCurrentIter);
-
- Dune::Timer assembleTimer(false);
- Dune::Timer solveTimer(false);
- Dune::Timer updateTimer(false);
-
- // tell the controller that we begin solving
- controller_->newtonBegin(uCurrentIter);
-
- // execute the method as long as the controller thinks
- // that we should do another iteration
- while (controller_->newtonProceed(uCurrentIter, controller_->newtonConverged()))
- {
- // notify the controller that we're about to start
- // a new timestep
- controller_->newtonBeginStep(uCurrentIter);
-
- // make the current solution to the old one
- if (controller_->newtonNumSteps() > 0)
- uLastIter = uCurrentIter;
-
- if (controller_->verbose()) {
- std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r";
- std::cout.flush();
- }
-
- ///////////////
- // assemble
- ///////////////
-
- // linearize the problem at the current solution
- assembleTimer.start();
- controller_->assembleLinearSystem(*assembler_, uCurrentIter);
- assembleTimer.stop();
-
- ///////////////
- // linear solve
- ///////////////
-
- // Clear the current line using an ansi escape
- // sequence. for an explanation see
- // http://en.wikipedia.org/wiki/ANSI_escape_code
- const char clearRemainingLine[] = { 0x1b, '[', 'K', 0 };
-
- if (controller_->verbose()) {
- std::cout << "\rSolve: M deltax^k = r";
- std::cout << clearRemainingLine;
- std::cout.flush();
- }
-
- // solve the resulting linear equation system
- solveTimer.start();
-
- // set the delta vector to zero before solving the linear system!
- deltaU = 0;
- // ask the controller to solve the linearized system
- controller_->solveLinearSystem(*linearSolver_,
- assembler_->jacobian(),
- deltaU,
- assembler_->residual());
- solveTimer.stop();
-
- ///////////////
- // update
- ///////////////
- if (controller_->verbose()) {
- std::cout << "\rUpdate: x^(k+1) = x^k - deltax^k";
- std::cout << clearRemainingLine;
- std::cout.flush();
- }
-
- updateTimer.start();
- // update the current solution (i.e. uOld) with the delta
- // (i.e. u). The result is stored in u
- controller_->newtonUpdate(*assembler_, uCurrentIter, uLastIter, deltaU);
- updateTimer.stop();
-
- // tell the controller that we're done with this iteration
- controller_->newtonEndStep(*assembler_, uCurrentIter, uLastIter);
-
- // if a convergence writer was specified compute residual and write output
- if (convWriter)
- {
- assembler_->assembleResidual(uCurrentIter);
- convWriter->write(uLastIter, deltaU, assembler_->residual());
- }
- }
-
- // tell controller we are done
- controller_->newtonEnd();
-
- // reset state if newton failed
- if (!controller_->newtonConverged())
- {
- controller_->newtonFail(*assembler_, uCurrentIter);
- return false;
- }
-
- // tell controller we converged successfully
- controller_->newtonSucceed();
-
- if (controller_->verbose()) {
- const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
- std::cout << "Assemble/solve/update time: "
- << assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
- << solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
- << updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
- << "\n";
- }
- return true;
-
- }
- catch (const NumericalProblem &e)
- {
- if (controller_->verbose())
- std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
- controller_->newtonFail(*assembler_, uCurrentIter);
- return false;
- }
- }
-
-private:
- std::shared_ptr<NewtonController> controller_;
- std::shared_ptr<JacobianAssembler> assembler_;
- std::shared_ptr<LinearSolver> linearSolver_;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/nonlinear/staggerednewtoncontroller.hh b/dumux/nonlinear/staggerednewtoncontroller.hh
deleted file mode 100644
index bedcf6d4ad..0000000000
--- a/dumux/nonlinear/staggerednewtoncontroller.hh
+++ /dev/null
@@ -1,299 +0,0 @@
-// -*- 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 2 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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Nonlinear
- * \ingroup StaggeredDiscretization
- * \brief A newton controller for staggered finite volume schemes
- */
-#ifndef DUMUX_STAGGERED_NEWTON_CONTROLLER_HH
-#define DUMUX_STAGGERED_NEWTON_CONTROLLER_HH
-
-#include <dune/common/indices.hh>
-#include <dune/common/hybridutilities.hh>
-#include <dune/common/fvector.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dumux/common/exceptions.hh>
-#include <dumux/nonlinear/newtoncontroller.hh>
-#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
-#include <dumux/linear/matrixconverter.hh>
-#include <dune/common/deprecated.hh>
-
-#warning "This file is deprecated. Use NewtonSolver instead."
-
-namespace Dumux {
-
-/*!
- * \ingroup Nonlinear
- * \ingroup StaggeredDiscretization
- * \brief A newton controller for staggered finite volume schemes
- */
-
-template <class Scalar,
- class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
-class DUNE_DEPRECATED_MSG("Use NewtonSolver instead.")
-StaggeredNewtonController : public NewtonController<Scalar, Comm>
-{
- using ParentType = NewtonController<Scalar, Comm>;
-
-public:
- using ParentType::ParentType;
-
- /*!
- * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
- *
- * Throws Dumux::NumericalProblem if the linear solver didn't
- * converge.
- *
- * If the linear solver doesn't accept multitype matrices we copy the matrix
- * into a 1x1 block BCRS matrix for solving.
- *
- * \param ls the linear solver
- * \param A The matrix of the linear system of equations
- * \param x The vector which solves the linear system
- * \param b The right hand side of the linear system
- */
- template<class LinearSolver, class JacobianMatrix, class SolutionVector>
- void solveLinearSystem(LinearSolver& ls,
- JacobianMatrix& A,
- SolutionVector& x,
- SolutionVector& b)
- {
- // check matrix sizes
- assert(checkMatrix_(A) && "Sub blocks of MultiType matrix have wrong sizes!");
-
- try
- {
- if (this->numSteps_ == 0)
- {
- Scalar norm2 = b.two_norm2();
- if (this->comm().size() > 1)
- norm2 = this->comm().sum(norm2);
-
- using std::sqrt;
- this->initialResidual_ = sqrt(norm2);
- }
-
- // solve by calling the appropriate implementation depending on whether the linear solver
- // is capable of handling MultiType matrices or not
- bool converged = solveLinearSystem_(ls, A, x, b,
- std::integral_constant<bool, linearSolverAcceptsMultiTypeMatrix<LinearSolver>::value>());
-
- // make sure all processes converged
- int convergedRemote = converged;
- if (this->comm().size() > 1)
- convergedRemote = this->comm().min(converged);
-
- if (!converged) {
- DUNE_THROW(NumericalProblem,
- "Linear solver did not converge");
- }
- else if (!convergedRemote) {
- DUNE_THROW(NumericalProblem,
- "Linear solver did not converge on a remote process");
- }
- }
- catch (const Dune::Exception &e) {
- // make sure all processes converged
- int converged = 0;
- if (this->comm().size() > 1)
- converged = this->comm().min(converged);
-
- NumericalProblem p;
- p.message(e.what());
- throw p;
- }
- }
-
-
-
- /*!
- * \brief Update the current solution with a delta vector.
- *
- * The error estimates required for the newtonConverged() and
- * newtonProceed() methods should be updated inside this method.
- *
- * Different update strategies, such as line search and chopped
- * updates can be implemented. The default behavior is just to
- * subtract deltaU from uLastIter, i.e.
- * \f[ u^{k+1} = u^k - \Delta u^k \f]
- *
- * \param assembler The assembler for Jacobian and residual
- * \param uCurrentIter The solution vector after the current iteration
- * \param uLastIter The solution vector after the last iteration
- * \param deltaU The delta as calculated from solving the linear
- * system of equations. This parameter also stores
- * the updated solution.
- */
- template<class JacobianAssembler, class SolutionVector>
- void newtonUpdate(JacobianAssembler& assembler,
- SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- if (this->enableShiftCriterion_)
- this->newtonUpdateShift(uLastIter, deltaU);
-
- if (this->useLineSearch_)
- {
- this->lineSearchUpdate_(assembler, uCurrentIter, uLastIter, deltaU);
- }
- else {
- using namespace Dune::Hybrid;
- forEach(integralRange(Dune::Hybrid::size(uLastIter)), [&](const auto dofTypeIdx)
- {
- for (unsigned int i = 0; i < uLastIter[dofTypeIdx].size(); ++i)
- {
- uCurrentIter[dofTypeIdx][i] = uLastIter[dofTypeIdx][i];
- uCurrentIter[dofTypeIdx][i] -= deltaU[dofTypeIdx][i];
- }
- });
-
- if (this->enableResidualCriterion_)
- {
- this->residualNorm_ = assembler.residualNorm(uCurrentIter);
- this->reduction_ = this->residualNorm_;
- this->reduction_ /= this->initialResidual_;
- }
- }
-
- // update the variables class to the new solution
- assembler.gridVariables().update(uCurrentIter);
- }
-
- /*!
- * \brief Update the maximum relative shift of the solution compared to
- * the previous iteration.
- *
- * \param uLastIter The current iterative solution
- * \param deltaU The difference between the current and the next solution
- */
- template<class SolutionVector>
- void newtonUpdateShift(const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- this->shift_ = 0;
-
- using namespace Dune::Hybrid;
- forEach(integralRange(Dune::Hybrid::size(uLastIter)), [&](const auto dofTypeIdx)
- {
- for (int i = 0; i < int(uLastIter[dofTypeIdx].size()); ++i)
- {
- auto uNewI = uLastIter[dofTypeIdx][i];
- uNewI -= deltaU[dofTypeIdx][i];
-
- Scalar shiftAtDof = this->relativeShiftAtDof_(uLastIter[dofTypeIdx][i], uNewI);
- this->shift_ = std::max(this->shift_, shiftAtDof);
- }
- });
-
- if (this->comm().size() > 1)
- this->shift_ = this->comm().max(this->shift_);
- }
-
-private:
- /*!
- * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
- *
- * Throws Dumux::NumericalProblem if the linear solver didn't
- * converge.
- *
- * Specialization for linear solvers that can handle MultiType matrices.
- *
- */
- template<class LinearSolver, class JacobianMatrix, class SolutionVector>
- bool solveLinearSystem_(LinearSolver& ls,
- JacobianMatrix& A,
- SolutionVector& x,
- SolutionVector& b,
- std::true_type)
- {
- // TODO: automatically derive the precondBlockLevel
- return ls.template solve</*precondBlockLevel=*/2>(A, x, b);
- }
-
- /*!
- * \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
- *
- * Throws Dumux::NumericalProblem if the linear solver didn't
- * converge.
- *
- * Specialization for linear solvers that cannot handle MultiType matrices.
- * We copy the matrix into a 1x1 block BCRS matrix before solving.
- *
- */
- template<class LinearSolver, class JacobianMatrix, class SolutionVector>
- bool solveLinearSystem_(LinearSolver& ls,
- JacobianMatrix& A,
- SolutionVector& x,
- SolutionVector& b,
- std::false_type)
- {
- // create the bcrs matrix the IterativeSolver backend can handle
- const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
-
- // get the new matrix sizes
- const std::size_t numRows = M.N();
- assert(numRows == M.M());
-
- // create the vector the IterativeSolver backend can handle
- const auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
- assert(bTmp.size() == numRows);
-
- // create a blockvector to which the linear solver writes the solution
- using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
- using BlockVector = typename Dune::BlockVector<VectorBlock>;
- BlockVector y(numRows);
-
- // solve
- const bool converged = ls.solve(M, y, bTmp);
-
- // copy back the result y into x
- if(converged)
- VectorConverter<SolutionVector>::retrieveValues(x, y);
-
- return converged;
- }
-
- //! helper method to assure the MultiType matrix's sub blocks have the correct sizes
- template<class JacobianMatrix>
- bool checkMatrix_(const JacobianMatrix& A)
- {
- bool matrixHasCorrectSize = true;
- using namespace Dune::Hybrid;
- using namespace Dune::Indices;
- forEach(A, [&matrixHasCorrectSize](const auto& rowOfMultiTypeMatrix)
- {
- const auto numRowsLeftMostBlock = rowOfMultiTypeMatrix[_0].N();
-
- forEach(rowOfMultiTypeMatrix, [&matrixHasCorrectSize, &numRowsLeftMostBlock](const auto& subBlock)
- {
- if (subBlock.N() != numRowsLeftMostBlock)
- matrixHasCorrectSize = false;
- });
- });
- return matrixHasCorrectSize;
- }
-
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/2pnc/CMakeLists.txt b/dumux/porousmediumflow/2pnc/CMakeLists.txt
index 612765aa72..5e3c930c89 100644
--- a/dumux/porousmediumflow/2pnc/CMakeLists.txt
+++ b/dumux/porousmediumflow/2pnc/CMakeLists.txt
@@ -3,7 +3,6 @@
install(FILES
indices.hh
model.hh
-newtoncontroller.hh
primaryvariableswitch.hh
volumevariables.hh
vtkoutputfields.hh
diff --git a/dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh b/dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh
deleted file mode 100644
index 553414dbbc..0000000000
--- a/dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh
+++ /dev/null
@@ -1,220 +0,0 @@
-// -*- 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 2 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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup PorousmediumCompositional
- * \brief Reference implementation of a controller class for the Newton solver.
- *
- * Usually this controller should be sufficient.
- */
-#ifndef DUMUX_PRIVARSWITCH_NEWTON_CONTROLLER_HH
-#define DUMUX_PRIVARSWITCH_NEWTON_CONTROLLER_HH
-
-#include <memory>
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-#include <dumux/discretization/methods.hh>
-#include <dumux/discretization/elementsolution.hh>
-#include <dumux/nonlinear/newtoncontroller.hh>
-#include <dune/common/deprecated.hh>
-
-#warning "This file is deprecated. Use PriVarSwitchNewtonSolver instead."
-
-namespace Dumux {
-
-/*!
- * \ingroup PorousmediumCompositional
- * \brief A newton controller that handles primary variable switches
- * \todo make this independent of TypeTag by making PrimaryVariableSwitch a template argument
- * and extracting everything model specific from there
- * \todo Implement for volume variable caching enabled
- */
-template <class TypeTag>
-class DUNE_DEPRECATED_MSG("Use PriVarSwitchNewtonSolver instead.")
-PriVarSwitchNewtonController : public NewtonController<typename GET_PROP_TYPE(TypeTag, Scalar)>
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using ParentType = NewtonController<Scalar>;
- using PrimaryVariableSwitch = typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
-
- static constexpr bool isBox = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::box;
-
-public:
- using ParentType::ParentType;
-
- /*!
- * \brief Returns true if the error of the solution is below the
- * tolerance.
- */
- bool newtonConverged() const
- {
- if (switchedInLastIteration_)
- return false;
-
- return ParentType::newtonConverged();
- }
-
- /*!
- *
- * \brief Called before the Newton method is applied to an
- * non-linear system of equations.
- *
- * \param u The initial solution
- */
- template<class SolutionVector>
- void newtonBegin(const SolutionVector &u)
- {
- ParentType::newtonBegin(u);
- priVarSwitch_ = std::make_unique<PrimaryVariableSwitch>(u.size());
- }
-
- /*!
- * \brief Indicates that one Newton iteration was finished.
- *
- * \param assembler The jacobian assembler
- * \param uCurrentIter The solution after the current Newton iteration
- * \param uLastIter The solution at the beginning of the current Newton iteration
- */
- template<class JacobianAssembler, class SolutionVector>
- void newtonEndStep(JacobianAssembler& assembler,
- SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter)
- {
- ParentType::newtonEndStep(assembler, uCurrentIter, uLastIter);
-
- // update the variable switch (returns true if the pri vars at at least one dof were switched)
- // for disabled grid variable caching
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- const auto& problem = assembler.problem();
- auto& gridVariables = assembler.gridVariables();
-
- // invoke the primary variable switch
- switchedInLastIteration_ = priVarSwitch_->update(uCurrentIter, gridVariables,
- problem, fvGridGeometry);
-
- if(switchedInLastIteration_)
- {
- for (const auto& element : elements(fvGridGeometry.gridView()))
- {
- // if the volume variables are cached globally, we need to update those where the primary variables have been switched
- updateSwitchedVolVars_(std::integral_constant<bool, GET_PROP_VALUE(TypeTag, EnableGridVolumeVariablesCache)>(),
- element, assembler, uCurrentIter, uLastIter);
-
- // if the flux variables are cached globally, we need to update those where the primary variables have been switched
- // (not needed for box discretization)
- updateSwitchedFluxVarsCache_(std::integral_constant<bool, (GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache) && !isBox)>(),
- element, assembler, uCurrentIter, uLastIter);
- }
- }
- }
-
- /*!
- * \brief Called if the Newton method ended
- * (not known yet if we failed or succeeded)
- */
- void newtonEnd()
- {
- ParentType::newtonEnd();
-
- // in any way, we have to reset the switch flag
- switchedInLastIteration_ = false;
- // free some memory
- priVarSwitch_.release();
- }
-
-private:
-
- /*!
- * \brief Update the volume variables whose primary variables were
- switched. Required when volume variables are cached globally.
- */
- template<class Element, class JacobianAssembler, class SolutionVector>
- void updateSwitchedVolVars_(std::true_type,
- const Element& element,
- JacobianAssembler& assembler,
- const SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter)
- {
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- const auto& problem = assembler.problem();
- auto& gridVariables = assembler.gridVariables();
-
- // make sure FVElementGeometry is bound to the element
- auto fvGeometry = localView(fvGridGeometry);
- fvGeometry.bindElement(element);
-
- // update the secondary variables if global caching is enabled
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto dofIdxGlobal = scv.dofIndex();
- if (priVarSwitch_->wasSwitched(dofIdxGlobal))
- {
- const auto elemSol = elementSolution(element, uCurrentIter, fvGridGeometry);
- auto& volVars = gridVariables.curGridVolVars().volVars(scv);
- volVars.update(elemSol, problem, element, scv);
- }
- }
- }
-
- /*!
- * \brief Update the fluxVars cache for dof whose primary variables were
- switched. Required when flux variables are cached globally.
- */
- template<class Element, class JacobianAssembler, class SolutionVector>
- void updateSwitchedFluxVarsCache_(std::true_type,
- const Element& element,
- JacobianAssembler& assembler,
- const SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter)
- {
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& gridVariables = assembler.gridVariables();
-
- // update the flux variables if global caching is enabled
- const auto dofIdxGlobal = fvGridGeometry.dofMapper().index(element);
-
- if (priVarSwitch_->wasSwitched(dofIdxGlobal))
- {
- // make sure FVElementGeometry and the volume variables are bound
- auto fvGeometry = localView(fvGridGeometry);
- fvGeometry.bind(element);
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, uCurrentIter);
- gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
- }
- }
-
- //! brief Do nothing when volume variables are not cached globally.
- template <typename... Args>
- void updateSwitchedVolVars_(std::false_type, Args&&... args) const {}
-
- //! brief Do nothing when flux variables are not cached globally.
- template <typename... Args>
- void updateSwitchedFluxVarsCache_(std::false_type, Args&&... args) const {}
-
- //! the class handling the primary variable switch
- std::unique_ptr<PrimaryVariableSwitch> priVarSwitch_;
- //! if we switched primary variables in the last iteration
- bool switchedInLastIteration_ = false;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/nonequilibrium/CMakeLists.txt b/dumux/porousmediumflow/nonequilibrium/CMakeLists.txt
index 31fa8eb04f..09252e0484 100644
--- a/dumux/porousmediumflow/nonequilibrium/CMakeLists.txt
+++ b/dumux/porousmediumflow/nonequilibrium/CMakeLists.txt
@@ -4,8 +4,8 @@ install(FILES
localresidual.hh
gridvariables.hh
indices.hh
-newtoncontroller.hh
volumevariables.hh
vtkoutputfields.hh
model.hh
+newtonsolver.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/nonequilibrium)
diff --git a/dumux/porousmediumflow/nonequilibrium/newtoncontroller.hh b/dumux/porousmediumflow/nonequilibrium/newtoncontroller.hh
deleted file mode 100644
index b1bfec283b..0000000000
--- a/dumux/porousmediumflow/nonequilibrium/newtoncontroller.hh
+++ /dev/null
@@ -1,66 +0,0 @@
-// -*- 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 2 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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup PorousmediumNonEquilibriumModel
- * \brief A MpNc specific controller for the newton solver.
- * This controller calls the velocity averaging in the model after each iteration.
- */
-#ifndef DUMUX_NONEQUILIBRIUM_NEWTON_CONTROLLER_HH
-#define DUMUX_NONEQUILIBRIUM_NEWTON_CONTROLLER_HH
-
-#include <algorithm>
-
-#include <dumux/nonlinear/newtoncontroller.hh>
-
-namespace Dumux {
-/*!
- * \ingroup PorousmediumNonEquilibriumModel
- * \brief A nonequilibrium specific controller for the newton solver.
- * This controller calls the velocity averaging in the problem after each iteration.
- */
-template <class Scalar,
- class Comm = Dune::CollectiveCommunication<Dune::MPIHelper::MPICommunicator> >
-class NonEquilibriumNewtonController : public NewtonController<Scalar, Comm>
-{
- using ParentType = NewtonController<Scalar, Comm>;
-
-public:
- using ParentType::ParentType;
-
- template<class JacobianAssembler, class SolutionVector>
- void newtonUpdate(JacobianAssembler& assembler,
- SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- ParentType::newtonUpdate(assembler,
- uCurrentIter,
- uLastIter,
- deltaU);
-
- auto& gridVariables = assembler.gridVariables();
- // Averages the face velocities of a vertex. Implemented in the model.
- // The velocities are stored in the model.
- gridVariables.calcVelocityAverage(uCurrentIter);
- }
-};
-
-} // end namespace Dumux
-#endif // DUMUX_VELO_PROB_AVERAGE_NEWTON_CONTROLLER_HH
diff --git a/dumux/porousmediumflow/richards/CMakeLists.txt b/dumux/porousmediumflow/richards/CMakeLists.txt
index 2929d226a4..199fc1455f 100644
--- a/dumux/porousmediumflow/richards/CMakeLists.txt
+++ b/dumux/porousmediumflow/richards/CMakeLists.txt
@@ -4,7 +4,8 @@ install(FILES
indices.hh
localresidual.hh
model.hh
-newtoncontroller.hh
+newtonsolver.hh
+privarswitchnewtonsolver.hh
primaryvariableswitch.hh
volumevariables.hh
vtkoutputfields.hh
diff --git a/dumux/porousmediumflow/richards/newtoncontroller.hh b/dumux/porousmediumflow/richards/newtoncontroller.hh
deleted file mode 100644
index e6a7cccd48..0000000000
--- a/dumux/porousmediumflow/richards/newtoncontroller.hh
+++ /dev/null
@@ -1,122 +0,0 @@
-// -*- 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 2 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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup RichardsModel
- * \brief A Richards model specific controller for the newton solver.
- */
-#ifndef DUMUX_RICHARDS_NEWTON_CONTROLLER_HH
-#define DUMUX_RICHARDS_NEWTON_CONTROLLER_HH
-
-#include <dumux/common/properties.hh>
-#include <dumux/nonlinear/newtoncontroller.hh>
-#include <dumux/discretization/elementsolution.hh>
-#include <dune/common/deprecated.hh>
-
-#warning "This file is deprecated. Use RichardsNewtonSolver instead."
-
-namespace Dumux {
-/*!
- * \ingroup RichardsModel
- * \brief A Richards model specific controller for the newton solver.
- *
- * This controller 'knows' what a 'physically meaningful' solution is
- * and can thus do update smarter than the plain Newton controller.
- *
- * \todo make this typetag independent by extracting anything model specific from assembler
- * or from possible ModelTraits.
- */
-template <class TypeTag>
-class DUNE_DEPRECATED_MSG("Use RichardsNewtonSolver instead.")
-RichardsNewtonController : public NewtonController<typename GET_PROP_TYPE(TypeTag, Scalar)>
-{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using ParentType = NewtonController<Scalar>;
-
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- enum { pressureIdx = Indices::pressureIdx };
-
-public:
- using ParentType::ParentType;
-
- /*!
- * \brief Update the current solution of the newton method
- *
- * This is basically the step
- * \f[ u^{k+1} = u^k - \Delta u^k \f]
- *
- * \param assembler The Jacobian assembler
- * \param uCurrentIter The solution after the current Newton iteration \f$ u^{k+1} \f$
- * \param uLastIter The solution after the last Newton iteration \f$ u^k \f$
- * \param deltaU The vector of differences between the last
- * iterative solution and the next one \f$ \Delta u^k \f$
- */
- template<class JacobianAssembler, class SolutionVector>
- void newtonUpdate(JacobianAssembler& assembler,
- SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter,
- const SolutionVector &deltaU)
- {
- ParentType::newtonUpdate(assembler, uCurrentIter, uLastIter, deltaU);
- if (!this->useLineSearch_ && getParamFromGroup<bool>(this->paramGroup(), "Newton.EnableChop"))
- {
- // do not clamp anything after 5 iterations
- if (this->numSteps_ > 4)
- return;
-
- // clamp saturation change to at most 20% per iteration
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- for (const auto& element : elements(fvGridGeometry.gridView()))
- {
- auto fvGeometry = localView(fvGridGeometry);
- fvGeometry.bindElement(element);
-
- for (auto&& scv : scvs(fvGeometry))
- {
- auto dofIdxGlobal = scv.dofIndex();
-
- // calculate the old wetting phase saturation
- const auto& spatialParams = assembler.problem().spatialParams();
- const auto elemSol = elementSolution(element, uCurrentIter, fvGridGeometry);
- const auto& materialLawParams = spatialParams.materialLawParams(element, scv, elemSol);
- const Scalar pcMin = MaterialLaw::pc(materialLawParams, 1.0);
- const Scalar pw = uLastIter[dofIdxGlobal][pressureIdx];
- using std::max;
- const Scalar pn = max(assembler.problem().nonWettingReferencePressure(), pw + pcMin);
- const Scalar pcOld = pn - pw;
- const Scalar SwOld = max(0.0, MaterialLaw::sw(materialLawParams, pcOld));
-
- // convert into minimum and maximum wetting phase
- // pressures
- const Scalar pwMin = pn - MaterialLaw::pc(materialLawParams, SwOld - 0.2);
- const Scalar pwMax = pn - MaterialLaw::pc(materialLawParams, SwOld + 0.2);
-
- // clamp the result
- using std::min; using std::max;
- uCurrentIter[dofIdxGlobal][pressureIdx] = max(pwMin, min(uCurrentIter[dofIdxGlobal][pressureIdx], pwMax));
- }
- }
- }
- }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/richardsnc/model.hh b/dumux/porousmediumflow/richardsnc/model.hh
index 659a42feb0..91720a3e1c 100644
--- a/dumux/porousmediumflow/richardsnc/model.hh
+++ b/dumux/porousmediumflow/richardsnc/model.hh
@@ -70,7 +70,6 @@
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/compositional/localresidual.hh>
-#include <dumux/porousmediumflow/richards/newtoncontroller.hh>
#include <dumux/material/spatialparams/fv1p.hh>
#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
--
GitLab