Loading...
Searching...
No Matches
richardsnewtonmethod.hh
Go to the documentation of this file.
// -*- 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.
*/
#ifndef EWOMS_RICHARDS_NEWTON_METHOD_HH
#define EWOMS_RICHARDS_NEWTON_METHOD_HH
#include <dune/common/fvector.hh>
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/discretization/common/fvbasenewtonmethod.hh>
#include <opm/models/richards/richardsproperties.hh>
#include <algorithm>
namespace Opm {
template <class TypeTag>
class RichardsNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMethod>
{
using ParentType = GetPropType<TypeTag, Properties::DiscNewtonMethod>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using EqVector = GetPropType<TypeTag, Properties::EqVector>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Linearizer = GetPropType<TypeTag, Properties::Linearizer>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
enum { pressureWIdx = Indices::pressureWIdx };
enum { numPhases = FluidSystem::numPhases };
enum { liquidPhaseIdx = getPropValue<TypeTag, Properties::LiquidPhaseIndex>() };
enum { gasPhaseIdx = getPropValue<TypeTag, Properties::GasPhaseIndex>() };
using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
public:
explicit RichardsNewtonMethod(Simulator& simulator)
: ParentType(simulator)
{}
protected:
friend NewtonMethod<TypeTag>;
friend ParentType;
void updatePrimaryVariables_(unsigned globalDofIdx,
PrimaryVariables& nextValue,
const PrimaryVariables& currentValue,
const EqVector& update,
const EqVector&)
{
// normal Newton-Raphson update
nextValue = currentValue;
nextValue -= update;
// do not clamp anything after 4 iterations
if (this->numIterations_ > 4) {
return;
}
const auto& problem = this->simulator_.problem();
// calculate the old wetting phase saturation
const MaterialLawParams& matParams =
problem.materialLawParams(globalDofIdx, /*timeIdx=*/0);
ImmiscibleFluidState<Scalar, FluidSystem> fs;
// set the temperature
const Scalar T = problem.temperature(globalDofIdx, /*timeIdx=*/0);
fs.setTemperature(T);
// calculate the phase pressures of the previous iteration
// first, we have to find the minimum capillary pressure
// (i.e. Sw = 0)
fs.setSaturation(liquidPhaseIdx, 1.0);
fs.setSaturation(gasPhaseIdx, 0.0);
PhaseVector pC;
MaterialLaw::capillaryPressures(pC, matParams, fs);
// non-wetting pressure can be larger than the
// reference pressure if the medium is fully
// saturated by the wetting phase
const Scalar pWOld = currentValue[pressureWIdx];
const Scalar pNOld =
std::max(problem.referencePressure(globalDofIdx, /*timeIdx=*/0),
pWOld + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
// find the saturations of the previous iteration
fs.setPressure(liquidPhaseIdx, pWOld);
fs.setPressure(gasPhaseIdx, pNOld);
PhaseVector satOld;
MaterialLaw::saturations(satOld, matParams, fs);
satOld[liquidPhaseIdx] = std::max(Scalar{0.0}, satOld[liquidPhaseIdx]);
// find the wetting phase pressures which
// corrospond to a 20% increase and a 20% decrease
// of the wetting saturation
fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] - 0.2);
fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] - 0.2));
MaterialLaw::capillaryPressures(pC, matParams, fs);
const Scalar pwMin = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] + 0.2);
fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] + 0.2));
MaterialLaw::capillaryPressures(pC, matParams, fs);
const Scalar pwMax = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
// clamp the result to the minimum and the maximum
// pressures we just calculated
const Scalar pW = std::clamp(nextValue[pressureWIdx], pwMin, pwMax);
nextValue[pressureWIdx] = pW;
}
};
} // namespace Opm
#endif
Opm::NewtonMethod
The multi-dimensional Newton method.
Definition newtonmethod.hh:100
Opm::RichardsNewtonMethod::updatePrimaryVariables_
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables &currentValue, const EqVector &update, const EqVector &)
Update a single primary variables object.
Definition richardsnewtonmethod.hh:82
fvbasenewtonmethod.hh
A Newton method for models using a finite volume discretization.
multiphasebaseproperties.hh
Defines the common properties required by the porous medium multi-phase models.
Opm
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilbioeffectsmodules.hh:43
Opm::GetPropType
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:233
Opm::getPropValue
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:240
richardsproperties.hh
Contains the property declarations for the Richards model.