diff --git a/ebos/eclnewtonmethod.hh b/ebos/eclnewtonmethod.hh new file mode 100644 index 000000000..fff1577c7 --- /dev/null +++ b/ebos/eclnewtonmethod.hh @@ -0,0 +1,232 @@ +// -*- 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(EclNewtonSumTolerance); +NEW_PROP_TAG(EclNewtonStrictIterations); +NEW_PROP_TAG(EclNewtonRelaxedVolumeFraction); +NEW_PROP_TAG(EclNewtonSumToleranceExponent); +NEW_PROP_TAG(EclNewtonRelaxedTolerance); + +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; + 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); + } + + /*! + * \brief Register all run-time parameters for the Newton method. + */ + static void registerParameters() + { + ParentType::registerParameters(); + + 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, 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' " + "region during a strict iteration."); + } + + /*! + * \brief Returns true if the error of the solution is below the + * tolerance. + */ + bool converged() const + { + if (errorPvFraction_ < relaxedMaxPvFraction_) + return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ; + else if (this->numIterations() > numStrictIterations_) + 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_); + + // 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 + 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 relaxedMaxPvFraction_; + + Scalar sumTolerance_; + + int numStrictIterations_; +}; +} // namespace Ewoms + +#endif diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index b86a8ec3d..faed62b7b 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,35 @@ 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); +// the default for the allowed volumetric error for oil per second +SET_SCALAR_PROP(EclBaseProblem, NewtonRawTolerance, 1e-2); -// 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); +// the tolerated amount of "incorrect" amount of oil per time step for the complete +// 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-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); + +// 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 @@ -281,6 +303,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