From 07e1b4ecde6d31a022e91db2fa131c3aee8c1954 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 14 Nov 2018 12:13:28 +0100 Subject: [PATCH 1/6] ebos: introduce an EclNewtonMethod this calculates the error and convergence differently from the standard Newton method. --- ebos/eclnewtonmethod.hh | 201 ++++++++++++++++++++++++++++++++++++++++ ebos/eclproblem.hh | 17 ++-- 2 files changed, 210 insertions(+), 8 deletions(-) create mode 100644 ebos/eclnewtonmethod.hh diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh new file mode 100644 index 000000000..c45a3dfa8 --- /dev/null +++ b/ebos/eclnewtonmethod.hh @@ -0,0 +1,201 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Ewoms::EclNewtonMethod + */ +#ifndef EWOMS_ECL_NEWTON_METHOD_HH +#define EWOMS_ECL_NEWTON_METHOD_HH + +#include +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(NewtonSumTolerance); + +END_PROPERTIES + +namespace Ewoms { + +/*! + * \brief A newton solver which is ebos specific. + */ +template +class EclNewtonMethod : public BlackOilNewtonMethod +{ + typedef BlackOilNewtonMethod ParentType; + typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) DiscNewtonMethod; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + static const unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq); + + static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx; + static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx; + static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx; + + friend NewtonMethod; + friend DiscNewtonMethod; + friend ParentType; + +public: + EclNewtonMethod(Simulator& simulator) : ParentType(simulator) + { + errorPvFraction_ = 1.0; + sumTolerance_ = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonSumTolerance); + relaxedTolerance_ = 1e9; + } + + /*! + * \brief Register all run-time parameters for the Newton method. + */ + static void registerParameters() + { + ParentType::registerParameters(); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, NewtonSumTolerance, + "The maximum error tolerated by the Newton" + "method for considering a solution to be " + "converged"); + } + + /*! + * \brief Returns true if the error of the solution is below the + * tolerance. + */ + bool converged() const + { + if (errorPvFraction_ < 0.01) + return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; + else if (this->numIterations() > 8) + return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; + + return this->error_ <= this->tolerance() && errorSum_ <= sumTolerance_; + } + + void preSolve_(const SolutionVector& currentSolution OPM_UNUSED, + const GlobalEqVector& currentResidual) + { + const auto& constraintsMap = this->model().linearizer().constraintsMap(); + this->lastError_ = this->error_; + Scalar newtonMaxError = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxError); + + // calculate the error as the maximum weighted tolerance of + // the solution's residual + this->error_ = 0.0; + Dune::FieldVector componentSumError; + std::fill(componentSumError.begin(), componentSumError.end(), 0.0); + Scalar sumPv = 0.0; + errorPvFraction_ = 0.0; + const Scalar dt = this->simulator_.timeStepSize(); + for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) { + // do not consider auxiliary DOFs for the error + if (dofIdx >= this->model().numGridDof() + || this->model().dofTotalVolume(dofIdx) <= 0.0) + continue; + + if (!this->model().isLocalDof(dofIdx)) + continue; + + // also do not consider DOFs which are constraint + if (this->enableConstraints_()) { + if (constraintsMap.count(dofIdx) > 0) + continue; + } + + const auto& r = currentResidual[dofIdx]; + const double pvValue = + this->simulator_.problem().porosity(dofIdx) + * this->model().dofTotalVolume(dofIdx); + sumPv += pvValue; + bool cnvViolated = false; + + for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) { + Scalar tmpError = r[eqIdx] * dt * this->model().eqWeight(dofIdx, eqIdx) / pvValue; + Scalar tmpError2 = r[eqIdx] * this->model().eqWeight(dofIdx, eqIdx); + + this->error_ = Opm::max(std::abs(tmpError), this->error_); + + if (std::abs(tmpError) > this->tolerance_) + cnvViolated = true; + + componentSumError[eqIdx] += std::abs(tmpError2); + } + if (cnvViolated) + errorPvFraction_ += pvValue; + } + + // take the other processes into account + this->error_ = this->comm_.max(this->error_); + componentSumError = this->comm_.sum(componentSumError); + sumPv = this->comm_.sum(sumPv); + errorPvFraction_ = this->comm_.sum(errorPvFraction_); + + componentSumError /= sumPv; + componentSumError *= dt; + + errorPvFraction_ /= sumPv; + + errorSum_ = 0; + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) + errorSum_ = std::max(std::abs(componentSumError[eqIdx]), errorSum_); + + // make sure that the error never grows beyond the maximum + // allowed one + if (this->error_ > newtonMaxError) + throw Opm::NumericalIssue("Newton: Error "+std::to_string(double(this->error_)) + +" is larger than maximum allowed error of " + +std::to_string(double(newtonMaxError))); + + // make sure that the error never grows beyond the maximum + // allowed one + if (errorSum_ > newtonMaxError) + throw Opm::NumericalIssue("Newton: Sum of the error "+std::to_string(double(errorSum_)) + +" is larger than maximum allowed error of " + +std::to_string(double(newtonMaxError))); + } + +private: + Scalar errorPvFraction_; + Scalar errorSum_; + + Scalar relaxedTolerance_; + + Scalar sumTolerance_; +}; +} // namespace Ewoms + +#endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index b86a8ec3d..c9b8a33d4 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -60,6 +60,7 @@ #include "ecldummygradientcalculator.hh" #include "eclfluxmodule.hh" #include "eclbaseaquifermodel.hh" +#include "eclnewtonmethod.hh" #include "ecltracermodel.hh" #include "vtkecltracermodule.hh" @@ -235,14 +236,10 @@ SET_SCALAR_PROP(EclBaseProblem, EndTime, 1e100); // not millions of trillions of years, that is...) SET_SCALAR_PROP(EclBaseProblem, InitialTimeStepSize, 1e100); -// increase the default raw tolerance for the newton solver to 10^-4 because this is what -// everone else seems to be doing... -SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-4); - -// reduce the maximum allowed Newton error to 0.1 kg/(m^3 s). The rationale is that if -// the error is above that limit, the time step is unlikely to succeed anyway and we can -// thus abort the futile attempt early. -SET_SCALAR_PROP(EclBaseProblem, NewtonMaxError, 0.1); +// increase the default raw tolerance for the newton solver because this is what everone +// else seems to be doing... +SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-2); +SET_SCALAR_PROP(EclBaseProblem, NewtonSumTolerance, 1e-3); // set the maximum number of Newton iterations to 14 because the likelyhood that a time // step succeeds at more than 14 Newton iteration is rather small @@ -281,6 +278,10 @@ SET_TYPE_PROP(EclBaseProblem, FluxModule, Ewoms::EclTransFluxModule); // Use the dummy gradient calculator in order not to do unnecessary work. SET_TYPE_PROP(EclBaseProblem, GradientCalculator, Ewoms::EclDummyGradientCalculator); +// Use a custom Newton-Raphson method class for ebos in order to attain more +// sophisticated update and error computation mechanisms +SET_TYPE_PROP(EclBaseProblem, NewtonMethod, Ewoms::EclNewtonMethod); + // The frequency of writing restart (*.ers) files. This is the number of time steps // between writing restart files SET_INT_PROP(EclBaseProblem, RestartWritingInterval, 0xffffff); // disable From 85a9b75076f52851aacbab441bfb4209a3782ecb Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 14 Nov 2018 12:13:28 +0100 Subject: [PATCH 2/6] ebos: allow larger errors for larger reservoirs albeit, we scale the error only to the cube root of the pore volume. the rationale is that the same amount of mass can get lost "along" a line for each timestep. maybe it would be a good idea to do something like this for time step size as well because taking multiple small time steps currently allows a much larger error in the result than doing it in one big step. --- ebos/eclnewtonmethod.hh | 10 ++++++++++ ebos/eclproblem.hh | 9 ++++++--- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index c45a3dfa8..62d9547e2 100644 --- a/ebos/eclnewtonmethod.hh +++ b/ebos/eclnewtonmethod.hh @@ -173,6 +173,16 @@ public: for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) errorSum_ = std::max(std::abs(componentSumError[eqIdx]), errorSum_); + // update the sum tolerance: larger reservoirs can tolerate a higher amount of + // mass lost per time step than smaller ones! since this is not linear, we use + // the cube root of the overall pore volume, i.e., the value specified by the + // NewtonSumTolerance parameter is the "incorrect" mass per timestep for an + // reservoir that exhibits 1 m^3 of pore volume. A reservoir with a total pore + // volume of 10^3 m^3 will tolerate 10 times as much. + sumTolerance_ = + EWOMS_GET_PARAM(TypeTag, Scalar, NewtonSumTolerance) + * std::cbrt(sumPv); + // make sure that the error never grows beyond the maximum // allowed one if (this->error_ > newtonMaxError) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index c9b8a33d4..ffd2c5876 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -236,10 +236,13 @@ SET_SCALAR_PROP(EclBaseProblem, EndTime, 1e100); // not millions of trillions of years, that is...) SET_SCALAR_PROP(EclBaseProblem, InitialTimeStepSize, 1e100); -// increase the default raw tolerance for the newton solver because this is what everone -// else seems to be doing... +// set the tolerated amount of "incorrect" mass to ~1e-6 kg of oil per time step for a +// reservoir that exhibits a pore volume of 1 m^3. larger reservoirs will tolerate larger +// residuals. +SET_SCALAR_PROP(EclBaseProblem, NewtonSumTolerance, 1e-6); + +// the default for the volumetric error for oil per second is 10^-2 kg/(m^3 * s). SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-2); -SET_SCALAR_PROP(EclBaseProblem, NewtonSumTolerance, 1e-3); // set the maximum number of Newton iterations to 14 because the likelyhood that a time // step succeeds at more than 14 Newton iteration is rather small From c08f0008adaff0986e710cf79a215181692231cc Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 14 Nov 2018 13:05:02 +0100 Subject: [PATCH 3/6] EclNewtonMethod: tweak the parameters a bit --- ebos/eclnewtonmethod.hh | 2 +- ebos/eclproblem.hh | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index 62d9547e2..86931402b 100644 --- a/ebos/eclnewtonmethod.hh +++ b/ebos/eclnewtonmethod.hh @@ -98,7 +98,7 @@ public: */ bool converged() const { - if (errorPvFraction_ < 0.01) + if (errorPvFraction_ < 0.03) return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; else if (this->numIterations() > 8) return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index ffd2c5876..176224e3b 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -236,10 +236,10 @@ SET_SCALAR_PROP(EclBaseProblem, EndTime, 1e100); // not millions of trillions of years, that is...) SET_SCALAR_PROP(EclBaseProblem, InitialTimeStepSize, 1e100); -// set the tolerated amount of "incorrect" mass to ~1e-6 kg of oil per time step for a +// set the tolerated amount of "incorrect" mass to ~1e-4 kg of oil per time step for a // reservoir that exhibits a pore volume of 1 m^3. larger reservoirs will tolerate larger // residuals. -SET_SCALAR_PROP(EclBaseProblem, NewtonSumTolerance, 1e-6); +SET_SCALAR_PROP(EclBaseProblem, NewtonSumTolerance, 1e-4); // the default for the volumetric error for oil per second is 10^-2 kg/(m^3 * s). SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-2); From ec391f529da8791cf8a39f2becce32994e74a2b2 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 7 Jan 2019 13:26:10 +0100 Subject: [PATCH 4/6] make the ebos-Newton specific parameters setable from the command line also, tweak them a bit: increase the sum tolerance before scaling to 1e-3 and reduce the default number of strict iterations to 4. --- ebos/eclnewtonmethod.hh | 35 ++++++++++++++++++++++++++++------- ebos/eclproblem.hh | 26 ++++++++++++++++++++------ 2 files changed, 48 insertions(+), 13 deletions(-) diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index 86931402b..ec40842a9 100644 --- a/ebos/eclnewtonmethod.hh +++ b/ebos/eclnewtonmethod.hh @@ -35,7 +35,10 @@ BEGIN_PROPERTIES -NEW_PROP_TAG(NewtonSumTolerance); +NEW_PROP_TAG(EclNewtonSumTolerance); +NEW_PROP_TAG(EclNewtonStrictIterations); +NEW_PROP_TAG(EclNewtonRelaxedVolumeFraction); +NEW_PROP_TAG(EclNewtonRelaxedTolerance); END_PROPERTIES @@ -75,8 +78,12 @@ public: EclNewtonMethod(Simulator& simulator) : ParentType(simulator) { errorPvFraction_ = 1.0; - sumTolerance_ = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonSumTolerance); - relaxedTolerance_ = 1e9; + relaxedMaxPvFraction_ = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonRelaxedVolumeFraction); + + sumTolerance_ = 0.0; // this gets determined in the error calculation proceedure + relaxedTolerance_ = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonRelaxedTolerance); + + numStrictIterations_ = EWOMS_GET_PARAM(TypeTag, int, EclNewtonStrictIterations); } /*! @@ -86,10 +93,21 @@ public: { ParentType::registerParameters(); - EWOMS_REGISTER_PARAM(TypeTag, Scalar, NewtonSumTolerance, + EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonSumTolerance, "The maximum error tolerated by the Newton" "method for considering a solution to be " "converged"); + EWOMS_REGISTER_PARAM(TypeTag, int, EclNewtonStrictIterations, + "The number of Newton iterations where the" + " volumetric error is considered."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonRelaxedVolumeFraction, + "The fraction of the pore volume of the reservoir " + "where the volumetric error may be voilated during " + "strict Newton iterations."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonRelaxedTolerance, + "The maximum error which the volumetric residual " + "may exhibit if it is in a 'relaxed' " + "region during a strict iteration."); } /*! @@ -98,9 +116,9 @@ public: */ bool converged() const { - if (errorPvFraction_ < 0.03) + if (errorPvFraction_ < relaxedMaxPvFraction_) return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; - else if (this->numIterations() > 8) + else if (this->numIterations() > numStrictIterations_) return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; return this->error_ <= this->tolerance() && errorSum_ <= sumTolerance_; @@ -180,7 +198,7 @@ public: // reservoir that exhibits 1 m^3 of pore volume. A reservoir with a total pore // volume of 10^3 m^3 will tolerate 10 times as much. sumTolerance_ = - EWOMS_GET_PARAM(TypeTag, Scalar, NewtonSumTolerance) + EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumTolerance) * std::cbrt(sumPv); // make sure that the error never grows beyond the maximum @@ -203,8 +221,11 @@ private: Scalar errorSum_; Scalar relaxedTolerance_; + Scalar relaxedMaxPvFraction_; Scalar sumTolerance_; + + int numStrictIterations_; }; } // namespace Ewoms diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 176224e3b..195871199 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -236,14 +236,28 @@ SET_SCALAR_PROP(EclBaseProblem, EndTime, 1e100); // not millions of trillions of years, that is...) SET_SCALAR_PROP(EclBaseProblem, InitialTimeStepSize, 1e100); -// set the tolerated amount of "incorrect" mass to ~1e-4 kg of oil per time step for a -// reservoir that exhibits a pore volume of 1 m^3. larger reservoirs will tolerate larger -// residuals. -SET_SCALAR_PROP(EclBaseProblem, NewtonSumTolerance, 1e-4); - -// the default for the volumetric error for oil per second is 10^-2 kg/(m^3 * s). +// the default for the allowed volumetric error for oil per second SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-2); +// the tolerated amount of "incorrect" amount of oil per time step for the complete +// reservoi. this is scaled by the pore volume of the reservoir, i.e., larger reservoirs +// will tolerate larger residuals. +SET_SCALAR_PROP(EclBaseProblem, EclNewtonSumTolerance, 1e-3); + +// set number of Newton iterations where the volumetric residual is considered for +// convergence +SET_INT_PROP(EclBaseProblem, EclNewtonStrictIterations, 4); + +// set fraction of the pore volume where the volumetric residual may be violated during +// strict Newton iterations +SET_SCALAR_PROP(EclBaseProblem, EclNewtonRelaxedVolumeFraction, 0.03); + +// the maximum volumetric error of a cell in the relaxed region +SET_SCALAR_PROP(EclBaseProblem, EclNewtonRelaxedTolerance, 1e9); + +// Ignore the maximum error mass for early termination of the newton method. +SET_SCALAR_PROP(EclBaseProblem, NewtonMaxError, 10e9); + // set the maximum number of Newton iterations to 14 because the likelyhood that a time // step succeeds at more than 14 Newton iteration is rather small SET_INT_PROP(EclBaseProblem, NewtonMaxIterations, 14); From 0492796e852c7bb87f8fc89a2ce5b5cbf12db891 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 11 Jan 2019 11:05:18 +0100 Subject: [PATCH 5/6] address review comments --- ebos/eclproblem.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 195871199..091061a15 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -240,13 +240,13 @@ SET_SCALAR_PROP(EclBaseProblem, InitialTimeStepSize, 1e100); SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-2); // the tolerated amount of "incorrect" amount of oil per time step for the complete -// reservoi. this is scaled by the pore volume of the reservoir, i.e., larger reservoirs +// reservoir. this is scaled by the pore volume of the reservoir, i.e., larger reservoirs // will tolerate larger residuals. -SET_SCALAR_PROP(EclBaseProblem, EclNewtonSumTolerance, 1e-3); +SET_SCALAR_PROP(EclBaseProblem, EclNewtonSumTolerance, 1e-4); // set number of Newton iterations where the volumetric residual is considered for // convergence -SET_INT_PROP(EclBaseProblem, EclNewtonStrictIterations, 4); +SET_INT_PROP(EclBaseProblem, EclNewtonStrictIterations, 8); // set fraction of the pore volume where the volumetric residual may be violated during // strict Newton iterations From c2377e7b1070b754fd52fbc1cb6d9af95f389c78 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Fri, 11 Jan 2019 11:18:06 +0100 Subject: [PATCH 6/6] EclNewtonMethod: make the exponent of the tolerance scaling settable by a parameter the parameter is called `EclNewtonSumToleranceExponent`. if it is set to 1, the specified tolerance will be used directly. (this is not desireable in the general case though, because at the same result quality, the sum error for large reservoirs can be larger than for small ones.) --- ebos/eclnewtonmethod.hh | 18 +++++++++--------- ebos/eclproblem.hh | 8 ++++++++ 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh index ec40842a9..fff1577c7 100644 --- a/ebos/eclnewtonmethod.hh +++ b/ebos/eclnewtonmethod.hh @@ -38,6 +38,7 @@ BEGIN_PROPERTIES NEW_PROP_TAG(EclNewtonSumTolerance); NEW_PROP_TAG(EclNewtonStrictIterations); NEW_PROP_TAG(EclNewtonRelaxedVolumeFraction); +NEW_PROP_TAG(EclNewtonSumToleranceExponent); NEW_PROP_TAG(EclNewtonRelaxedTolerance); END_PROPERTIES @@ -104,6 +105,9 @@ public: "The fraction of the pore volume of the reservoir " "where the volumetric error may be voilated during " "strict Newton iterations."); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonSumToleranceExponent, + "The the exponent used to scale the sum tolerance by " + "the total pore volume of the reservoir."); EWOMS_REGISTER_PARAM(TypeTag, Scalar, EclNewtonRelaxedTolerance, "The maximum error which the volumetric residual " "may exhibit if it is in a 'relaxed' " @@ -191,15 +195,11 @@ public: for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) errorSum_ = std::max(std::abs(componentSumError[eqIdx]), errorSum_); - // update the sum tolerance: larger reservoirs can tolerate a higher amount of - // mass lost per time step than smaller ones! since this is not linear, we use - // the cube root of the overall pore volume, i.e., the value specified by the - // NewtonSumTolerance parameter is the "incorrect" mass per timestep for an - // reservoir that exhibits 1 m^3 of pore volume. A reservoir with a total pore - // volume of 10^3 m^3 will tolerate 10 times as much. - sumTolerance_ = - EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumTolerance) - * std::cbrt(sumPv); + // scale the tolerance for the total error with the pore volume. by default, the + // exponent is 1/3, i.e., cubic root. + Scalar x = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumTolerance); + Scalar y = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumToleranceExponent); + sumTolerance_ = x*std::pow(sumPv, y); // make sure that the error never grows beyond the maximum // allowed one diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 091061a15..faed62b7b 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -244,6 +244,14 @@ SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-2); // will tolerate larger residuals. SET_SCALAR_PROP(EclBaseProblem, EclNewtonSumTolerance, 1e-4); +// set the exponent for the volume scaling of the sum tolerance: larger reservoirs can +// tolerate a higher amount of mass lost per time step than smaller ones! since this is +// not linear, we use the cube root of the overall pore volume by default, i.e., the +// value specified by the NewtonSumTolerance parameter is the "incorrect" mass per +// timestep for an reservoir that exhibits 1 m^3 of pore volume. A reservoir with a total +// pore volume of 10^3 m^3 will tolerate 10 times as much. +SET_SCALAR_PROP(EclBaseProblem, EclNewtonSumToleranceExponent, 1.0/3.0); + // set number of Newton iterations where the volumetric residual is considered for // convergence SET_INT_PROP(EclBaseProblem, EclNewtonStrictIterations, 8);