From 07e1b4ecde6d31a022e91db2fa131c3aee8c1954 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 14 Nov 2018 12:13:28 +0100 Subject: [PATCH] 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