Loading...
Searching...
No Matches
fvbasefdlocallinearizer.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_FD_LOCAL_LINEARIZER_HH
#define EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/matrix.hh>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/linalg/linalgproperties.hh>
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <memory>
namespace Opm {
// forward declaration
template<class TypeTag>
class FvBaseFdLocalLinearizer;
} // namespace Opm
namespace Opm::Properties {
// declare the property tags required for the finite differences local linearizer
namespace TTag {
struct FiniteDifferenceLocalLinearizer {};
} // namespace TTag
template<class TypeTag, class MyTypeTag>
struct BaseEpsilon { using type = UndefinedProperty; };
// set the properties to be spliced in
template<class TypeTag>
struct LocalLinearizer<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
{ using type = FvBaseFdLocalLinearizer<TypeTag>; };
template<class TypeTag>
struct Evaluation<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
{ using type = GetPropType<TypeTag, Properties::Scalar>; };
template<class TypeTag>
struct BaseEpsilon<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
{
using type = GetPropType<TypeTag, Properties::Scalar>;
static constexpr type value = std::max<type>(0.9123e-10,
std::numeric_limits<type>::epsilon() * 1.23e3);
};
} // namespace Opm::Properties
namespace Opm::Parameters {
struct NumericDifferenceMethod { static constexpr int value = +1; };
} // namespace Opm::Parameters
namespace Opm {
template<class TypeTag>
class FvBaseFdLocalLinearizer
{
private:
using Implementation = GetPropType<TypeTag, Properties::LocalLinearizer>;
using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using Model = GetPropType<TypeTag, Properties::Model>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Element = typename GridView::template Codim<0>::Entity;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
// extract local matrices from jacobian matrix for consistency
using ScalarMatrixBlock = typename GetPropType<TypeTag, Properties::SparseMatrixAdapter>::MatrixBlock;
using ScalarVectorBlock = Dune::FieldVector<Scalar, numEq>;
using ScalarLocalBlockVector = Dune::BlockVector<ScalarVectorBlock>;
using ScalarLocalBlockMatrix = Dune::Matrix<ScalarMatrixBlock>;
using LocalEvalBlockVector = typename LocalResidual::LocalEvalBlockVector;
public:
FvBaseFdLocalLinearizer() = default;
// Disable copying for efficiency reasons
FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&) = delete;
static void registerParameters()
{
Parameters::Register<Parameters::NumericDifferenceMethod>
("The method used for numeric differentiation (-1: backward "
"differences, 0: central differences, 1: forward differences)");
}
void init(Simulator& simulator)
{
simulatorPtr_ = &simulator;
internalElemContext_ = std::make_unique<ElementContext>(simulator);
}
void linearize(const Element& element)
{
linearize(*internalElemContext_, element);
}
void linearize(ElementContext& elemCtx, const Element& elem)
{
elemCtx.updateAll(elem);
// update the weights of the primary variables for the context
model_().updatePVWeights(elemCtx);
resize_(elemCtx);
reset_(elemCtx);
// calculate the local residual
localResidual_.eval(residual_, elemCtx);
// calculate the local jacobian matrix
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
for (auto dofIdx = 0 * numPrimaryDof; dofIdx < numPrimaryDof; ++dofIdx) {
for (auto pvIdx = 0 * numEq; pvIdx < numEq; ++pvIdx) {
asImp_().evalPartialDerivative_(elemCtx, dofIdx, pvIdx);
// incorporate the partial derivatives into the local Jacobian matrix
updateLocalJacobian_(elemCtx, dofIdx, pvIdx);
}
}
}
static Scalar baseEpsilon()
{ return getPropValue<TypeTag, Properties::BaseEpsilon>(); }
Scalar numericEpsilon(const ElementContext& elemCtx,
unsigned dofIdx,
unsigned pvIdx) const
{
const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
const Scalar pvWeight = elemCtx.model().primaryVarWeight(globalIdx, pvIdx);
assert(pvWeight > 0 && std::isfinite(pvWeight));
Valgrind::CheckDefined(pvWeight);
return baseEpsilon() / pvWeight;
}
LocalResidual& localResidual()
{ return localResidual_; }
const LocalResidual& localResidual() const
{ return localResidual_; }
const ScalarMatrixBlock& jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
{ return jacobian_[domainScvIdx][rangeScvIdx]; }
const ScalarVectorBlock& residual(unsigned dofIdx) const
{ return residual_[dofIdx]; }
protected:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation& asImp_() const
{ return *static_cast<const Implementation*>(this); }
const Simulator& simulator_() const
{ return *simulatorPtr_; }
const Problem& problem_() const
{ return simulatorPtr_->problem(); }
const Model& model_() const
{ return simulatorPtr_->model(); }
static int numericDifferenceMethod_()
{
static int diff = Parameters::Get<Parameters::NumericDifferenceMethod>();
return diff;
}
void resize_(const ElementContext& elemCtx)
{
const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
residual_.resize(numDof);
jacobian_.setSize(numDof, numPrimaryDof);
derivResidual_.resize(numDof);
}
void reset_(const ElementContext& elemCtx)
{
const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
for (unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++primaryDofIdx) {
for (unsigned dof2Idx = 0; dof2Idx < numDof; ++dof2Idx) {
jacobian_[dof2Idx][primaryDofIdx] = 0.0;
}
}
for (unsigned primaryDofIdx = 0; primaryDofIdx < numDof; ++primaryDofIdx) {
residual_[primaryDofIdx] = 0.0;
}
}
void evalPartialDerivative_(ElementContext& elemCtx,
unsigned dofIdx,
unsigned pvIdx)
{
// save all quantities which depend on the specified primary
// variable at the given sub control volume
elemCtx.stashIntensiveQuantities(dofIdx);
PrimaryVariables priVars(elemCtx.primaryVars(dofIdx, /*timeIdx=*/0));
const Scalar eps = asImp_().numericEpsilon(elemCtx, dofIdx, pvIdx);
Scalar delta = 0.0;
if (numericDifferenceMethod_() >= 0) {
// we are not using backward differences, i.e. we need to
// calculate f(x + \epsilon)
// deflect primary variables
priVars[pvIdx] += eps;
delta += eps;
// calculate the deflected residual
elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
elemCtx.updateAllExtensiveQuantities();
localResidual_.eval(derivResidual_, elemCtx);
}
else {
// we are using backward differences, i.e. we don't need
// to calculate f(x + \epsilon) and we can recycle the
// (already calculated) residual f(x)
derivResidual_ = residual_;
}
if (numericDifferenceMethod_() <= 0) {
// we are not using forward differences, i.e. we don't
// need to calculate f(x - \epsilon)
// deflect the primary variables
priVars[pvIdx] -= delta + eps;
delta += eps;
// calculate the deflected residual again, this time we use the local
// residual's internal storage.
elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
elemCtx.updateAllExtensiveQuantities();
localResidual_.eval(elemCtx);
derivResidual_ -= localResidual_.residual();
}
else {
// we are using forward differences, i.e. we don't need to
// calculate f(x - \epsilon) and we can recycle the
// (already calculated) residual f(x)
derivResidual_ -= residual_;
}
assert(delta > 0);
// divide difference in residuals by the magnitude of the
// deflections between the two function evaluation
derivResidual_ /= delta;
// restore the original state of the element's volume
// variables
elemCtx.restoreIntensiveQuantities(dofIdx);
#ifndef NDEBUG
for (unsigned i = 0; i < derivResidual_.size(); ++i) {
Valgrind::CheckDefined(derivResidual_[i]);
}
#endif
}
void updateLocalJacobian_(const ElementContext& elemCtx,
unsigned focusDofIdx,
unsigned pvIdx)
{
const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
// A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of the
// residual function 'eqIdx' for the degree of freedom 'dofIdx' with
// regard to the primary variable 'pvIdx' of the degree of freedom
// 'focusDofIdx'
jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
}
}
}
Simulator* simulatorPtr_{};
std::unique_ptr<ElementContext> internalElemContext_{};
LocalEvalBlockVector residual_{};
LocalEvalBlockVector derivResidual_{};
ScalarLocalBlockMatrix jacobian_{};
LocalResidual localResidual_{};
};
} // namespace Opm
#endif
Opm::FvBaseFdLocalLinearizer
Calculates the Jacobian of the local residual for finite volume spatial discretizations using a finit...
Definition fvbasefdlocallinearizer.hh:144
Opm::FvBaseFdLocalLinearizer::evalPartialDerivative_
void evalPartialDerivative_(ElementContext &elemCtx, unsigned dofIdx, unsigned pvIdx)
Compute the partial derivatives of a context's residual functions.
Definition fvbasefdlocallinearizer.hh:415
Opm::FvBaseFdLocalLinearizer::linearize
void linearize(const Element &element)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition fvbasefdlocallinearizer.hh:209
Opm::FvBaseFdLocalLinearizer::numericEpsilon
Scalar numericEpsilon(const ElementContext &elemCtx, unsigned dofIdx, unsigned pvIdx) const
Returns the epsilon value which is added and removed from the current solution.
Definition fvbasefdlocallinearizer.hh:272
Opm::FvBaseFdLocalLinearizer::baseEpsilon
static Scalar baseEpsilon()
Returns the unweighted epsilon value used to calculate the local derivatives.
Definition fvbasefdlocallinearizer.hh:258
Opm::FvBaseFdLocalLinearizer::updateLocalJacobian_
void updateLocalJacobian_(const ElementContext &elemCtx, unsigned focusDofIdx, unsigned pvIdx)
Updates the current local Jacobian matrix with the partial derivatives of all equations for primary v...
Definition fvbasefdlocallinearizer.hh:492
Opm::FvBaseFdLocalLinearizer::numericDifferenceMethod_
static int numericDifferenceMethod_()
Returns the numeric difference method which is applied.
Definition fvbasefdlocallinearizer.hh:334
Opm::FvBaseFdLocalLinearizer::jacobian
const ScalarMatrixBlock & jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
Returns the local Jacobian matrix of the residual of a sub-control volume.
Definition fvbasefdlocallinearizer.hh:304
Opm::FvBaseFdLocalLinearizer::registerParameters
static void registerParameters()
Register all run-time parameters for the local jacobian.
Definition fvbasefdlocallinearizer.hh:177
Opm::FvBaseFdLocalLinearizer::localResidual
const LocalResidual & localResidual() const
Return reference to the local residual.
Definition fvbasefdlocallinearizer.hh:293
Opm::FvBaseFdLocalLinearizer::init
void init(Simulator &simulator)
Initialize the local Jacobian object.
Definition fvbasefdlocallinearizer.hh:192
Opm::FvBaseFdLocalLinearizer::residual
const ScalarVectorBlock & residual(unsigned dofIdx) const
Returns the local residual of a sub-control volume.
Definition fvbasefdlocallinearizer.hh:312
Opm::FvBaseFdLocalLinearizer::linearize
void linearize(ElementContext &elemCtx, const Element &elem)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition fvbasefdlocallinearizer.hh:229
Opm::FvBaseFdLocalLinearizer::reset_
void reset_(const ElementContext &elemCtx)
Reset the all relevant internal attributes to 0.
Definition fvbasefdlocallinearizer.hh:358
Opm::FvBaseFdLocalLinearizer::localResidual
LocalResidual & localResidual()
Return reference to the local residual.
Definition fvbasefdlocallinearizer.hh:287
Opm::FvBaseFdLocalLinearizer::resize_
void resize_(const ElementContext &elemCtx)
Resize all internal attributes to the size of the element.
Definition fvbasefdlocallinearizer.hh:344
Opm::Simulator
Manages the initializing and running of time dependent problems.
Definition simulator.hh:84
fvbaseproperties.hh
Declare the properties used by the infrastructure code of the finite volume discretizations.
linalgproperties.hh
Declares the properties required by the black oil model.
Opm::Properties::TTag
The generic type tag for problems using the immiscible multi-phase model.
Definition blackoilmodel.hh:81
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.
Opm::Parameters::Register
void Register(const char *usageString)
Register a run-time parameter.
Definition parametersystem.hpp:292
Opm::Parameters::Get
auto Get(bool errorIfNotRegistered=true)
Retrieve a runtime parameter.
Definition parametersystem.hpp:187
propertysystem.hh
The Opm property system, traits with inheritance.
Opm::Parameters::NumericDifferenceMethod
Specify which kind of method should be used to numerically calculate the partial derivatives of the r...
Definition fvbasefdlocallinearizer.hh:100
Opm::Properties::BaseEpsilon
Definition fvbasefdlocallinearizer.hh:69
Opm::Properties::Evaluation
Representation of a function evaluation and all necessary derivatives with regard to the intensive qu...
Definition fvbaseproperties.hh:66
Opm::Properties::LocalLinearizer
The type of the local linearizer.
Definition fvbaseproperties.hh:97
Opm::Properties::TTag::FiniteDifferenceLocalLinearizer
Definition fvbasefdlocallinearizer.hh:65
Opm::Properties::UndefinedProperty
a tag to mark properties as undefined
Definition propertysystem.hh:38