From 7c21a630e5a219feaa825a56638f45dc7230ccbe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 2 Oct 2015 09:19:07 +0200 Subject: [PATCH 1/4] Rename NewtonSolver -> NonlinearSolver. --- CMakeLists_files.cmake | 4 +- .../{NewtonSolver.hpp => NonlinearSolver.hpp} | 47 ++++++------- ...lver_impl.hpp => NonlinearSolver_impl.hpp} | 69 +++++++++++-------- opm/autodiff/SimulatorBase.hpp | 1 - opm/autodiff/SimulatorBase_impl.hpp | 8 +-- .../SimulatorFullyImplicitBlackoil.hpp | 9 ++- .../SimulatorFullyImplicitBlackoilSolvent.hpp | 4 +- 7 files changed, 78 insertions(+), 64 deletions(-) rename opm/autodiff/{NewtonSolver.hpp => NonlinearSolver.hpp} (71%) rename opm/autodiff/{NewtonSolver_impl.hpp => NonlinearSolver_impl.hpp} (76%) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index b7c9af3c2..332a36f58 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -126,8 +126,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/NewtonIterationBlackoilInterleaved.hpp opm/autodiff/NewtonIterationBlackoilSimple.hpp opm/autodiff/NewtonIterationUtilities.hpp - opm/autodiff/NewtonSolver.hpp - opm/autodiff/NewtonSolver_impl.hpp + opm/autodiff/NonlinearSolver.hpp + opm/autodiff/NonlinearSolver_impl.hpp opm/autodiff/LinearisedBlackoilResidual.hpp opm/autodiff/ParallelDebugOutput.hpp opm/autodiff/RateConverter.hpp diff --git a/opm/autodiff/NewtonSolver.hpp b/opm/autodiff/NonlinearSolver.hpp similarity index 71% rename from opm/autodiff/NewtonSolver.hpp rename to opm/autodiff/NonlinearSolver.hpp index dfdb89697..0b7ee7319 100644 --- a/opm/autodiff/NewtonSolver.hpp +++ b/opm/autodiff/NonlinearSolver.hpp @@ -18,8 +18,8 @@ along with OPM. If not, see . */ -#ifndef OPM_NEWTONSOLVER_HEADER_INCLUDED -#define OPM_NEWTONSOLVER_HEADER_INCLUDED +#ifndef OPM_NONLINEARSOLVER_HEADER_INCLUDED +#define OPM_NONLINEARSOLVER_HEADER_INCLUDED #include #include @@ -28,9 +28,10 @@ namespace Opm { - /// A Newton solver class suitable for general fully-implicit models. + /// A nonlinear solver class suitable for general fully-implicit models, + /// as well as pressure, transport and sequential models. template - class NewtonSolver + class NonlinearSolver { public: // --------- Types and enums --------- @@ -38,18 +39,18 @@ namespace Opm { typedef ADB::V V; typedef ADB::M M; - // The Newton relaxation scheme type + // Available relaxation scheme types. enum RelaxType { DAMPEN, SOR }; - // Solver parameters controlling nonlinear Newton process. + // Solver parameters controlling nonlinear process. struct SolverParameters { enum RelaxType relax_type_; double relax_max_; double relax_increment_; double relax_rel_tol_; - int max_iter_; // max newton iterations - int min_iter_; // min newton iterations + int max_iter_; // max nonlinear iterations + int min_iter_; // min nonlinear iterations explicit SolverParameters( const parameter::ParameterGroup& param ); SolverParameters(); @@ -66,11 +67,11 @@ namespace Opm { /// Construct solver for a given model. /// /// The model is a std::unique_ptr because the object to which model points to is - /// not allowed to be deleted as long as the NewtonSolver object exists. + /// not allowed to be deleted as long as the NonlinearSolver object exists. /// - /// \param[in] param parameters controlling nonlinear Newton process + /// \param[in] param parameters controlling nonlinear process /// \param[in, out] model physical simulation model. - explicit NewtonSolver(const SolverParameters& param, + explicit NonlinearSolver(const SolverParameters& param, std::unique_ptr model); /// Take a single forward step, after which the states will be modified @@ -84,14 +85,14 @@ namespace Opm { ReservoirState& reservoir_state, WellState& well_state); - /// Number of Newton iterations used in all calls to step(). - unsigned int newtonIterations() const; + /// Number of nonlinear solver iterations used in all calls to step(). + unsigned int nonlinearIterations() const; /// Number of linear solver iterations used in all calls to step(). unsigned int linearIterations() const; - /// Number of linear solver iterations used in the last call to step(). - unsigned int newtonIterationsLastStep() const; + /// Number of nonlinear solver iterations used in the last call to step(). + unsigned int nonlinearIterationsLastStep() const; /// Number of linear solver iterations used in the last call to step(). unsigned int linearIterationsLastStep() const; @@ -103,9 +104,9 @@ namespace Opm { // --------- Data members --------- SolverParameters param_; std::unique_ptr model_; - unsigned int newtonIterations_; + unsigned int nonlinearIterations_; unsigned int linearIterations_; - unsigned int newtonIterationsLast_; + unsigned int nonlinearIterationsLast_; unsigned int linearIterationsLast_; // --------- Private methods --------- @@ -115,13 +116,13 @@ namespace Opm { double relaxRelTol() const { return param_.relax_rel_tol_; } double maxIter() const { return param_.max_iter_; } double minIter() const { return param_.min_iter_; } - void detectNewtonOscillations(const std::vector>& residual_history, - const int it, const double relaxRelTol, - bool& oscillate, bool& stagnate) const; - void stabilizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const; + void detectOscillations(const std::vector>& residual_history, + const int it, const double relaxRelTol, + bool& oscillate, bool& stagnate) const; + void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const; }; } // namespace Opm -#include "NewtonSolver_impl.hpp" +#include "NonlinearSolver_impl.hpp" -#endif // OPM_NEWTONSOLVER_HEADER_INCLUDED +#endif // OPM_NONLINEARSOLVER_HEADER_INCLUDED diff --git a/opm/autodiff/NewtonSolver_impl.hpp b/opm/autodiff/NonlinearSolver_impl.hpp similarity index 76% rename from opm/autodiff/NewtonSolver_impl.hpp rename to opm/autodiff/NonlinearSolver_impl.hpp index 095f7f8db..00b2e20ee 100644 --- a/opm/autodiff/NewtonSolver_impl.hpp +++ b/opm/autodiff/NonlinearSolver_impl.hpp @@ -20,31 +20,33 @@ along with OPM. If not, see . */ -#ifndef OPM_NEWTONSOLVER_IMPL_HEADER_INCLUDED -#define OPM_NEWTONSOLVER_IMPL_HEADER_INCLUDED +#ifndef OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED +#define OPM_NONLINEARSOLVER_IMPL_HEADER_INCLUDED -#include +#include namespace Opm { template - NewtonSolver::NewtonSolver(const SolverParameters& param, - std::unique_ptr model) + NonlinearSolver::NonlinearSolver(const SolverParameters& param, + std::unique_ptr model) : param_(param), model_(std::move(model)), - newtonIterations_(0), - linearIterations_(0) + nonlinearIterations_(0), + linearIterations_(0), + nonlinearIterationsLast_(0), + linearIterationsLast_(0) { } template - unsigned int NewtonSolver::newtonIterations () const + unsigned int NonlinearSolver::nonlinearIterations() const { - return newtonIterations_; + return nonlinearIterations_; } template - unsigned int NewtonSolver::linearIterations () const + unsigned int NonlinearSolver::linearIterations() const { return linearIterations_; } @@ -56,9 +58,22 @@ namespace Opm return *model_; } + template + unsigned int NonlinearSolver::nonlinearIterationsLastStep() const + { + return nonlinearIterationsLast_; + } + + template + unsigned int NonlinearSolver::linearIterationsLastStep() const + { + return linearIterationsLast_; + } + + template int - NewtonSolver:: + NonlinearSolver:: step(const double dt, ReservoirState& reservoir_state, WellState& well_state) @@ -74,7 +89,7 @@ namespace Opm model_->assemble(reservoir_state, well_state, true); residual_norms_history.push_back(model_->computeResidualNorms()); - // Set up for main Newton loop. + // Set up for main solver loop. double omega = 1.0; int iteration = 0; bool converged = model_->getConvergence(dt, iteration); @@ -85,16 +100,16 @@ namespace Opm const enum RelaxType relaxtype = relaxType(); int linIters = 0; - // ---------- Main Newton loop ---------- + // ---------- Main nonlinear solver loop ---------- while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) { - // Compute the Newton update to the primary variables. + // Compute the update to the primary variables. V dx = model_->solveJacobianSystem(); // Store number of linear iterations used. linIters += model_->linearIterationsLastSolve(); - // Stabilize the Newton update. - detectNewtonOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate); + // Stabilize the nonlinear update. + detectOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate); if (isOscillate) { omega -= relaxIncrement(); omega = std::max(omega, relaxMax()); @@ -102,7 +117,7 @@ namespace Opm std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl; } } - stabilizeNewton(dx, dxOld, omega, relaxtype); + stabilizeNonlinearUpdate(dx, dxOld, omega, relaxtype); // Apply the update, the model may apply model-dependent // limitations and chopping of the update. @@ -126,9 +141,9 @@ namespace Opm } linearIterations_ += linIters; - newtonIterations_ += iteration; + nonlinearIterations_ += iteration; linearIterationsLast_ = linIters; - newtonIterationsLast_ = iteration; + nonlinearIterationsLast_ = iteration; // Do model-specific post-step actions. model_->afterStep(dt, reservoir_state, well_state); @@ -139,7 +154,7 @@ namespace Opm template - void NewtonSolver::SolverParameters:: + void NonlinearSolver::SolverParameters:: reset() { // default values for the solver parameters @@ -152,7 +167,7 @@ namespace Opm } template - NewtonSolver::SolverParameters:: + NonlinearSolver::SolverParameters:: SolverParameters() { // set default values @@ -160,7 +175,7 @@ namespace Opm } template - NewtonSolver::SolverParameters:: + NonlinearSolver::SolverParameters:: SolverParameters( const parameter::ParameterGroup& param ) { // set default values @@ -183,9 +198,9 @@ namespace Opm template void - NewtonSolver::detectNewtonOscillations(const std::vector>& residual_history, - const int it, const double relaxRelTol_arg, - bool& oscillate, bool& stagnate) const + NonlinearSolver::detectOscillations(const std::vector>& residual_history, + const int it, const double relaxRelTol_arg, + bool& oscillate, bool& stagnate) const { // The detection of oscillation in two primary variable results in the report of the detection // of oscillation for the solver. @@ -220,8 +235,8 @@ namespace Opm template void - NewtonSolver::stabilizeNewton(V& dx, V& dxOld, const double omega, - const RelaxType relax_type) const + NonlinearSolver::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega, + const RelaxType relax_type) const { // The dxOld is updated with dx. // If omega is equal to 1., no relaxtion will be appiled. diff --git a/opm/autodiff/SimulatorBase.hpp b/opm/autodiff/SimulatorBase.hpp index 80069e74c..5e22d9326 100644 --- a/opm/autodiff/SimulatorBase.hpp +++ b/opm/autodiff/SimulatorBase.hpp @@ -26,7 +26,6 @@ #include #include -#include #include #include #include diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index 3154f73ed..72be25bf5 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -106,7 +106,7 @@ namespace Opm output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep ); } - unsigned int totalNewtonIterations = 0; + unsigned int totalNonlinearIterations = 0; unsigned int totalLinearIterations = 0; // Main simulation loop. @@ -167,8 +167,8 @@ namespace Opm // take time that was used to solve system for this reportStep solver_timer.stop(); - // accumulate the number of Newton and Linear Iterations - totalNewtonIterations += solver->newtonIterations(); + // accumulate the number of nonlinear and linear Iterations + totalNonlinearIterations += solver->nonlinearIterations(); totalLinearIterations += solver->linearIterations(); // Report timing. @@ -201,7 +201,7 @@ namespace Opm report.pressure_time = stime; report.transport_time = 0.0; report.total_time = total_timer.secsSinceStart(); - report.total_newton_iterations = totalNewtonIterations; + report.total_newton_iterations = totalNonlinearIterations; report.total_linear_iterations = totalLinearIterations; return report; } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index dca177e5f..ef6080741 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp @@ -1,5 +1,5 @@ /* - Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. Copyright 2015 Andreas Lauser This file is part of the Open Porous Media project (OPM). @@ -21,9 +21,8 @@ #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED -#include "SimulatorBase.hpp" - -#include "NewtonSolver.hpp" +#include +#include namespace Opm { @@ -38,7 +37,7 @@ struct SimulatorTraits > typedef BlackoilOutputWriter OutputWriter; typedef GridT Grid; typedef BlackoilModel Model; - typedef NewtonSolver Solver; + typedef NonlinearSolver Solver; }; /// a simulator for the blackoil model diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp index 1919bdaf4..823e11485 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp @@ -32,7 +32,7 @@ #include #include #include -#include +#include #include #include @@ -88,7 +88,7 @@ namespace Opm typedef BlackoilOutputWriter OutputWriter; typedef GridT Grid; typedef BlackoilSolventModel Model; - typedef NewtonSolver Solver; + typedef NonlinearSolver Solver; }; /// Class collecting all necessary components for a blackoil simulation with polymer From 29a1a891d2053e489a00b253e37c66edd9c03da7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 2 Oct 2015 13:51:40 +0200 Subject: [PATCH 2/4] Make nonlinearIteration() the only interface to the model. This means that details such as calling assemble(), solveJacobianSystem(), updateState() etc. are now left to the model class. This will make it easier to create new model classes with different behaviour (such as sequential models). --- opm/autodiff/BlackoilModelBase.hpp | 30 +++++++++++ opm/autodiff/BlackoilModelBase_impl.hpp | 51 +++++++++++++++++++ opm/autodiff/NonlinearSolver.hpp | 21 +++++--- opm/autodiff/NonlinearSolver_impl.hpp | 67 +++++++------------------ 4 files changed, 113 insertions(+), 56 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 785ded71e..61d32942d 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -76,6 +76,17 @@ namespace Opm { + /// Class used for reporting the outcome of a nonlinearIteration() call. + struct IterationReport + { + bool failed; + bool converged; + int linear_iterations; + }; + + + + /// Traits to encapsulate the types used by classes using or /// extending this model. Forward declared here, must be /// specialised for each concrete model class. @@ -155,6 +166,22 @@ namespace Opm { ReservoirState& reservoir_state, WellState& well_state); + /// Called once per nonlinear iteration. + /// This model will perform a Newton-Raphson update, changing reservoir_state + /// and well_state. It will also use the nonlinear_solver to do relaxation of + /// updates if necessary. + /// \param[in] iteration should be 0 for the first call of a new timestep + /// \param[in] dt time step size + /// \param[in] nonlinear_solver nonlinear solver used (for oscillation/relaxation control) + /// \param[in, out] reservoir_state reservoir state variables + /// \param[in, out] well_state well state variables + template + IterationReport nonlinearIteration(const int iteration, + const double dt, + NonlinearSolverType& nonlinear_solver, + ReservoirState& reservoir_state, + WellState& well_state); + /// Called once after each time step. /// In this class, this function does nothing. /// \param[in] dt time step size @@ -290,6 +317,9 @@ namespace Opm { std::vector primalVariable_; V pvdt_; std::vector material_name_; + std::vector> residual_norms_history_; + double current_relaxation_; + V dx_old_; // --------- Protected methods --------- diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 97073ec7a..4d3790cf7 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -180,6 +180,7 @@ namespace detail { { 1.1169, 1.0031, 0.0031 }} ) // the default magic numbers , terminal_output_ (terminal_output) , material_name_{ "Water", "Oil", "Gas" } + , current_relaxation_(1.0) { assert(numMaterials() == 3); // Due to the material_name_ init above. #if HAVE_MPI @@ -227,6 +228,56 @@ namespace detail { + template + template + IterationReport + BlackoilModelBase:: + nonlinearIteration(const int iteration, + const double dt, + NonlinearSolverType& nonlinear_solver, + ReservoirState& reservoir_state, + WellState& well_state) + { + if (iteration == 0) { + // For each iteration we store in a vector the norms of the residual of + // the mass balance for each active phase, the well flux and the well equations. + residual_norms_history_.clear(); + current_relaxation_ = 1.0; + dx_old_ = V::Zero(sizeNonLinear()); + } + assemble(reservoir_state, well_state, iteration == 0); + residual_norms_history_.push_back(computeResidualNorms()); + const bool converged = getConvergence(dt, iteration); + if (!converged) { + V dx = solveJacobianSystem(); + + // Stabilize the nonlinear update. + bool isOscillate = false; + bool isStagnate = false; + nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate); + if (isOscillate) { + current_relaxation_ -= nonlinear_solver.relaxIncrement(); + current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax()); + if (terminalOutputEnabled()) { + std::cout << " Oscillating behavior detected: Relaxation set to " + << current_relaxation_ << std::endl; + } + } + nonlinear_solver.stabilizeNonlinearUpdate(dx, dx_old_, current_relaxation_); + + // Apply the update, applying model-dependent + // limitations and chopping of the update. + updateState(dx, reservoir_state, well_state); + } + const bool failed = false; // Not needed in this model. + const int linear_iters = converged ? 0 : linearIterationsLastSolve(); + return IterationReport{ failed, converged, linear_iters }; + } + + + + + template void BlackoilModelBase:: diff --git a/opm/autodiff/NonlinearSolver.hpp b/opm/autodiff/NonlinearSolver.hpp index 0b7ee7319..134de8fc9 100644 --- a/opm/autodiff/NonlinearSolver.hpp +++ b/opm/autodiff/NonlinearSolver.hpp @@ -97,9 +97,22 @@ namespace Opm { /// Number of linear solver iterations used in the last call to step(). unsigned int linearIterationsLastStep() const; - /// return reference to physical model + /// Reference to physical model. const PhysicalModel& model() const; + /// Detect oscillation or stagnation in a given residual history. + void detectOscillations(const std::vector>& residual_history, + const int it, bool& oscillate, bool& stagnate) const; + + /// Apply a stabilization to dx, depending on dxOld and relaxation parameters. + void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const; + + /// The greatest relaxation factor (i.e. smallest factor) allowed. + double relaxMax() const { return param_.relax_max_; } + + /// The step-change size for the relaxation factor. + double relaxIncrement() const { return param_.relax_increment_; } + private: // --------- Data members --------- SolverParameters param_; @@ -111,15 +124,9 @@ namespace Opm { // --------- Private methods --------- enum RelaxType relaxType() const { return param_.relax_type_; } - double relaxMax() const { return param_.relax_max_; } - double relaxIncrement() const { return param_.relax_increment_; } double relaxRelTol() const { return param_.relax_rel_tol_; } double maxIter() const { return param_.max_iter_; } double minIter() const { return param_.min_iter_; } - void detectOscillations(const std::vector>& residual_history, - const int it, const double relaxRelTol, - bool& oscillate, bool& stagnate) const; - void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const; }; } // namespace Opm diff --git a/opm/autodiff/NonlinearSolver_impl.hpp b/opm/autodiff/NonlinearSolver_impl.hpp index 00b2e20ee..bd74a86ed 100644 --- a/opm/autodiff/NonlinearSolver_impl.hpp +++ b/opm/autodiff/NonlinearSolver_impl.hpp @@ -81,57 +81,27 @@ namespace Opm // Do model-specific once-per-step calculations. model_->prepareStep(dt, reservoir_state, well_state); - // For each iteration we store in a vector the norms of the residual of - // the mass balance for each active phase, the well flux and the well equations. - std::vector> residual_norms_history; + int iteration = 0; - // Assemble residual and Jacobian, store residual norms. - model_->assemble(reservoir_state, well_state, true); - residual_norms_history.push_back(model_->computeResidualNorms()); + // Let the model do one nonlinear iteration. // Set up for main solver loop. - double omega = 1.0; - int iteration = 0; - bool converged = model_->getConvergence(dt, iteration); - const int sizeNonLinear = model_->sizeNonLinear(); - V dxOld = V::Zero(sizeNonLinear); - bool isOscillate = false; - bool isStagnate = false; - const enum RelaxType relaxtype = relaxType(); int linIters = 0; + bool converged = false; // ---------- Main nonlinear solver loop ---------- - while ( (!converged && (iteration < maxIter())) || (minIter() > iteration)) { - // Compute the update to the primary variables. - V dx = model_->solveJacobianSystem(); - - // Store number of linear iterations used. - linIters += model_->linearIterationsLastSolve(); - - // Stabilize the nonlinear update. - detectOscillations(residual_norms_history, iteration, relaxRelTol(), isOscillate, isStagnate); - if (isOscillate) { - omega -= relaxIncrement(); - omega = std::max(omega, relaxMax()); - if (model_->terminalOutputEnabled()) { - std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl; - } + do { + IterationReport report = model_->nonlinearIteration(iteration, dt, *this, reservoir_state, well_state); + if (report.failed) { + OPM_THROW(Opm::NumericalProblem, "Failed to complete a nonlinear iteration."); } - stabilizeNonlinearUpdate(dx, dxOld, omega, relaxtype); - - // Apply the update, the model may apply model-dependent - // limitations and chopping of the update. - model_->updateState(dx, reservoir_state, well_state); - - // Assemble residual and Jacobian, store residual norms. - model_->assemble(reservoir_state, well_state, false); - residual_norms_history.push_back(model_->computeResidualNorms()); - - // increase iteration counter + if (report.converged) { + assert(report.linear_iterations == 0); + converged = true; + } + linIters += report.linear_iterations; ++iteration; - - converged = model_->getConvergence(dt, iteration); - } + } while ( (!converged && (iteration <= maxIter())) || (minIter() > iteration)); if (!converged) { if (model_->terminalOutputEnabled()) { @@ -141,7 +111,7 @@ namespace Opm } linearIterations_ += linIters; - nonlinearIterations_ += iteration; + nonlinearIterations_ += iteration - 1; // Since the last one will always be trivial. linearIterationsLast_ = linIters; nonlinearIterationsLast_ = iteration; @@ -199,7 +169,7 @@ namespace Opm template void NonlinearSolver::detectOscillations(const std::vector>& residual_history, - const int it, const double relaxRelTol_arg, + const int it, bool& oscillate, bool& stagnate) const { // The detection of oscillation in two primary variable results in the report of the detection @@ -222,7 +192,7 @@ namespace Opm const double d1 = std::abs((F0[p] - F2[p]) / F0[p]); const double d2 = std::abs((F0[p] - F1[p]) / F0[p]); - oscillatePhase += (d1 < relaxRelTol_arg) && (relaxRelTol_arg < d2); + oscillatePhase += (d1 < relaxRelTol()) && (relaxRelTol() < d2); // Process is 'stagnate' unless at least one phase // exhibits significant residual change. @@ -235,8 +205,7 @@ namespace Opm template void - NonlinearSolver::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega, - const RelaxType relax_type) const + NonlinearSolver::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const { // The dxOld is updated with dx. // If omega is equal to 1., no relaxtion will be appiled. @@ -244,7 +213,7 @@ namespace Opm const V tempDxOld = dxOld; dxOld = dx; - switch (relax_type) { + switch (relaxType()) { case DAMPEN: if (omega == 1.) { return; From 87330984aa1fc6a1d4525a122b30d7cfa2780386 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 12 Nov 2015 17:55:13 +0100 Subject: [PATCH 3/4] Adapt to renaming of class. --- opm/autodiff/NonlinearSolver_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/NonlinearSolver_impl.hpp b/opm/autodiff/NonlinearSolver_impl.hpp index bd74a86ed..37c6db3c0 100644 --- a/opm/autodiff/NonlinearSolver_impl.hpp +++ b/opm/autodiff/NonlinearSolver_impl.hpp @@ -52,7 +52,7 @@ namespace Opm } template - const PhysicalModel& NewtonSolver::model() const + const PhysicalModel& NonlinearSolver::model() const { assert( model_ ); return *model_; From 88fbe5e9b6d47d5797bc33106c148abcab068f0c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 13 Nov 2015 15:08:00 +0100 Subject: [PATCH 4/4] Ensure non-null model in constructor. --- opm/autodiff/NonlinearSolver_impl.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/NonlinearSolver_impl.hpp b/opm/autodiff/NonlinearSolver_impl.hpp index 37c6db3c0..984b54d42 100644 --- a/opm/autodiff/NonlinearSolver_impl.hpp +++ b/opm/autodiff/NonlinearSolver_impl.hpp @@ -37,6 +37,9 @@ namespace Opm nonlinearIterationsLast_(0), linearIterationsLast_(0) { + if (!model_) { + OPM_THROW(std::logic_error, "Must provide a non-null model argument for NonlinearSolver."); + } } template @@ -54,7 +57,6 @@ namespace Opm template const PhysicalModel& NonlinearSolver::model() const { - assert( model_ ); return *model_; }