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; }
217 const char *clearRemainingLine =
"\n";
218 if (isatty(fileno(stdout))) {
219 static const char blubb[] = { 0x1b,
'[',
'K',
'\r', 0 };
220 clearRemainingLine = blubb;
224 prePostProcessTimer_.
halt();
225 linearizeTimer_.
halt();
229 SolutionVector& nextSolution =
model().solution(0);
230 SolutionVector currentSolution(nextSolution);
231 GlobalEqVector solutionUpdate(nextSolution.size());
233 Linearizer& linearizer =
model().linearizer();
235 TimerGuard prePostProcessTimerGuard(prePostProcessTimer_);
238 prePostProcessTimer_.
start();
239 asImp_().begin_(nextSolution);
240 prePostProcessTimer_.
stop();
243 TimerGuard innerPrePostProcessTimerGuard(prePostProcessTimer_);
244 TimerGuard linearizeTimerGuard(linearizeTimer_);
255 prePostProcessTimer_.
start();
256 asImp_().beginIteration_();
257 prePostProcessTimer_.
stop();
260 currentSolution = nextSolution;
263 std::cout <<
"Linearize: r(x^k) = dS/dt + div F - q; M = grad r"
264 << clearRemainingLine
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();
286 asImp_().preSolve_(currentSolution, residual);
290 if (asImp_().
verbose_() && isatty(fileno(stdout)))
291 std::cout << clearRemainingLine
295 prePostProcessTimer_.
start();
296 asImp_().endIteration_(nextSolution, currentSolution);
297 prePostProcessTimer_.
stop();
304 std::cout <<
"Solve: M deltax^k = r"
305 << clearRemainingLine
312 linearSolver_.setMatrix(jacobian);
313 solutionUpdate = 0.0;
314 bool converged = linearSolver_.solve(solutionUpdate);
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"
332 << clearRemainingLine
338 updateTimer_.
start();
339 asImp_().postSolve_(currentSolution,
342 asImp_().update_(nextSolution, currentSolution, solutionUpdate, residual);
345 if (asImp_().
verbose_() && isatty(fileno(stdout)))
347 std::cout << clearRemainingLine
351 prePostProcessTimer_.
start();
352 asImp_().endIteration_(nextSolution, currentSolution);
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();
368 catch (
const NumericalProblem& e)
371 std::cout <<
"Newton method caught exception: \""
372 << e.what() <<
"\"\n" << std::flush;
374 prePostProcessTimer_.
start();
376 prePostProcessTimer_.
stop();
382 if (asImp_().
verbose_() && isatty(fileno(stdout)))
383 std::cout << clearRemainingLine
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_;
440 Scalar nextDt = std::max(
problem().minTimeStepSize(),
441 oldDt / (Scalar{1.0} + percent));
445 Scalar percent = Scalar(params_.targetIterations_ - numIterations_) / params_.targetIterations_;
446 Scalar nextDt = std::max(
problem().minTimeStepSize(),
447 oldDt*(Scalar{1.0} + percent / Scalar{1.2}));
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();
521 bool succeeded =
true;
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;
533 succeeded = comm.min(succeeded);
536 throw NumericalProblem(
"pre processing of the problem failed");
547 model().linearizer().linearizeDomain();
550 void linearizeAuxiliaryEquations_()
552 model().linearizer().linearizeAuxiliaryEquations();
553 model().linearizer().finalize();
556 void preSolve_(
const SolutionVector&,
557 const GlobalEqVector& currentResidual)
559 const auto& constraintsMap =
model().linearizer().constraintsMap();
561 Scalar newtonMaxError = params_.maxError_;
566 for (
unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
568 if (dofIdx >=
model().numGridDof() ||
model().dofTotalVolume(dofIdx) <= 0.0)
572 if (enableConstraints_()) {
573 if (constraintsMap.count(dofIdx) > 0)
577 const auto& r = currentResidual[dofIdx];
578 for (
unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx)
579 error_ = max(std::abs(r[eqIdx] *
model().eqWeight(dofIdx, eqIdx)), error_);
583 error_ = comm_.max(error_);
587 if (error_ > newtonMaxError)
588 throw NumericalProblem(
"Newton: Error "+std::to_string(
double(error_))
589 +
" is larger than maximum allowed error of "
590 + std::to_string(
double(newtonMaxError)));
606 const GlobalEqVector&,
607 GlobalEqVector& solutionUpdate)
611 auto&
model = simulator_.model();
612 const auto& comm = simulator_.gridView().comm();
613 for (
unsigned i = 0; i <
model.numAuxiliaryModules(); ++i) {
614 auto& auxMod = *
model.auxiliaryModule(i);
616 bool succeeded =
true;
618 auxMod.postSolve(solutionUpdate);
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;
628 succeeded = comm.min(succeeded);
631 throw NumericalProblem(
"post processing of an auxilary equation failed");
650 const SolutionVector& currentSolution,
651 const GlobalEqVector& solutionUpdate,
652 const GlobalEqVector& currentResidual)
654 const auto& constraintsMap =
model().linearizer().constraintsMap();
658 asImp_().writeConvergence_(currentSolution, solutionUpdate);
661 if (!std::isfinite(solutionUpdate.one_norm()))
662 throw NumericalProblem(
"Non-finite update!");
664 size_t numGridDof =
model().numGridDof();
665 for (
unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
666 if (enableConstraints_()) {
667 if (constraintsMap.count(dofIdx) > 0) {
668 const auto& constraints = constraintsMap.at(dofIdx);
669 asImp_().updateConstraintDof_(dofIdx,
670 nextSolution[dofIdx],
674 asImp_().updatePrimaryVariables_(dofIdx,
675 nextSolution[dofIdx],
676 currentSolution[dofIdx],
677 solutionUpdate[dofIdx],
678 currentResidual[dofIdx]);
681 asImp_().updatePrimaryVariables_(dofIdx,
682 nextSolution[dofIdx],
683 currentSolution[dofIdx],
684 solutionUpdate[dofIdx],
685 currentResidual[dofIdx]);
689 size_t numDof =
model().numTotalDof();
690 for (
size_t dofIdx = numGridDof; dofIdx < numDof; ++dofIdx) {
691 nextSolution[dofIdx] = currentSolution[dofIdx];
692 nextSolution[dofIdx] -= solutionUpdate[dofIdx];
700 PrimaryVariables& nextValue,
701 const Constraints& constraints)
702 { nextValue = constraints; }
708 PrimaryVariables& nextValue,
709 const PrimaryVariables& currentValue,
710 const EqVector& update,
713 nextValue = currentValue;
724 const GlobalEqVector& solutionUpdate)
726 if (params_.writeConvergence_) {
727 convergenceWriter_.beginIteration();
728 convergenceWriter_.writeFields(currentSolution, solutionUpdate);
729 convergenceWriter_.endIteration();
740 const SolutionVector&)
744 const auto& comm = simulator_.gridView().comm();
745 bool succeeded =
true;
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;
757 succeeded = comm.min(succeeded);
760 throw NumericalProblem(
"post processing of the problem failed");
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_;
833 NewtonMethodParams<Scalar> params_;
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
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
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
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 nullconvergencewriter.hh:39
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.