27 #ifndef EWOMS_NEWTON_METHOD_HH
28 #define EWOMS_NEWTON_METHOD_HH
37 #include <opm/material/densead/Math.hpp>
39 #include <opm/common/Unused.hpp>
40 #include <opm/common/Exceptions.hpp>
41 #include <opm/common/ErrorMacros.hpp>
43 #include <dune/istl/istlexception.hh>
44 #include <dune/common/classname.hh>
45 #include <dune/common/version.hh>
46 #include <dune/common/parallel/mpihelper.hh>
55 template <
class TypeTag>
61 namespace Properties {
169 template <
class TypeTag>
172 typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) Implementation;
173 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
174 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
175 typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
176 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
178 typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
179 typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
180 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
181 typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
182 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
183 typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
184 typedef typename GET_PROP_TYPE(TypeTag, JacobianMatrix) JacobianMatrix;
185 typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) LinearSolverBackend;
186 typedef typename GET_PROP_TYPE(TypeTag, NewtonConvergenceWriter) ConvergenceWriter;
188 typedef typename Dune::MPIHelper::MPICommunicator Communicator;
189 typedef Dune::CollectiveCommunication<Communicator> CollectiveCommunication;
192 NewtonMethod(Simulator& simulator)
193 : simulator_(simulator)
194 , endIterMsgStream_(std::ostringstream::out)
195 , linearSolver_(simulator)
196 , comm_(Dune::MPIHelper::getCommunicator())
197 , convergenceWriter_(asImp_())
211 LinearSolverBackend::registerParameters();
214 "Specify whether the Newton method should inform "
215 "the user about its progress or not");
217 "Write the convergence behaviour of the Newton "
218 "method to a VTK file");
220 "The 'optimum' number of Newton iterations per "
223 "The maximum number of Newton iterations per time "
226 "The maximum raw error tolerated by the Newton"
227 "method for considering a solution to be "
230 "The maximum error tolerated by the Newton "
231 "method to which does not cause an abort");
245 {
return simulator_.
problem(); }
251 {
return simulator_.
problem(); }
257 {
return simulator_.
model(); }
263 {
return simulator_.
model(); }
270 {
return numIterations_; }
280 { numIterations_ = value; }
287 {
return tolerance_; }
294 { tolerance_ = value; }
307 const char *clearRemainingLine =
"\n";
308 if (isatty(fileno(stdout))) {
309 static const char blubb[] = { 0x1b,
'[',
'K',
'\r', 0 };
310 clearRemainingLine = blubb;
314 prePostProcessTimer_.
halt();
315 linearizeTimer_.
halt();
319 SolutionVector& nextSolution =
model().solution(0);
320 SolutionVector currentSolution(nextSolution);
321 GlobalEqVector solutionUpdate(nextSolution.size());
323 Linearizer& linearizer =
model().linearizer();
328 prePostProcessTimer_.
start();
329 asImp_().begin_(nextSolution);
330 prePostProcessTimer_.
stop();
345 prePostProcessTimer_.
start();
346 asImp_().beginIteration_();
347 prePostProcessTimer_.
stop();
350 currentSolution = nextSolution;
353 std::cout <<
"Linearize: r(x^k) = dS/dt + div F - q; M = grad r"
354 << clearRemainingLine
359 linearizeTimer_.
start();
360 asImp_().linearize_();
361 linearizeTimer_.
stop();
366 updateTimer_.
start();
367 auto& M = linearizer.matrix();
368 auto& b = linearizer.residual();
369 linearSolver_.prepareRhs(M, b);
370 asImp_().preSolve_(currentSolution, b);
374 if (asImp_().
verbose_() && isatty(fileno(stdout)))
375 std::cout << clearRemainingLine
379 prePostProcessTimer_.
start();
380 asImp_().endIteration_(nextSolution, currentSolution);
381 prePostProcessTimer_.
stop();
388 std::cout <<
"Solve: M deltax^k = r"
389 << clearRemainingLine
395 linearSolver_.prepareMatrix(M);
396 bool converged = linearSolver_.solve(solutionUpdate);
402 std::cout <<
"Newton: Linear solver did not converge\n" << std::flush;
404 prePostProcessTimer_.
start();
406 prePostProcessTimer_.
stop();
413 std::cout <<
"Update: x^(k+1) = x^k - deltax^k"
414 << clearRemainingLine
420 updateTimer_.
start();
421 asImp_().postSolve_(nextSolution,
425 asImp_().update_(nextSolution, currentSolution, solutionUpdate, b);
428 if (asImp_().
verbose_() && isatty(fileno(stdout)))
430 std::cout << clearRemainingLine
434 prePostProcessTimer_.
start();
435 asImp_().endIteration_(nextSolution, currentSolution);
436 prePostProcessTimer_.
stop();
439 catch (
const Dune::Exception& e)
442 std::cout <<
"Newton method caught exception: \""
443 << e.what() <<
"\"\n" << std::flush;
445 prePostProcessTimer_.
start();
447 prePostProcessTimer_.
stop();
451 catch (
const Opm::NumericalProblem& e)
454 std::cout <<
"Newton method caught exception: \""
455 << e.what() <<
"\"\n" << std::flush;
457 prePostProcessTimer_.
start();
459 prePostProcessTimer_.
stop();
465 if (asImp_().
verbose_() && isatty(fileno(stdout)))
466 std::cout << clearRemainingLine
470 prePostProcessTimer_.
start();
472 prePostProcessTimer_.
stop();
480 std::cout <<
"Linearization/solve/update time: "
487 <<
"\n" << std::flush;
493 prePostProcessTimer_.
start();
495 prePostProcessTimer_.
stop();
500 prePostProcessTimer_.
start();
501 asImp_().succeeded_();
502 prePostProcessTimer_.
stop();
521 if (numIterations_ > targetIterations_()) {
522 Scalar percent = Scalar(numIterations_ - targetIterations_())
523 / targetIterations_();
524 return oldTimeStep / (1.0 + percent);
527 Scalar percent = Scalar(targetIterations_() - numIterations_)
528 / targetIterations_();
529 return oldTimeStep * (1.0 + percent / 1.2);
537 {
return endIterMsgStream_; }
544 { linearSolver_.eraseMatrix(); }
550 {
return linearSolver_; }
556 {
return linearSolver_; }
559 {
return prePostProcessTimer_; }
562 {
return linearizeTimer_; }
565 {
return solveTimer_; }
568 {
return updateTimer_; }
576 return EWOMS_GET_PARAM(TypeTag,
bool, NewtonVerbose) && (comm_.rank() == 0);
585 void begin_(
const SolutionVector& u OPM_UNUSED)
590 convergenceWriter_.beginTimeStep();
606 {
model().linearizer().linearize(); }
608 void preSolve_(
const SolutionVector& currentSolution OPM_UNUSED,
609 const GlobalEqVector& currentResidual)
611 const auto& constraintsMap =
model().linearizer().constraintsMap();
617 for (
unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
619 if (dofIdx >=
model().numGridDof() ||
model().dofTotalVolume(dofIdx) <= 0.0)
623 if (enableConstraints_()) {
624 if (constraintsMap.count(dofIdx) > 0)
628 const auto& r = currentResidual[dofIdx];
629 for (
unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx)
630 error_ = Opm::max(std::abs(r[eqIdx] *
model().eqWeight(dofIdx, eqIdx)), error_);
634 error_ = comm_.max(error_);
639 OPM_THROW(Opm::NumericalProblem,
640 "Newton: Error " << error_
641 <<
" is larger than maximum allowed error of "
658 void postSolve_(
const SolutionVector& nextSolution OPM_UNUSED,
659 const SolutionVector& currentSolution OPM_UNUSED,
660 const GlobalEqVector& currentResidual OPM_UNUSED,
661 const GlobalEqVector& solutionUpdate OPM_UNUSED)
679 const SolutionVector& currentSolution,
680 const GlobalEqVector& solutionUpdate,
681 const GlobalEqVector& currentResidual)
683 const auto& constraintsMap =
model().linearizer().constraintsMap();
687 asImp_().writeConvergence_(currentSolution, solutionUpdate);
690 if (!std::isfinite(solutionUpdate.one_norm()))
691 OPM_THROW(Opm::NumericalProblem,
"Non-finite update!");
693 size_t numGridDof =
model().numGridDof();
694 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
695 if (enableConstraints_()) {
696 if (constraintsMap.count(dofIdx) > 0) {
697 const auto& constraints = constraintsMap.at(dofIdx);
698 asImp_().updateConstraintDof_(dofIdx,
699 nextSolution[dofIdx],
703 asImp_().updatePrimaryVariables_(dofIdx,
704 nextSolution[dofIdx],
705 currentSolution[dofIdx],
706 solutionUpdate[dofIdx],
707 currentResidual[dofIdx]);
710 asImp_().updatePrimaryVariables_(dofIdx,
711 nextSolution[dofIdx],
712 currentSolution[dofIdx],
713 solutionUpdate[dofIdx],
714 currentResidual[dofIdx]);
718 size_t numDof =
model().numTotalDof();
719 for (
size_t dofIdx = numGridDof; dofIdx < numDof; ++dofIdx) {
720 nextSolution[dofIdx] = currentSolution[dofIdx];
721 nextSolution[dofIdx] -= solutionUpdate[dofIdx];
729 PrimaryVariables& nextValue,
730 const Constraints& constraints)
731 { nextValue = constraints; }
737 PrimaryVariables& nextValue,
738 const PrimaryVariables& currentValue,
739 const EqVector& update,
740 const EqVector& currentResidual OPM_UNUSED)
742 nextValue = currentValue;
753 const GlobalEqVector& solutionUpdate)
756 convergenceWriter_.beginIteration();
757 convergenceWriter_.writeFields(currentSolution, solutionUpdate);
758 convergenceWriter_.endIteration();
769 const SolutionVector& currentSolution OPM_UNUSED)
775 std::cout <<
"Newton iteration " << numIterations_ <<
""
776 <<
" error: " << error_
780 endIterMsgStream_.str(
"");
795 else if (asImp_().
numIterations() >= asImp_().maxIterations_()) {
800 return error_ * 4.0 < lastError_;
813 convergenceWriter_.endTimeStep();
822 { numIterations_ = targetIterations_() * 2; }
833 int targetIterations_()
const
836 int maxIterations_()
const
839 static bool enableConstraints_()
842 Simulator& simulator_;
849 std::ostringstream endIterMsgStream_;
859 LinearSolverBackend linearSolver_;
863 CollectiveCommunication comm_;
867 ConvergenceWriter convergenceWriter_;
870 Implementation& asImp_()
871 {
return *
static_cast<Implementation *
>(
this); }
872 const Implementation& asImp_()
const
873 {
return *
static_cast<const Implementation *
>(
this); }
void beginIteration_()
Indicates the beginning of a Newton iteration.
Definition: newtonmethod.hh:596
bool verbose_() const
Returns true if the Newton method ought to be chatty.
Definition: newtonmethod.hh:574
A simple class which makes sure that a timer gets stopped if an exception is thrown.
Definition: timerguard.hh:40
void updatePrimaryVariables_(unsigned globalDofIdx OPM_UNUSED, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector ¤tResidual OPM_UNUSED)
Update a single primary variables object.
Definition: newtonmethod.hh:736
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
LinearSolverBackend & linearSolver()
Returns the linear solver backend object for external use.
Definition: newtonmethod.hh:549
void start()
Start counting the time resources used by the simulation.
Definition: timer.hh:62
void setIterationIndex(int value)
Set the index of current iteration.
Definition: newtonmethod.hh:279
void begin_(const SolutionVector &u OPM_UNUSED)
Called before the Newton method is applied to an non-linear system of equations.
Definition: newtonmethod.hh:585
void succeeded_()
Called if the Newton method was successful.
Definition: newtonmethod.hh:829
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
std::ostringstream & endIterMsg()
Message that should be printed for the user after the end of an iteration.
Definition: newtonmethod.hh:536
void failed_()
Called if the Newton method broke down.
Definition: newtonmethod.hh:821
const LinearSolverBackend & linearSolver() const
Returns the linear solver backend object for external use.
Definition: newtonmethod.hh:555
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
void linearize_()
Linearize the global non-linear system of equations.
Definition: newtonmethod.hh:605
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition: simulator.hh:202
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
void halt()
Stop the measurement reset all timing values.
Definition: timer.hh:99
void writeConvergence_(const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate)
Write the convergence behaviour of the newton method to disk.
Definition: newtonmethod.hh:752
void end_()
Indicates that we're done solving the non-linear system of equations.
Definition: newtonmethod.hh:810
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
bool apply()
Run the Newton method.
Definition: newtonmethod.hh:302
void postSolve_(const SolutionVector &nextSolution OPM_UNUSED, const SolutionVector ¤tSolution OPM_UNUSED, const GlobalEqVector ¤tResidual OPM_UNUSED, const GlobalEqVector &solutionUpdate OPM_UNUSED)
Update the error of the solution given the previous iteration.
Definition: newtonmethod.hh:658
const Problem & problem() const
Returns a reference to the object describing the current physical problem.
Definition: newtonmethod.hh:250
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
This file provides the infrastructure to retrieve run-time parameters.
Scalar tolerance() const
Return the current tolerance at which the Newton method considers itself to be converged.
Definition: newtonmethod.hh:286
A convergence writer for the Newton method which does nothing.
Provides an encapsulation to measure the system time.
Definition: timer.hh:48
void updateConstraintDof_(unsigned globalDofIdx OPM_UNUSED, PrimaryVariables &nextValue, const Constraints &constraints)
Update the primary variables for a degree of freedom which is constraint.
Definition: newtonmethod.hh:728
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
int numIterations() const
Returns the number of iterations done since the Newton method was invoked.
Definition: newtonmethod.hh:269
The multi-dimensional Newton method.
Definition: newtonmethod.hh:56
double realTimeElapsed() const
Return the real time [s] elapsed during the periods the timer was active since the last reset...
Definition: timer.hh:121
static void registerParameters()
Register all run-time parameters for the Newton method.
Definition: newtonmethod.hh:209
const Model & model() const
Returns a reference to the numeric model.
Definition: newtonmethod.hh:262
Provides an encapsulation to measure the system time.
void endIteration_(const SolutionVector &nextSolution OPM_UNUSED, const SolutionVector ¤tSolution OPM_UNUSED)
Indicates that one Newton iteration was finished.
Definition: newtonmethod.hh:768
bool proceed_() const
Returns true iff another Newton iteration should be done.
Definition: newtonmethod.hh:786
Scalar suggestTimeStepSize(Scalar oldTimeStep) const
Suggest a new time-step size based on the old time-step size.
Definition: newtonmethod.hh:515
Provides the magic behind the eWoms property system.
bool converged() const
Returns true if the error of the solution is below the tolerance.
Definition: newtonmethod.hh:238
A simple class which makes sure that a timer gets stopped if an exception is thrown.
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual)
Update the current solution with a delta vector.
Definition: newtonmethod.hh:678
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
Problem & problem()
Returns a reference to the object describing the current physical problem.
Definition: newtonmethod.hh:244
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
Model & model()
Returns a reference to the numeric model.
Definition: newtonmethod.hh:256
void eraseMatrix()
Causes the solve() method to discared the structure of the linear system of equations the next time i...
Definition: newtonmethod.hh:543
A convergence writer for the Newton method which does nothing.
Definition: nullconvergencewriter.hh:51
void setTolerance(Scalar value)
Set the current tolerance at which the Newton method considers itself to be converged.
Definition: newtonmethod.hh:293
double stop()
Stop counting the time resources.
Definition: timer.hh:73
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394