28#ifndef EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
29#define EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
35#include <opm/material/common/MathToolbox.hpp>
36#include <opm/material/common/Valgrind.hpp>
38#include <dune/istl/bvector.hh>
39#include <dune/istl/matrix.hh>
41#include <dune/common/fvector.hh>
42#include <dune/common/fmatrix.hh>
48template<
class TypeTag>
49class FvBaseFdLocalLinearizer;
53namespace Opm::Properties {
61template<
class TypeTag,
class MyTypeTag>
65template<
class TypeTag>
69template<
class TypeTag>
70struct Evaluation<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
74template<
class TypeTag>
75struct BaseEpsilon<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
78 static constexpr type value = std::max<type>(0.9123e-10, std::numeric_limits<type>::epsilon()*1.23e3);
83namespace Opm::Parameters {
134template<
class TypeTag>
147 using Element =
typename GridView::template Codim<0>::Entity;
153 using ScalarVectorBlock = Dune::FieldVector<Scalar, numEq>;
155 using ScalarLocalBlockVector = Dune::BlockVector<ScalarVectorBlock>;
156 using ScalarLocalBlockMatrix = Dune::Matrix<ScalarMatrixBlock>;
158 using LocalEvalBlockVector =
typename LocalResidual::LocalEvalBlockVector;
160#if __GNUC__ == 4 && __GNUC_MINOR__ <= 6
165 : internalElemContext_(0)
175 : internalElemContext_(0)
179 {
delete internalElemContext_; }
186 Parameters::Register<Parameters::NumericDifferenceMethod>
187 (
"The method used for numeric differentiation (-1: backward "
188 "differences, 0: central differences, 1: forward differences)");
199 void init(Simulator& simulator)
201 simulatorPtr_ = &simulator;
202 delete internalElemContext_;
203 internalElemContext_ =
new ElementContext(simulator);
219 linearize(*internalElemContext_, element);
236 void linearize(ElementContext& elemCtx,
const Element& elem)
238 elemCtx.updateAll(elem);
241 model_().updatePVWeights(elemCtx);
247 localResidual_.eval(residual_, elemCtx);
250 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
251 for (
unsigned dofIdx = 0; dofIdx < numPrimaryDof; dofIdx++) {
252 for (
unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
253 asImp_().evalPartialDerivative_(elemCtx, dofIdx, pvIdx);
281 unsigned pvIdx)
const
283 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
284 Scalar pvWeight = elemCtx.model().primaryVarWeight(globalIdx, pvIdx);
285 assert(pvWeight > 0 && std::isfinite(pvWeight));
286 Valgrind::CheckDefined(pvWeight);
295 {
return localResidual_; }
301 {
return localResidual_; }
311 const ScalarMatrixBlock&
jacobian(
unsigned domainScvIdx,
unsigned rangeScvIdx)
const
312 {
return jacobian_[domainScvIdx][rangeScvIdx]; }
319 const ScalarVectorBlock&
residual(
unsigned dofIdx)
const
320 {
return residual_[dofIdx]; }
323 Implementation& asImp_()
324 {
return *
static_cast<Implementation*
>(
this); }
325 const Implementation& asImp_()
const
326 {
return *
static_cast<const Implementation*
>(
this); }
328 const Simulator& simulator_()
const
329 {
return *simulatorPtr_; }
330 const Problem& problem_()
const
331 {
return simulatorPtr_->problem(); }
332 const Model& model_()
const
333 {
return simulatorPtr_->model(); }
340 static int diff = Parameters::Get<Parameters::NumericDifferenceMethod>();
350 size_t numDof = elemCtx.numDof(0);
351 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
353 residual_.resize(numDof);
354 jacobian_.setSize(numDof, numPrimaryDof);
356 derivResidual_.resize(numDof);
362 void reset_(
const ElementContext& elemCtx)
364 size_t numDof = elemCtx.numDof(0);
365 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
366 for (
unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx)
367 for (
unsigned dof2Idx = 0; dof2Idx < numDof; ++ dof2Idx)
368 jacobian_[dof2Idx][primaryDofIdx] = 0.0;
370 for (
unsigned primaryDofIdx = 0; primaryDofIdx < numDof; ++ primaryDofIdx)
371 residual_[primaryDofIdx] = 0.0;
422 elemCtx.stashIntensiveQuantities(dofIdx);
424 PrimaryVariables priVars(elemCtx.primaryVars(dofIdx, 0));
425 Scalar eps = asImp_().numericEpsilon(elemCtx, dofIdx, pvIdx);
433 priVars[pvIdx] += eps;
437 elemCtx.updateIntensiveQuantities(priVars, dofIdx, 0);
438 elemCtx.updateAllExtensiveQuantities();
439 localResidual_.eval(derivResidual_, elemCtx);
445 derivResidual_ = residual_;
453 priVars[pvIdx] -= delta + eps;
458 elemCtx.updateIntensiveQuantities(priVars, dofIdx, 0);
459 elemCtx.updateAllExtensiveQuantities();
460 localResidual_.eval(elemCtx);
462 derivResidual_ -= localResidual_.residual();
468 derivResidual_ -= residual_;
475 derivResidual_ /= delta;
479 elemCtx.restoreIntensiveQuantities(dofIdx);
482 for (
unsigned i = 0; i < derivResidual_.size(); ++i)
483 Valgrind::CheckDefined(derivResidual_[i]);
493 unsigned focusDofIdx,
496 size_t numDof = elemCtx.numDof(0);
497 for (
unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
498 for (
unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
503 jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
504 Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
512 ElementContext *internalElemContext_;
514 LocalEvalBlockVector residual_;
515 LocalEvalBlockVector derivResidual_;
516 ScalarLocalBlockMatrix jacobian_;
518 LocalResidual localResidual_;
Calculates the Jacobian of the local residual for finite volume spatial discretizations using a finit...
Definition fvbasefdlocallinearizer.hh:136
void evalPartialDerivative_(ElementContext &elemCtx, unsigned dofIdx, unsigned pvIdx)
Compute the partial derivatives of a context's residual functions.
Definition fvbasefdlocallinearizer.hh:416
void linearize(const Element &element)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition fvbasefdlocallinearizer.hh:217
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:279
static Scalar baseEpsilon()
Returns the unweighted epsilon value used to calculate the local derivatives.
Definition fvbasefdlocallinearizer.hh:265
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
static int numericDifferenceMethod_()
Returns the numeric difference method which is applied.
Definition fvbasefdlocallinearizer.hh:338
const ScalarMatrixBlock & jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
Returns the local Jacobian matrix of the residual of a sub-control volume.
Definition fvbasefdlocallinearizer.hh:311
static void registerParameters()
Register all run-time parameters for the local jacobian.
Definition fvbasefdlocallinearizer.hh:184
const LocalResidual & localResidual() const
Return reference to the local residual.
Definition fvbasefdlocallinearizer.hh:300
void init(Simulator &simulator)
Initialize the local Jacobian object.
Definition fvbasefdlocallinearizer.hh:199
const ScalarVectorBlock & residual(unsigned dofIdx) const
Returns the local residual of a sub-control volume.
Definition fvbasefdlocallinearizer.hh:319
void linearize(ElementContext &elemCtx, const Element &elem)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition fvbasefdlocallinearizer.hh:236
void reset_(const ElementContext &elemCtx)
Reset the all relevant internal attributes to 0.
Definition fvbasefdlocallinearizer.hh:362
LocalResidual & localResidual()
Return reference to the local residual.
Definition fvbasefdlocallinearizer.hh:294
void resize_(const ElementContext &elemCtx)
Resize all internal attributes to the size of the element.
Definition fvbasefdlocallinearizer.hh:348
Manages the initializing and running of time dependent problems.
Definition simulator.hh:92
Declare the properties used by the infrastructure code of the finite volume discretizations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
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:235
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
Specify which kind of method should be used to numerically calculate the partial derivatives of the r...
Definition fvbasefdlocallinearizer.hh:92
Definition fvbasefdlocallinearizer.hh:62
Representation of a function evaluation and all necessary derivatives with regard to the intensive qu...
Definition fvbaseproperties.hh:66
The type of the local linearizer.
Definition fvbaseproperties.hh:97
Definition fvbasefdlocallinearizer.hh:58
a tag to mark properties as undefined
Definition propertysystem.hh:40