27#ifndef EWOMS_NEWTON_METHOD_HH
28#define EWOMS_NEWTON_METHOD_HH
30#include <dune/istl/istlexception.hh>
31#include <dune/common/classname.hh>
32#include <dune/common/parallel/mpihelper.hh>
34#include <opm/common/Exceptions.hpp>
36#include <opm/material/densead/Math.hpp>
40#include <opm/models/nonlinear/newtonmethodparams.hpp>
41#include <opm/models/nonlinear/newtonmethodproperties.hh>
56template <
class TypeTag>
64namespace Opm::Properties {
75template<
class TypeTag>
77template<
class TypeTag>
90template <
class TypeTag>
108 using Communicator =
typename Dune::MPIHelper::MPICommunicator;
109 using CollectiveCommunication =
typename Dune::Communication<typename Dune::MPIHelper::MPICommunicator>;
113 : simulator_(simulator)
114 , endIterMsgStream_(std::ostringstream::out)
115 , linearSolver_(simulator)
116 , comm_(Dune::MPIHelper::getCommunicator())
117 , convergenceWriter_(asImp_())
131 LinearSolverBackend::registerParameters();
155 {
return simulator_.problem(); }
161 {
return simulator_.problem(); }
167 {
return simulator_.model(); }
173 {
return simulator_.model(); }
180 {
return numIterations_; }
190 { numIterations_ = value; }
197 {
return params_.tolerance_; }
204 { params_.tolerance_ = value; }
219 static const char blubb[] = { 0x1b,
'[',
'K',
'\r', 0 };
224 prePostProcessTimer_.
halt();
225 linearizeTimer_.
halt();
233 Linearizer& linearizer =
model().linearizer();
238 prePostProcessTimer_.
start();
240 prePostProcessTimer_.
stop();
255 prePostProcessTimer_.
start();
256 asImp_().beginIteration_();
257 prePostProcessTimer_.
stop();
263 std::cout <<
"Linearize: r(x^k) = dS/dt + div F - q; M = grad r"
269 linearizeTimer_.
start();
270 asImp_().linearizeDomain_();
271 asImp_().linearizeAuxiliaryEquations_();
272 linearizeTimer_.
stop();
275 auto& residual = linearizer.residual();
276 const auto& jacobian = linearizer.jacobian();
277 linearSolver_.prepare(jacobian, residual);
278 linearSolver_.setResidual(residual);
279 linearSolver_.getResidual(residual);
285 updateTimer_.
start();
295 prePostProcessTimer_.
start();
297 prePostProcessTimer_.
stop();
304 std::cout <<
"Solve: M deltax^k = r"
312 linearSolver_.setMatrix(jacobian);
320 std::cout <<
"Newton: Linear solver did not converge\n" << std::flush;
322 prePostProcessTimer_.
start();
324 prePostProcessTimer_.
stop();
331 std::cout <<
"Update: x^(k+1) = x^k - deltax^k"
338 updateTimer_.
start();
351 prePostProcessTimer_.
start();
353 prePostProcessTimer_.
stop();
356 catch (
const Dune::Exception&
e)
359 std::cout <<
"Newton method caught exception: \""
360 <<
e.what() <<
"\"\n" << std::flush;
362 prePostProcessTimer_.
start();
364 prePostProcessTimer_.
stop();
371 std::cout <<
"Newton method caught exception: \""
372 <<
e.what() <<
"\"\n" << std::flush;
374 prePostProcessTimer_.
start();
376 prePostProcessTimer_.
stop();
387 prePostProcessTimer_.
start();
389 prePostProcessTimer_.
stop();
397 std::cout <<
"Linearization/solve/update time: "
404 <<
"\n" << std::flush;
410 prePostProcessTimer_.
start();
412 prePostProcessTimer_.
stop();
417 prePostProcessTimer_.
start();
418 asImp_().succeeded_();
419 prePostProcessTimer_.
stop();
438 if (numIterations_ > params_.targetIterations_) {
439 Scalar
percent = Scalar(numIterations_ - params_.targetIterations_) / params_.targetIterations_;
445 Scalar
percent = Scalar(params_.targetIterations_ - numIterations_) / params_.targetIterations_;
456 {
return endIterMsgStream_; }
463 { linearSolver_.eraseMatrix(); }
469 {
return linearSolver_; }
475 {
return linearSolver_; }
477 const Timer& prePostProcessTimer()
const
478 {
return prePostProcessTimer_; }
480 const Timer& linearizeTimer()
const
481 {
return linearizeTimer_; }
483 const Timer& solveTimer()
const
484 {
return solveTimer_; }
486 const Timer& updateTimer()
const
487 {
return updateTimer_; }
495 return params_.verbose_ && (comm_.rank() == 0);
508 if (params_.writeConvergence_) {
509 convergenceWriter_.beginTimeStep();
519 endIterMsgStream_.str(
"");
520 const auto& comm = simulator_.gridView().comm();
525 catch (
const std::exception&
e) {
528 std::cout <<
"rank " << simulator_.gridView().comm().rank()
529 <<
" caught an exception while pre-processing the problem:" <<
e.what()
530 <<
"\n" << std::flush;
547 model().linearizer().linearizeDomain();
550 void linearizeAuxiliaryEquations_()
552 model().linearizer().linearizeAuxiliaryEquations();
553 model().linearizer().finalize();
556 void preSolve_(
const SolutionVector&,
559 const auto& constraintsMap =
model().linearizer().constraintsMap();
572 if (enableConstraints_()) {
573 if (constraintsMap.count(
dofIdx) > 0)
583 error_ = comm_.max(error_);
589 +
" is larger than maximum allowed error of "
606 const GlobalEqVector&,
611 auto&
model = simulator_.model();
612 const auto& comm = simulator_.gridView().comm();
613 for (
unsigned i = 0; i <
model.numAuxiliaryModules(); ++i) {
620 catch (
const std::exception&
e) {
623 std::cout <<
"rank " << simulator_.gridView().comm().rank()
624 <<
" caught an exception while post processing an auxiliary module:" <<
e.what()
625 <<
"\n" << std::flush;
654 const auto& constraintsMap =
model().linearizer().constraintsMap();
664 size_t numGridDof =
model().numGridDof();
666 if (enableConstraints_()) {
667 if (constraintsMap.count(
dofIdx) > 0) {
668 const auto& constraints = constraintsMap.at(
dofIdx);
669 asImp_().updateConstraintDof_(
dofIdx,
674 asImp_().updatePrimaryVariables_(
dofIdx,
681 asImp_().updatePrimaryVariables_(
dofIdx,
689 size_t numDof =
model().numTotalDof();
701 const Constraints& constraints)
710 const EqVector& update,
726 if (params_.writeConvergence_) {
727 convergenceWriter_.beginIteration();
729 convergenceWriter_.endIteration();
740 const SolutionVector&)
744 const auto& comm = simulator_.gridView().comm();
749 catch (
const std::exception&
e) {
752 std::cout <<
"rank " << simulator_.gridView().comm().rank()
753 <<
" caught an exception while letting the problem post-process:" <<
e.what()
754 <<
"\n" << std::flush;
763 std::cout <<
"Newton iteration " << numIterations_ <<
""
764 <<
" error: " << error_
781 else if (asImp_().
numIterations() >= params_.maxIterations_) {
786 return error_ * 4.0 < lastError_;
798 if (params_.writeConvergence_) {
799 convergenceWriter_.endTimeStep();
809 { numIterations_ = params_.targetIterations_ * 2; }
819 static bool enableConstraints_()
822 Simulator& simulator_;
824 Timer prePostProcessTimer_;
825 Timer linearizeTimer_;
829 std::ostringstream endIterMsgStream_;
839 LinearSolverBackend linearSolver_;
843 CollectiveCommunication comm_;
847 ConvergenceWriter convergenceWriter_;
850 Implementation& asImp_()
851 {
return *
static_cast<Implementation *
>(
this); }
852 const Implementation& asImp_()
const
853 {
return *
static_cast<const Implementation *
>(
this); }
The multi-dimensional Newton method.
Definition newtonmethod.hh:92
bool verbose_() const
Returns true if the Newton method ought to be chatty.
Definition newtonmethod.hh:493
void end_()
Indicates that we're done solving the non-linear system of equations.
Definition newtonmethod.hh:796
bool proceed_() const
Returns true iff another Newton iteration should be done.
Definition newtonmethod.hh:772
int numIterations() const
Returns the number of iterations done since the Newton method was invoked.
Definition newtonmethod.hh:179
const LinearSolverBackend & linearSolver() const
Returns the linear solver backend object for external use.
Definition newtonmethod.hh:474
void updatePrimaryVariables_(unsigned, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector &)
Update a single primary variables object.
Definition newtonmethod.hh:707
void setTolerance(Scalar value)
Set the current tolerance at which the Newton method considers itself to be converged.
Definition newtonmethod.hh:203
std::ostringstream & endIterMsg()
Message that should be printed for the user after the end of an iteration.
Definition newtonmethod.hh:455
static void registerParameters()
Register all run-time parameters for the Newton method.
Definition newtonmethod.hh:129
void postSolve_(const SolutionVector &, const GlobalEqVector &, GlobalEqVector &solutionUpdate)
Update the error of the solution given the previous iteration.
Definition newtonmethod.hh:605
bool apply()
Run the Newton method.
Definition newtonmethod.hh:212
void writeConvergence_(const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate)
Write the convergence behaviour of the newton method to disk.
Definition newtonmethod.hh:723
void begin_(const SolutionVector &)
Called before the Newton method is applied to an non-linear system of equations.
Definition newtonmethod.hh:504
void failed_()
Called if the Newton method broke down.
Definition newtonmethod.hh:808
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual)
Update the current solution with a delta vector.
Definition newtonmethod.hh:649
const Model & model() const
Returns a reference to the numeric model.
Definition newtonmethod.hh:172
void succeeded_()
Called if the Newton method was successful.
Definition newtonmethod.hh:816
Scalar tolerance() const
Return the current tolerance at which the Newton method considers itself to be converged.
Definition newtonmethod.hh:196
LinearSolverBackend & linearSolver()
Returns the linear solver backend object for external use.
Definition newtonmethod.hh:468
void linearizeDomain_()
Linearize the global non-linear system of equations associated with the spatial domain.
Definition newtonmethod.hh:545
const Problem & problem() const
Returns a reference to the object describing the current physical problem.
Definition newtonmethod.hh:160
Scalar suggestTimeStepSize(Scalar oldDt) const
Suggest a new time-step size based on the old time-step size.
Definition newtonmethod.hh:432
bool converged() const
Returns true if the error of the solution is below the tolerance.
Definition newtonmethod.hh:148
Problem & problem()
Returns a reference to the object describing the current physical problem.
Definition newtonmethod.hh:154
void eraseMatrix()
Causes the solve() method to discared the structure of the linear system of equations the next time i...
Definition newtonmethod.hh:462
void updateConstraintDof_(unsigned, PrimaryVariables &nextValue, const Constraints &constraints)
Update the primary variables for a degree of freedom which is constraint.
Definition newtonmethod.hh:699
void finishInit()
Finialize the construction of the object.
Definition newtonmethod.hh:141
Model & model()
Returns a reference to the numeric model.
Definition newtonmethod.hh:166
void beginIteration_()
Indicates the beginning of a Newton iteration.
Definition newtonmethod.hh:516
void endIteration_(const SolutionVector &, const SolutionVector &)
Indicates that one Newton iteration was finished.
Definition newtonmethod.hh:739
void setIterationIndex(int value)
Set the index of current iteration.
Definition newtonmethod.hh:189
A convergence writer for the Newton method which does nothing.
Definition nullconvergencewriter.hh:51
A simple class which makes sure that a timer gets stopped if an exception is thrown.
Definition timerguard.hh:41
Provides an encapsulation to measure the system time.
Definition timer.hpp:46
void start()
Start counting the time resources used by the simulation.
Definition timer.cpp:46
void halt()
Stop the measurement reset all timing values.
Definition timer.cpp:75
double realTimeElapsed() const
Return the real time [s] elapsed during the periods the timer was active since the last reset.
Definition timer.cpp:90
double stop()
Stop counting the time resources.
Definition timer.cpp:52
Declare the properties used by the infrastructure code of the finite volume discretizations.
Declares the properties required by the black oil model.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
A convergence writer for the Newton method which does nothing.
static void registerParameters()
Registers the parameters in parameter system.
Definition newtonmethodparams.cpp:36
Specifies the type of the class which writes out the Newton convergence.
Definition newtonmethodproperties.hh:40
Specifies the type of the actual Newton method.
Definition newtonmethodproperties.hh:32
The type tag on which the default properties for the Newton method are attached.
Definition newtonmethod.hh:70
Provides an encapsulation to measure the system time.
A simple class which makes sure that a timer gets stopped if an exception is thrown.