Loading...
Searching...
No Matches
fvbaseadlocallinearizer.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_AD_LOCAL_LINEARIZER_HH
#define EWOMS_FV_BASE_AD_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/Valgrind.hpp>
#include <opm/material/densead/Math.hpp>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/simulators/linalg/linalgproperties.hh>
#include <cstddef>
#include <memory>
namespace Opm {
// forward declaration
template<class TypeTag>
class FvBaseAdLocalLinearizer;
}
namespace Opm::Properties {
// declare the property tags required for the finite differences local linearizer
namespace TTag {
struct AutoDiffLocalLinearizer {};
} // namespace TTag
// set the properties to be spliced in
template<class TypeTag>
struct LocalLinearizer<TypeTag, TTag::AutoDiffLocalLinearizer>
{ using type = FvBaseAdLocalLinearizer<TypeTag>; };
template<class TypeTag>
struct Evaluation<TypeTag, TTag::AutoDiffLocalLinearizer>
{
private:
static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = DenseAd::Evaluation<Scalar, numEq>;
};
} // namespace Opm::Properties
namespace Opm {
template<class TypeTag>
class FvBaseAdLocalLinearizer
{
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>() };
using ScalarVectorBlock = Dune::FieldVector<Scalar, numEq>;
// extract local matrices from jacobian matrix for consistency
using ScalarMatrixBlock = typename GetPropType<TypeTag, Properties::SparseMatrixAdapter>::MatrixBlock;
using ScalarLocalBlockVector = Dune::BlockVector<ScalarVectorBlock>;
using ScalarLocalBlockMatrix = Dune::Matrix<ScalarMatrixBlock>;
public:
FvBaseAdLocalLinearizer() = default;
// copying local linearizer objects around is a very bad idea, so we explicitly
// prevent it...
FvBaseAdLocalLinearizer(const FvBaseAdLocalLinearizer&) = delete;
static void registerParameters()
{}
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.updateStencil(elem);
elemCtx.updateAllIntensiveQuantities();
// update the weights of the primary variables for the context
model_().updatePVWeights(elemCtx);
// resize the internal arrays of the linearizer
resize_(elemCtx);
// compute the local residual and its Jacobian
const unsigned numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
for (unsigned focusDofIdx = 0; focusDofIdx < numPrimaryDof; ++focusDofIdx) {
elemCtx.setFocusDofIndex(focusDofIdx);
elemCtx.updateAllExtensiveQuantities();
// calculate the local residual
localResidual_.eval(elemCtx);
// convert the local Jacobian matrix and the right hand side from the data
// structures used by the automatic differentiation code to the conventional
// ones used by the linear solver.
updateLocalLinearization_(elemCtx, focusDofIdx);
}
}
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(); }
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);
if (jacobian_.N() != numDof || jacobian_.M() != numPrimaryDof) {
jacobian_.setSize(numDof, numPrimaryDof);
}
}
void reset_(const ElementContext&)
{
residual_ = 0.0;
jacobian_ = 0.0;
}
void updateLocalLinearization_(const ElementContext& elemCtx,
unsigned focusDofIdx)
{
const auto& resid = localResidual_.residual();
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
residual_[focusDofIdx][eqIdx] = resid[focusDofIdx][eqIdx].value();
}
const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
// A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of
// the residual function 'eqIdx' for the degree of freedom 'dofIdx'
// with regard to the focus variable 'pvIdx' of the degree of freedom
// 'focusDofIdx'
jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = resid[dofIdx][eqIdx].derivative(pvIdx);
Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
}
}
}
}
Simulator* simulatorPtr_{};
std::unique_ptr<ElementContext> internalElemContext_{};
LocalResidual localResidual_{};
ScalarLocalBlockVector residual_{};
ScalarLocalBlockMatrix jacobian_{};
};
} // namespace Opm
#endif
Opm::FvBaseAdLocalLinearizer
Calculates the local residual and its Jacobian for a single element of the grid.
Definition fvbaseadlocallinearizer.hh:93
Opm::FvBaseAdLocalLinearizer::residual
const ScalarVectorBlock & residual(unsigned dofIdx) const
Returns the local residual of a sub-control volume.
Definition fvbaseadlocallinearizer.hh:228
Opm::FvBaseAdLocalLinearizer::linearize
void linearize(const Element &element)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition fvbaseadlocallinearizer.hh:153
Opm::FvBaseAdLocalLinearizer::localResidual
LocalResidual & localResidual()
Return reference to the local residual.
Definition fvbaseadlocallinearizer.hh:203
Opm::FvBaseAdLocalLinearizer::updateLocalLinearization_
void updateLocalLinearization_(const ElementContext &elemCtx, unsigned focusDofIdx)
Updates the current local Jacobian matrix with the partial derivatives of all equations for the degre...
Definition fvbaseadlocallinearizer.hh:275
Opm::FvBaseAdLocalLinearizer::reset_
void reset_(const ElementContext &)
Reset the all relevant internal attributes to 0.
Definition fvbaseadlocallinearizer.hh:265
Opm::FvBaseAdLocalLinearizer::resize_
void resize_(const ElementContext &elemCtx)
Resize all internal attributes to the size of the element.
Definition fvbaseadlocallinearizer.hh:251
Opm::FvBaseAdLocalLinearizer::localResidual
const LocalResidual & localResidual() const
Return reference to the local residual.
Definition fvbaseadlocallinearizer.hh:209
Opm::FvBaseAdLocalLinearizer::init
void init(Simulator &simulator)
Initialize the local Jacobian object.
Definition fvbaseadlocallinearizer.hh:136
Opm::FvBaseAdLocalLinearizer::registerParameters
static void registerParameters()
Register all run-time parameters for the local jacobian.
Definition fvbaseadlocallinearizer.hh:125
Opm::FvBaseAdLocalLinearizer::linearize
void linearize(ElementContext &elemCtx, const Element &elem)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition fvbaseadlocallinearizer.hh:173
Opm::FvBaseAdLocalLinearizer::jacobian
const ScalarMatrixBlock & jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
Returns the local Jacobian matrix of the residual of a sub-control volume.
Definition fvbaseadlocallinearizer.hh:220
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
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::AutoDiffLocalLinearizer
Definition fvbaseadlocallinearizer.hh:58