diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index aa7dfb5eb60fcd2142f6f2717565cc4fc2017bc4..ee493a2b21049f24e0759e716d30e9f6d388afdb 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -39,7 +39,6 @@
#include <dumux/common/parameters.hh>
#include <dumux/common/exceptions.hh>
-#include <dumux/common/timeloop.hh>
#include <dumux/common/typetraits/vector.hh>
#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
@@ -110,38 +109,6 @@ public:
}
}
- /*!
- * \brief Constructor for instationary problems
- */
- NewtonSolver(std::shared_ptr<Assembler> assembler,
- std::shared_ptr<LinearSolver> linearSolver,
- std::shared_ptr<TimeLoop<Scalar>> timeLoop,
- const Communication& comm = Dune::MPIHelper::getCollectiveCommunication(),
- const std::string& paramGroup = "")
- : endIterMsgStream_(std::ostringstream::out)
- , assembler_(assembler)
- , linearSolver_(linearSolver)
- , comm_(comm)
- , timeLoop_(timeLoop)
- , paramGroup_(paramGroup)
- {
- initParams_(paramGroup);
-
- // 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
- linearSolver_->setResidualReduction(getParamFromGroup<Scalar>(paramGroup, "LinearSolver.ResidualReduction", 1e-6));
-
- // initialize the partial reassembler
- if (enablePartialReassembly_)
- {
- partialReassembler_ = std::make_unique<Reassembler>(assembler_->fvGridGeometry());
- distanceFromLastLinearization_.resize(assembler_->fvGridGeometry().numDofs(), 0);
- }
- }
-
//! the communicator for parallel runs
const Communication& comm() const
{ return comm_; }
@@ -198,137 +165,62 @@ public:
/*!
* \brief Run the newton method to solve a non-linear system.
- * The controller is responsible for all the strategic decisions.
+ * Does time step control when the newton fails to converge
*/
- bool solve(SolutionVector& uCurrentIter, const std::unique_ptr<ConvergenceWriter>& convWriter = nullptr)
+ template<class TimeLoop>
+ void solve(SolutionVector& uCurrentIter, TimeLoop& timeLoop,
+ const std::unique_ptr<ConvergenceWriter>& convWriter = nullptr)
{
- try
- {
- // the given solution is the initial guess
- SolutionVector uLastIter(uCurrentIter);
- SolutionVector deltaU(uCurrentIter);
+ if (assembler_->isStationaryProblem())
+ DUNE_THROW(Dune::InvalidStateException, "Using time step control with stationary problem makes no sense!");
- Dune::Timer assembleTimer(false);
- Dune::Timer solveTimer(false);
- Dune::Timer updateTimer(false);
+ // try solving the non-linear system
+ for (std::size_t i = 0; i <= maxTimeStepDivisions_; ++i)
+ {
+ // linearize & solve
+ const bool converged = solve_(uCurrentIter, convWriter);
- if (enablePartialReassembly_)
- {
- partialReassembler_->resetColors();
- std::fill(distanceFromLastLinearization_.begin(),
- distanceFromLastLinearization_.end(), 0.0);
- }
- newtonBegin(uCurrentIter);
+ if (converged)
+ return;
- // execute the method as long as the controller thinks
- // that we should do another iteration
- while (newtonProceed(uCurrentIter, newtonConverged()))
+ else if (!converged && i < maxTimeStepDivisions_)
{
- // notify the controller that we're about to start
- // a new timestep
- newtonBeginStep(uCurrentIter);
-
- // make the current solution to the old one
- if (numSteps_ > 0)
- uLastIter = uCurrentIter;
+ // set solution to previous solution
+ uCurrentIter = assembler_->prevSol();
- if (verbose_) {
- std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r"
- << std::flush;
- }
+ // reset the grid variables to the previous solution
+ assembler_->gridVariables().resetTimeStep(uCurrentIter);
- ///////////////
- // assemble
- ///////////////
+ if (verbose_)
+ std::cout << "Newton solver did not converge with dt = "
+ << timeLoop.timeStepSize() << " seconds. Retrying with time step of "
+ << timeLoop.timeStepSize()/2 << " seconds\n";
- // linearize the problem at the current solution
- assembleTimer.start();
- assembleLinearSystem(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 (verbose_) {
- std::cout << "\rSolve: M deltax^k = r";
- std::cout << clearRemainingLine
- << std::flush;
- }
-
- // solve the resulting linear equation system
- solveTimer.start();
-
- // set the delta vector to zero before solving the linear system!
- deltaU = 0;
-
- solveLinearSystem(deltaU);
- solveTimer.stop();
-
- ///////////////
- // update
- ///////////////
- if (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
- newtonUpdate(uCurrentIter, uLastIter, deltaU);
- updateTimer.stop();
-
- // tell the controller that we're done with this iteration
- newtonEndStep(uCurrentIter, uLastIter);
-
- // if a convergence writer was specified compute residual and write output
- if (convWriter)
- {
- assembler_->assembleResidual(uCurrentIter);
- convWriter->write(uLastIter, deltaU, assembler_->residual());
- }
+ // try again with dt = dt/2
+ timeLoop.setTimeStepSize(timeLoop.timeStepSize()/2);
}
- // tell controller we are done
- newtonEnd();
-
- // reset state if newton failed
- if (!newtonConverged())
+ else
{
- newtonFail(uCurrentIter);
- return false;
+ DUNE_THROW(NumericalProblem, "Newton solver didn't converge after "
+ << maxTimeStepDivisions_ << " time-step divisions. dt="
+ << timeLoop.timeStepSize() << '\n');
}
-
- // tell controller we converged successfully
- newtonSucceed();
-
- if (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 (verbose_)
- std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
- newtonFail(uCurrentIter);
- return false;
}
}
+ /*!
+ * \brief Run the newton method to solve a non-linear system.
+ * The controller is responsible for all the strategic decisions.
+ */
+ void solve(SolutionVector& uCurrentIter, const std::unique_ptr<ConvergenceWriter>& convWriter = nullptr)
+ {
+ const bool converged = solve_(uCurrentIter, convWriter);
+ if (!converged)
+ DUNE_THROW(NumericalProblem, "Newton solver didn't converge after "
+ << numSteps_ << " iterations.\n");
+ }
+
/*!
* \brief Called before the Newton method is applied to an
* non-linear system of equations.
@@ -626,29 +518,7 @@ public:
* \brief Called if the Newton method broke down.
* This method is called _after_ newtonEnd()
*/
- virtual void newtonFail(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");
- }
+ virtual void newtonFail(SolutionVector& u) {}
/*!
* \brief Called if the Newton method ended succcessfully
@@ -722,9 +592,6 @@ protected:
Assembler& assembler()
{ return *assembler_; }
- const TimeLoop<Scalar>& timeLoop() const
- { return *timeLoop_; }
-
//! optimal number of iterations we want to achieve
int targetSteps_;
//! maximum number of iterations we do before giving up
@@ -748,6 +615,140 @@ protected:
private:
+ /*!
+ * \brief Run the newton method to solve a non-linear system.
+ * The controller is responsible for all the strategic decisions.
+ */
+ bool solve_(SolutionVector& uCurrentIter, const std::unique_ptr<ConvergenceWriter>& convWriter = nullptr)
+ {
+ // 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);
+
+ if (enablePartialReassembly_)
+ {
+ partialReassembler_->resetColors();
+ std::fill(distanceFromLastLinearization_.begin(),
+ distanceFromLastLinearization_.end(), 0.0);
+ }
+
+ try
+ {
+ newtonBegin(uCurrentIter);
+
+ // execute the method as long as the controller thinks
+ // that we should do another iteration
+ while (newtonProceed(uCurrentIter, newtonConverged()))
+ {
+ // notify the controller that we're about to start
+ // a new timestep
+ newtonBeginStep(uCurrentIter);
+
+ // make the current solution to the old one
+ if (numSteps_ > 0)
+ uLastIter = uCurrentIter;
+
+ if (verbose_) {
+ std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r"
+ << std::flush;
+ }
+
+ ///////////////
+ // assemble
+ ///////////////
+
+ // linearize the problem at the current solution
+ assembleTimer.start();
+ assembleLinearSystem(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 (verbose_) {
+ std::cout << "\rSolve: M deltax^k = r";
+ std::cout << clearRemainingLine
+ << std::flush;
+ }
+
+ // solve the resulting linear equation system
+ solveTimer.start();
+
+ // set the delta vector to zero before solving the linear system!
+ deltaU = 0;
+
+ solveLinearSystem(deltaU);
+ solveTimer.stop();
+
+ ///////////////
+ // update
+ ///////////////
+ if (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
+ newtonUpdate(uCurrentIter, uLastIter, deltaU);
+ updateTimer.stop();
+
+ // tell the controller that we're done with this iteration
+ newtonEndStep(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
+ newtonEnd();
+
+ // reset state if newton failed
+ if (!newtonConverged())
+ {
+ newtonFail(uCurrentIter);
+ return false;
+ }
+
+ // tell controller we converged successfully
+ newtonSucceed();
+
+ if (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 (verbose_)
+ std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
+ newtonFail(uCurrentIter);
+ return false;
+ }
+ }
+
//! assembleLinearSystem_ for assemblers that support partial reassembly
template<class A = Assembler>
auto assembleLinearSystem_(const SolutionVector& uCurrentIter)
@@ -1017,6 +1018,8 @@ private:
reassemblyMaxThreshold_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyMaxThreshold", 1e2*shiftTolerance_);
reassemblyShiftWeight_ = getParamFromGroup<Scalar>(group, "Newton.ReassemblyShiftWeight", 1e-3);
+ maxTimeStepDivisions_ = getParamFromGroup<std::size_t>(group, "Newton.MaxTimeStepDivisions", 10);
+
verbose_ = comm_.rank() == 0;
numSteps_ = 0;
}
@@ -1073,9 +1076,6 @@ private:
//! The communication object
Communication comm_;
- //! The time loop for stationary simulations
- std::shared_ptr<TimeLoop<Scalar>> timeLoop_;
-
//! switches on/off verbosity
bool verbose_;
@@ -1083,6 +1083,9 @@ private:
Scalar reductionTolerance_;
Scalar residualTolerance_;
+ // time step control
+ std::size_t maxTimeStepDivisions_;
+
// further parameters
bool useLineSearch_;
bool useChop_;
diff --git a/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc b/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc
index 0109b9ccab82e954b06e47a3ede3965bd62bf6e8..9712deaf8ecf1a6f19ca420258dc19e153e48bed 100644
--- a/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc
+++ b/test/porousmediumflow/1pnc/implicit/test_1p2c_fv.cc
@@ -107,7 +107,6 @@
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
- auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
// intialize the vtk output module
@@ -129,7 +128,7 @@
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
- NewtonSolver<Assembler, LinearSolver> nonLinearSolver(assembler, linearSolver, timeLoop);
+ NewtonSolver<Assembler, LinearSolver> nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
@@ -137,24 +136,8 @@
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
- // try solving the non-linear system
- for (int i = 0; i < maxDivisions; ++i)
- {
- // linearize & solve
- auto converged = nonLinearSolver.solve(x);
-
- if (converged)
- break;
-
- if (!converged && i == maxDivisions-1)
- DUNE_THROW(Dune::MathError,
- "Newton solver didn't converge after "
- << maxDivisions
- << " time-step divisions. dt="
- << timeLoop->timeStepSize()
- << ".\nThe solutions of the current and the previous time steps "
- << "have been saved to restart files.");
- }
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
diff --git a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
index 36fdc3307aa0042ce4c3f7b7dab36d22594bf99a..ad74e0d38dcb404a03cd090a422909e235c94c8d 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
+++ b/test/porousmediumflow/2p/implicit/incompressible/test_2p_fv.cc
@@ -134,7 +134,6 @@ int main(int argc, char** argv) try
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -163,7 +162,7 @@ int main(int argc, char** argv) try
// the non-linear solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
- NewtonSolver nonLinearSolver(assembler, linearSolver, timeLoop);
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
@@ -171,24 +170,8 @@ int main(int argc, char** argv) try
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
- // try solving the non-linear system
- for (int i = 0; i < maxDivisions; ++i)
- {
- // linearize & solve
- auto converged = nonLinearSolver.solve(x);
-
- if (converged)
- break;
-
- if (!converged && i == maxDivisions-1)
- DUNE_THROW(Dune::MathError,
- "Newton solver didn't converge after "
- << maxDivisions
- << " time-step divisions. dt="
- << timeLoop->timeStepSize()
- << ".\nThe solutions of the current and the previous time steps "
- << "have been saved to restart files.");
- }
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
diff --git a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
index 48414c9f14d4142bacce3b9909f9894ff9bae0b1..efac6ec2c78eb6e0fcbd2a5e7eda5fae2950e080 100644
--- a/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
+++ b/test/porousmediumflow/2p2c/implicit/test_2p2c_fv.cc
@@ -103,7 +103,6 @@ int main(int argc, char** argv) try
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -140,24 +139,8 @@ int main(int argc, char** argv) try
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
- // try solving the non-linear system
- for (int i = 0; i < maxDivisions; ++i)
- {
- // linearize & solve
- auto converged = nonLinearSolver.solve(x);
-
- if (converged)
- break;
-
- if (!converged && i == maxDivisions-1)
- DUNE_THROW(Dune::MathError,
- "Newton solver didn't converge after "
- << maxDivisions
- << " time-step divisions. dt="
- << timeLoop->timeStepSize()
- << ".\nThe solutions of the current and the previous time steps "
- << "have been saved to restart files.");
- }
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;
diff --git a/test/porousmediumflow/richards/implicit/test_richardslens_fv.cc b/test/porousmediumflow/richards/implicit/test_richardslens_fv.cc
index c108a3c94e90da4bcdae6cdfbb6b736fa34e0519..9d260e610e4532a6cdc41a9902cce78c46d5f890 100644
--- a/test/porousmediumflow/richards/implicit/test_richardslens_fv.cc
+++ b/test/porousmediumflow/richards/implicit/test_richardslens_fv.cc
@@ -127,7 +127,6 @@ int main(int argc, char** argv) try
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
- const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
@@ -156,7 +155,7 @@ int main(int argc, char** argv) try
// the non-linear solver
using NewtonSolver = Dumux::RichardsNewtonSolver<TypeTag, Assembler, LinearSolver>;
- NewtonSolver nonLinearSolver(assembler, linearSolver, timeLoop);
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
@@ -164,24 +163,8 @@ int main(int argc, char** argv) try
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
- // try solving the non-linear system
- for (int i = 0; i < maxDivisions; ++i)
- {
- // linearize & solve
- auto converged = nonLinearSolver.solve(x);
-
- if (converged)
- break;
-
- if (!converged && i == maxDivisions-1)
- DUNE_THROW(Dune::MathError,
- "Newton solver didn't converge after "
- << maxDivisions
- << " time-step divisions. dt="
- << timeLoop->timeStepSize()
- << ".\nThe solutions of the current and the previous time steps "
- << "have been saved to restart files.");
- }
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(x, *timeLoop);
// make the new solution the old solution
xOld = x;