OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-23 09:46:24 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #428 from andlaus/EclNewtonMethod

Browse Source 
ebos: introduce an EclNewtonMethod
This commit is contained in:
Tor Harald Sandve 2019-01-11 13:06:20 +01:00 committed by GitHub
commit f971caeae9
2 changed files with 265 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

232
ebos/eclnewtonmethod.hh Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
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 <ewoms/models/blackoil/blackoilnewtonmethod.hh>
#include <ewoms/common/signum.hh>
#include <opm/material/common/Unused.hpp>
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 TypeTag>
class EclNewtonMethod : public BlackOilNewtonMethod<TypeTag>
{
typedef BlackOilNewtonMethod<TypeTag> 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<TypeTag>;
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<Scalar, numEq> 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

40
ebos/eclproblem.hh
View File

@ -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<TypeTag>);
// Use the dummy gradient calculator in order not to do unnecessary work.
SET_TYPE_PROP(EclBaseProblem, GradientCalculator, Ewoms::EclDummyGradientCalculator<TypeTag>);
// 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<TypeTag>);
// 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