28 #ifndef EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
29 #define EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
34 #include <opm/material/common/MathToolbox.hpp>
35 #include <opm/common/Valgrind.hpp>
37 #include <dune/istl/bvector.hh>
38 #include <dune/istl/matrix.hh>
40 #include <dune/common/fvector.hh>
41 #include <dune/common/fmatrix.hh>
47 template<
class TypeTag>
50 namespace Properties {
71 SET_TYPE_PROP(FiniteDifferenceLocalLinearizer, LocalLinearizer,
84 SET_INT_PROP(FiniteDifferenceLocalLinearizer, NumericDifferenceMethod, +1);
89 std::max<Scalar>(0.9123e-10, std::numeric_limits<Scalar>::epsilon()*1.23e3));
128 template<
class TypeTag>
132 typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizer) Implementation;
133 typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
135 typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
136 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
137 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
138 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
139 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
140 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
141 typedef typename GridView::template Codim<0>::Entity Element;
144 typedef Dune::FieldMatrix<Scalar, numEq, numEq> ScalarMatrixBlock;
145 typedef Dune::FieldVector<Scalar, numEq> ScalarVectorBlock;
147 typedef Dune::BlockVector<ScalarVectorBlock> ScalarLocalBlockVector;
148 typedef Dune::Matrix<ScalarMatrixBlock> ScalarLocalBlockMatrix;
150 typedef typename LocalResidual::LocalEvalBlockVector LocalEvalBlockVector;
152 #if __GNUC__ == 4 && __GNUC_MINOR__ <= 6
156 FvBaseFdLocalLinearizer(
const FvBaseFdLocalLinearizer&)
157 : internalElemContext_(0)
163 FvBaseFdLocalLinearizer(
const FvBaseFdLocalLinearizer&) =
delete;
166 FvBaseFdLocalLinearizer()
167 : internalElemContext_(0)
170 ~FvBaseFdLocalLinearizer()
171 {
delete internalElemContext_; }
179 "The method used for numeric differentiation (-1: backward "
180 "differences, 0: central differences, 1: forward differences)");
193 simulatorPtr_ = &simulator;
194 delete internalElemContext_;
195 internalElemContext_ =
new ElementContext(simulator);
211 linearize(*internalElemContext_, element);
228 void linearize(ElementContext& elemCtx,
const Element& elem)
230 elemCtx.updateAll(elem);
233 model_().updatePVWeights(elemCtx);
239 localResidual_.eval(residual_, elemCtx);
242 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
243 for (
unsigned dofIdx = 0; dofIdx < numPrimaryDof; dofIdx++) {
244 for (
unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
245 asImp_().evalPartialDerivative_(elemCtx, dofIdx, pvIdx);
273 unsigned pvIdx)
const
275 unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
276 Scalar pvWeight = elemCtx.model().primaryVarWeight(globalIdx, pvIdx);
277 assert(pvWeight > 0 && std::isfinite(pvWeight));
278 Opm::Valgrind::CheckDefined(pvWeight);
287 {
return localResidual_; }
293 {
return localResidual_; }
303 const ScalarMatrixBlock&
jacobian(
unsigned domainScvIdx,
unsigned rangeScvIdx)
const
304 {
return jacobian_[domainScvIdx][rangeScvIdx]; }
311 const ScalarVectorBlock&
residual(
unsigned dofIdx)
const
312 {
return residual_[dofIdx]; }
315 Implementation& asImp_()
316 {
return *
static_cast<Implementation*
>(
this); }
317 const Implementation& asImp_()
const
318 {
return *
static_cast<const Implementation*
>(
this); }
320 const Simulator& simulator_()
const
321 {
return *simulatorPtr_; }
322 const Problem& problem_()
const
323 {
return simulatorPtr_->
problem(); }
324 const Model& model_()
const
325 {
return simulatorPtr_->
model(); }
339 size_t numDof = elemCtx.numDof(0);
340 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
342 residual_.resize(numDof);
343 jacobian_.setSize(numDof, numPrimaryDof);
345 derivResidual_.resize(numDof);
351 void reset_(
const ElementContext& elemCtx)
353 size_t numDof = elemCtx.numDof(0);
354 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
355 for (
unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx)
356 for (
unsigned dof2Idx = 0; dof2Idx < numDof; ++ dof2Idx)
357 jacobian_[dof2Idx][primaryDofIdx] = 0.0;
359 for (
unsigned primaryDofIdx = 0; primaryDofIdx < numDof; ++ primaryDofIdx)
360 residual_[primaryDofIdx] = 0.0;
411 elemCtx.stashIntensiveQuantities(dofIdx);
413 PrimaryVariables priVars(elemCtx.primaryVars(dofIdx, 0));
414 Scalar eps = asImp_().numericEpsilon(elemCtx, dofIdx, pvIdx);
422 priVars[pvIdx] += eps;
426 elemCtx.updateIntensiveQuantities(priVars, dofIdx, 0);
427 elemCtx.updateAllExtensiveQuantities();
428 localResidual_.eval(derivResidual_, elemCtx);
434 derivResidual_ = residual_;
442 priVars[pvIdx] -= delta + eps;
447 elemCtx.updateIntensiveQuantities(priVars, dofIdx, 0);
448 elemCtx.updateAllExtensiveQuantities();
449 localResidual_.eval(elemCtx);
451 derivResidual_ -= localResidual_.residual();
457 derivResidual_ -= residual_;
464 derivResidual_ /= delta;
468 elemCtx.restoreIntensiveQuantities(dofIdx);
471 for (
unsigned i = 0; i < derivResidual_.size(); ++i)
472 Opm::Valgrind::CheckDefined(derivResidual_[i]);
482 unsigned focusDofIdx,
485 size_t numDof = elemCtx.numDof(0);
486 for (
unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
487 for (
unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
492 jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
493 Opm::Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
501 ElementContext *internalElemContext_;
503 LocalEvalBlockVector residual_;
504 LocalEvalBlockVector derivResidual_;
505 ScalarLocalBlockMatrix jacobian_;
507 LocalResidual localResidual_;
static void registerParameters()
Register all run-time parameters for the local jacobian.
Definition: fvbasefdlocallinearizer.hh:176
const ScalarVectorBlock & residual(unsigned dofIdx) const
Returns the local residual of a sub-control volume.
Definition: fvbasefdlocallinearizer.hh:311
void linearize(ElementContext &elemCtx, const Element &elem)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition: fvbasefdlocallinearizer.hh:228
void evalPartialDerivative_(ElementContext &elemCtx, unsigned dofIdx, unsigned pvIdx)
Compute the partial derivatives of a context's residual functions.
Definition: fvbasefdlocallinearizer.hh:405
const LocalResidual & localResidual() const
Return reference to the local residual.
Definition: fvbasefdlocallinearizer.hh:292
Model & model()
Return the physical model used in the simulation.
Definition: simulator.hh:189
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
static Scalar baseEpsilon()
Returns the unweighted epsilon value used to calculate the local derivatives.
Definition: fvbasefdlocallinearizer.hh:257
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:271
const ScalarMatrixBlock & jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
Returns the local Jacobian matrix of the residual of a sub-control volume.
Definition: fvbasefdlocallinearizer.hh:303
Problem & problem()
Return the object which specifies the pysical setup of the simulation.
Definition: simulator.hh:202
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Calculates the Jacobian of the local residual for finite volume spatial discretizations using a finit...
Definition: fvbasefdlocallinearizer.hh:48
#define SET_INT_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant integer value.
Definition: propertysystem.hh:345
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:481
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
This file provides the infrastructure to retrieve run-time parameters.
static int numericDifferenceMethod_()
Returns the numeric difference method which is applied.
Definition: fvbasefdlocallinearizer.hh:330
void resize_(const ElementContext &elemCtx)
Resize all internal attributes to the size of the element.
Definition: fvbasefdlocallinearizer.hh:337
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
void linearize(const Element &element)
Compute an element's local Jacobian matrix and evaluate its residual.
Definition: fvbasefdlocallinearizer.hh:209
LocalResidual & localResidual()
Return reference to the local residual.
Definition: fvbasefdlocallinearizer.hh:286
void init(Simulator &simulator)
Initialize the local Jacobian object.
Definition: fvbasefdlocallinearizer.hh:191
void reset_(const ElementContext &elemCtx)
Reset the all relevant internal attributes to 0.
Definition: fvbasefdlocallinearizer.hh:351
Provides the magic behind the eWoms property system.
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394