Loading...
Searching...
No Matches
fvbaselocalresidual.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_FV_BASE_LOCAL_RESIDUAL_HH
#define EWOMS_FV_BASE_LOCAL_RESIDUAL_HH
#include <dune/common/classname.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/geometry.hh>
#include <dune/istl/bvector.hh>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/utils/alignedallocator.hh>
#include <opm/models/utils/parametersystem.hpp>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <type_traits>
namespace Opm {
template<class TypeTag>
class FvBaseLocalResidual
{
private:
using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using EqVector = GetPropType<TypeTag, Properties::EqVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using BoundaryContext = GetPropType<TypeTag, Properties::BoundaryContext>;
static constexpr bool useVolumetricResidual = getPropValue<TypeTag, Properties::UseVolumetricResidual>();
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { extensiveStorageTerm = getPropValue<TypeTag, Properties::ExtensiveStorageTerm>() };
using Toolbox = MathToolbox<Evaluation>;
using EvalVector = Dune::FieldVector<Evaluation, numEq>;
public:
using LocalEvalBlockVector = Dune::BlockVector<EvalVector, aligned_allocator<EvalVector, alignof(EvalVector)>>;
FvBaseLocalResidual() = default;
// copying the local residual class is not a good idea
FvBaseLocalResidual(const FvBaseLocalResidual&) = delete;
static void registerParameters()
{}
const LocalEvalBlockVector& residual() const
{ return internalResidual_; }
const EvalVector& residual(unsigned dofIdx) const
{ return internalResidual_[dofIdx]; }
void eval(const Problem& problem, const Element& element)
{
ElementContext elemCtx(problem);
elemCtx.updateAll(element);
eval(elemCtx);
}
void eval(ElementContext& elemCtx)
{
std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
internalResidual_.resize(numDof);
asImp_().eval(internalResidual_, elemCtx);
}
void eval(LocalEvalBlockVector& residual,
ElementContext& elemCtx) const
{
assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));
residual = 0.0;
// evaluate the flux terms
asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
// evaluate the storage and the source terms
asImp_().evalVolumeTerms_(residual, elemCtx);
// evaluate the boundary conditions
asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
if (useVolumetricResidual) {
// make the residual volume specific (i.e., make it incorrect mass per cubic
// meter instead of total mass)
const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
if (elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0) > 0.0) {
// interior DOF
const Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0);
assert(std::isfinite(dofVolume));
Valgrind::CheckDefined(dofVolume);
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
residual[dofIdx][eqIdx] /= dofVolume;
}
}
}
}
}
void evalStorage(LocalEvalBlockVector& storage,
const ElementContext& elemCtx,
unsigned timeIdx) const
{
// the derivative of the storage term depends on the current primary variables;
// for time indices != 0, the storage term is constant (because these solutions
// are not changed by the Newton method!)
if (timeIdx == 0) {
// calculate the amount of conservation each quantity inside
// all primary sub control volumes
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
storage[dofIdx] = 0.0;
// the volume of the associated DOF
const Scalar alpha =
elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume() *
elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
// If the degree of freedom which we currently look at is the one at the
// center of attention, we need to consider the derivatives for the
// storage term, else the storage term is constant w.r.t. the primary
// variables of the focused DOF.
if (dofIdx == elemCtx.focusDofIndex()) {
asImp_().computeStorage(storage[dofIdx],
elemCtx,
dofIdx,
timeIdx);
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
storage[dofIdx][eqIdx] *= alpha;
}
}
else {
Dune::FieldVector<Scalar, numEq> tmp;
asImp_().computeStorage(tmp,
elemCtx,
dofIdx,
timeIdx);
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
storage[dofIdx][eqIdx] = tmp[eqIdx] * alpha;
}
}
}
}
else {
// for all previous solutions, the storage term does _not_ depend on the
// current primary variables, so we use scalars to store it.
if (elemCtx.enableStorageCache()) {
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx, timeIdx);
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
storage[dofIdx][eqIdx] = cachedStorage[eqIdx];
}
}
}
else {
// calculate the amount of conservation each quantity inside
// all primary sub control volumes
Dune::FieldVector<Scalar, numEq> tmp;
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
tmp = 0.0;
asImp_().computeStorage(tmp,
elemCtx,
dofIdx,
timeIdx);
tmp *= elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume() *
elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
storage[dofIdx][eqIdx] = tmp[eqIdx];
}
}
}
#ifndef NDEBUG
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
Valgrind::CheckDefined(storage[dofIdx][eqIdx]);
assert(isfinite(storage[dofIdx][eqIdx]));
}
}
#endif
}
void evalFluxes(LocalEvalBlockVector& residual,
const ElementContext& elemCtx,
unsigned timeIdx) const
{
RateVector flux;
const auto& stencil = elemCtx.stencil(timeIdx);
// calculate the mass flux over the sub-control volume faces
const std::size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; ++scvfIdx) {
const auto& face = stencil.interiorFace(scvfIdx);
const unsigned i = face.interiorIndex();
const unsigned j = face.exteriorIndex();
Valgrind::SetUndefined(flux);
asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx);
Valgrind::CheckDefined(flux);
#ifndef NDEBUG
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
assert(isfinite(flux[eqIdx]));
}
#endif
const Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor() * face.area();
Valgrind::CheckDefined(alpha);
assert(alpha > 0.0);
assert(isfinite(alpha));
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
flux[eqIdx] *= alpha;
}
// The balance equation for a finite volume is given by
//
// dStorage/dt + Flux = Source
//
// where the 'Flux' and the 'Source' terms represent the
// mass per second which leaves the finite
// volume. Re-arranging this, we get
//
// dStorage/dt + Flux - Source = 0
//
// Since the mass flux as calculated by computeFlux() goes out of sub-control
// volume i and into sub-control volume j, we need to add the flux to finite
// volume i and subtract it from finite volume j
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
assert(isfinite(flux[eqIdx]));
residual[i][eqIdx] += flux[eqIdx];
residual[j][eqIdx] -= flux[eqIdx];
}
}
#if !defined NDEBUG
// in debug mode, ensure that the residual is well-defined
const std::size_t numDof = elemCtx.numDof(timeIdx);
for (unsigned i = 0; i < numDof; ++i) {
for (unsigned j = 0; j < numEq; ++j) {
assert(isfinite(residual[i][j]));
Valgrind::CheckDefined(residual[i][j]);
}
}
#endif
}
// The following methods _must_ be overloaded by the actual flow
// models!
void computeStorage(EqVector&,
const ElementContext&,
unsigned,
unsigned) const
{
throw std::logic_error("Not implemented: The local residual " + Dune::className<Implementation>() +
" does not implement the required method 'computeStorage()'");
}
void computeFlux(RateVector&,
const ElementContext&,
unsigned,
unsigned) const
{
throw std::logic_error("Not implemented: The local residual " + Dune::className<Implementation>() +
" does not implement the required method 'computeFlux()'");
}
void computeSource(RateVector&,
const ElementContext&,
unsigned,
unsigned) const
{
throw std::logic_error("Not implemented: The local residual " + Dune::className<Implementation>() +
" does not implement the required method 'computeSource()'");
}
protected:
void evalBoundary_(LocalEvalBlockVector& residual,
const ElementContext& elemCtx,
unsigned timeIdx) const
{
if (!elemCtx.onBoundary()) {
return;
}
BoundaryContext boundaryCtx(elemCtx);
// move the iterator to the first boundary
if (boundaryCtx.intersection(0).neighbor()) {
boundaryCtx.increment();
}
// evaluate the boundary for all boundary faces of the current context
const std::size_t numBoundaryFaces = boundaryCtx.numBoundaryFaces(/*timeIdx=*/0);
for (unsigned faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx, boundaryCtx.increment()) {
// add the residual of all vertices of the boundary
// segment
evalBoundarySegment_(residual,
boundaryCtx,
faceIdx,
timeIdx);
}
#if !defined NDEBUG
// in debug mode, ensure that the residual and the storage terms are well-defined
const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
for (unsigned i = 0; i < numDof; ++i) {
for (unsigned j = 0; j < numEq; ++j) {
assert(isfinite(residual[i][j]));
Valgrind::CheckDefined(residual[i][j]);
}
}
#endif
}
void evalBoundarySegment_(LocalEvalBlockVector& residual,
const BoundaryContext& boundaryCtx,
unsigned boundaryFaceIdx,
unsigned timeIdx) const
{
BoundaryRateVector values;
Valgrind::SetUndefined(values);
boundaryCtx.problem().boundary(values, boundaryCtx, boundaryFaceIdx, timeIdx);
Valgrind::CheckDefined(values);
const auto& stencil = boundaryCtx.stencil(timeIdx);
const unsigned dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex();
const auto& insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx);
for (unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx) {
values[eqIdx] *= stencil.boundaryFace(boundaryFaceIdx).area() *
insideIntQuants.extrusionFactor();
Valgrind::CheckDefined(values[eqIdx]);
assert(isfinite(values[eqIdx]));
}
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
residual[dofIdx][eqIdx] += values[eqIdx];
}
}
void evalVolumeTerms_(LocalEvalBlockVector& residual,
ElementContext& elemCtx) const
{
EvalVector tmp;
EqVector tmp2;
RateVector sourceRate;
tmp = 0.0;
tmp2 = 0.0;
// evaluate the volumetric terms (storage + source terms)
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
const Scalar extrusionFactor =
elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).extrusionFactor();
Valgrind::CheckDefined(extrusionFactor);
assert(isfinite(extrusionFactor));
assert(extrusionFactor > 0.0);
const Scalar scvVolume =
elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume() * extrusionFactor;
Valgrind::CheckDefined(scvVolume);
assert(isfinite(scvVolume));
assert(scvVolume > 0.0);
// if the model uses extensive quantities in its storage term, and we use
// automatic differention and current DOF is also not the one we currently
// focus on, the storage term does not need any derivatives!
if (!extensiveStorageTerm &&
!std::is_same_v<Scalar, Evaluation> &&
dofIdx != elemCtx.focusDofIndex())
{
asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/0);
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
tmp[eqIdx] = tmp2[eqIdx];
}
}
else {
asImp_().computeStorage(tmp, elemCtx, dofIdx, /*timeIdx=*/0);
}
#ifndef NDEBUG
Valgrind::CheckDefined(tmp);
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
assert(isfinite(tmp[eqIdx]));
}
#endif
if (elemCtx.enableStorageCache()) {
const auto& model = elemCtx.model();
const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
if (model.newtonMethod().numIterations() == 0 &&
!elemCtx.haveStashedIntensiveQuantities())
{
if (!elemCtx.problem().recycleFirstIterationStorage()) {
// we re-calculate the storage term for the solution of the
// previous time step from scratch instead of using the one of
// the first iteration of the current time step.
tmp2 = 0.0;
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/1);
asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1);
}
else {
// if the storage term is cached and we're in the first iteration
// of the time step, use the storage term of the first iteration
// as the one as the solution of the last time step (this assumes
// that the initial guess for the solution at the end of the time
// step is the same as the solution at the beginning of the time
// step. This is usually true, but some fancy preprocessing
// scheme might invalidate that assumption.)
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]);
}
}
Valgrind::CheckDefined(tmp2);
model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2);
}
else {
// if the mass storage at the beginning of the time step is not cached,
// if the storage term is cached and we're not looking at the first
// iteration of the time step, we take the cached data.
tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1);
Valgrind::CheckDefined(tmp2);
}
}
else {
// if the mass storage at the beginning of the time step is not cached,
// we re-calculate it from scratch.
tmp2 = 0.0;
asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1);
Valgrind::CheckDefined(tmp2);
}
// Use the implicit Euler time discretization
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
const double dt = elemCtx.simulator().timeStepSize();
assert(dt > 0);
tmp[eqIdx] -= tmp2[eqIdx];
tmp[eqIdx] *= scvVolume / dt;
residual[dofIdx][eqIdx] += tmp[eqIdx];
}
Valgrind::CheckDefined(residual[dofIdx]);
// deal with the source term
asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
// if the model uses extensive quantities in its storage term, and we use
// automatic differention and current DOF is also not the one we currently
// focus on, the storage term does not need any derivatives!
if (!extensiveStorageTerm &&
!std::is_same_v<Scalar, Evaluation> &&
dofIdx != elemCtx.focusDofIndex())
{
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
residual[dofIdx][eqIdx] -= scalarValue(sourceRate[eqIdx])*scvVolume;
}
}
else {
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
sourceRate[eqIdx] *= scvVolume;
residual[dofIdx][eqIdx] -= sourceRate[eqIdx];
}
}
Valgrind::CheckDefined(residual[dofIdx]);
}
#if !defined NDEBUG
// in debug mode, ensure that the residual is well-defined
std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
for (unsigned i = 0; i < numDof; ++i) {
for (unsigned j = 0; j < numEq; ++j) {
assert(isfinite(residual[i][j]));
Valgrind::CheckDefined(residual[i][j]);
}
}
#endif
}
private:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation*>(this); }
LocalEvalBlockVector internalResidual_;
};
} // namespace Opm
#endif
alignedallocator.hh
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1....
Opm::FvBaseLocalResidual::evalVolumeTerms_
void evalVolumeTerms_(LocalEvalBlockVector &residual, ElementContext &elemCtx) const
Add the change in the storage terms and the source term to the local residual of all sub-control volu...
Definition fvbaselocalresidual.hh:495
Opm::FvBaseLocalResidual::registerParameters
static void registerParameters()
Register all run-time parameters for the local residual.
Definition fvbaselocalresidual.hh:99
Opm::FvBaseLocalResidual::computeSource
void computeSource(RateVector &, const ElementContext &, unsigned, unsigned) const
Calculate the source term of the equation.
Definition fvbaselocalresidual.hh:409
Opm::FvBaseLocalResidual::evalFluxes
void evalFluxes(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Add the flux term to a local residual.
Definition fvbaselocalresidual.hh:302
Opm::FvBaseLocalResidual::computeStorage
void computeStorage(EqVector &, const ElementContext &, unsigned, unsigned) const
Evaluate the amount all conservation quantities (e.g.
Definition fvbaselocalresidual.hh:378
Opm::FvBaseLocalResidual::evalStorage
void evalStorage(LocalEvalBlockVector &storage, const ElementContext &elemCtx, unsigned timeIdx) const
Calculate the amount of all conservation quantities stored in all element's sub-control volumes for a...
Definition fvbaselocalresidual.hh:205
Opm::FvBaseLocalResidual::eval
void eval(const Problem &problem, const Element &element)
Compute the local residual, i.e.
Definition fvbaselocalresidual.hh:128
Opm::FvBaseLocalResidual::evalBoundarySegment_
void evalBoundarySegment_(LocalEvalBlockVector &residual, const BoundaryContext &boundaryCtx, unsigned boundaryFaceIdx, unsigned timeIdx) const
Evaluate all boundary conditions for a single sub-control volume face to the local residual.
Definition fvbaselocalresidual.hh:463
Opm::FvBaseLocalResidual::eval
void eval(LocalEvalBlockVector &residual, ElementContext &elemCtx) const
Compute the local residual, i.e.
Definition fvbaselocalresidual.hh:158
Opm::FvBaseLocalResidual::computeFlux
void computeFlux(RateVector &, const ElementContext &, unsigned, unsigned) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume.
Definition fvbaselocalresidual.hh:394
Opm::FvBaseLocalResidual::evalBoundary_
void evalBoundary_(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Evaluate the boundary conditions of an element.
Definition fvbaselocalresidual.hh:422
Opm::FvBaseLocalResidual::residual
const LocalEvalBlockVector & residual() const
Return the result of the eval() call using internal storage.
Definition fvbaselocalresidual.hh:106
Opm::FvBaseLocalResidual::eval
void eval(ElementContext &elemCtx)
Compute the local residual, i.e.
Definition fvbaselocalresidual.hh:144
Opm::FvBaseLocalResidual::residual
const EvalVector & residual(unsigned dofIdx) const
Return the result of the eval() call using internal storage.
Definition fvbaselocalresidual.hh:115
Opm::aligned_allocator
Definition alignedallocator.hh:97
fvbaseproperties.hh
Declare the properties used by the infrastructure code of the finite volume discretizations.
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
parametersystem.hpp
This file provides the infrastructure to retrieve run-time parameters.