28#ifndef EWOMS_FV_BASE_LOCAL_RESIDUAL_HH
29#define EWOMS_FV_BASE_LOCAL_RESIDUAL_HH
36#include <opm/material/common/Valgrind.hpp>
38#include <dune/istl/bvector.hh>
39#include <dune/grid/common/geometry.hh>
41#include <dune/common/fvector.hh>
43#include <dune/common/classname.hh>
56template<
class TypeTag>
63 using Element =
typename GridView::template Codim<0>::Entity;
81 using EvalVector = Dune::FieldVector<Evaluation, numEq>;
88 using LocalEvalBlockVector = Dune::BlockVector<EvalVector,
aligned_allocator<EvalVector,
alignof(EvalVector)> >;
107 {
return internalResidual_; }
116 {
return internalResidual_[dofIdx]; }
128 void eval(
const Problem& problem,
const Element& element)
130 ElementContext elemCtx(problem);
131 elemCtx.updateAll(element);
144 void eval(ElementContext& elemCtx)
146 size_t numDof = elemCtx.numDof(0);
147 internalResidual_.resize(numDof);
148 asImp_().eval(internalResidual_, elemCtx);
159 ElementContext& elemCtx)
const
166 asImp_().evalFluxes(
residual, elemCtx, 0);
169 asImp_().evalVolumeTerms_(
residual, elemCtx);
172 asImp_().evalBoundary_(
residual, elemCtx, 0);
174 if (useVolumetricResidual) {
177 size_t numDof = elemCtx.numDof(0);
178 for (
unsigned dofIdx=0; dofIdx < numDof; ++dofIdx) {
179 if (elemCtx.dofTotalVolume(dofIdx, 0) > 0.0) {
181 Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, 0);
183 assert(std::isfinite(dofVolume));
184 Valgrind::CheckDefined(dofVolume);
205 const ElementContext& elemCtx,
214 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
215 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
220 elemCtx.stencil(
timeIdx).subControlVolume(dofIdx).volume()
221 * elemCtx.intensiveQuantities(dofIdx,
timeIdx).extrusionFactor();
227 if (dofIdx == elemCtx.focusDofIndex()) {
228 asImp_().computeStorage(
storage[dofIdx],
237 Dune::FieldVector<Scalar, numEq> tmp;
238 asImp_().computeStorage(tmp,
251 if (elemCtx.enableStorageCache()) {
252 size_t numPrimaryDof = elemCtx.numPrimaryDof(
timeIdx);
253 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
254 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx,
timeIdx);
255 const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx,
timeIdx);
263 Dune::FieldVector<Scalar, numEq> tmp;
264 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
265 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
267 asImp_().computeStorage(tmp,
272 elemCtx.stencil(
timeIdx).subControlVolume(dofIdx).volume()
273 * elemCtx.intensiveQuantities(dofIdx,
timeIdx).extrusionFactor();
282 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
283 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
300 const ElementContext& elemCtx,
305 const auto& stencil = elemCtx.stencil(
timeIdx);
307 size_t numInteriorFaces = elemCtx.numInteriorFaces(
timeIdx);
310 unsigned i =
face.interiorIndex();
311 unsigned j =
face.exteriorIndex();
313 Valgrind::SetUndefined(
flux);
315 Valgrind::CheckDefined(
flux);
321 Scalar alpha = elemCtx.extensiveQuantities(
scvfIdx,
timeIdx).extrusionFactor();
322 alpha *=
face.area();
323 Valgrind::CheckDefined(alpha);
352 size_t numDof = elemCtx.numDof(
timeIdx);
353 for (
unsigned i=0; i < numDof; i++) {
354 for (
unsigned j = 0; j < numEq; ++ j) {
356 Valgrind::CheckDefined(
residual[i][j]);
376 const ElementContext&,
380 throw std::logic_error(
"Not implemented: The local residual "+Dune::className<Implementation>()
381 +
" does not implement the required method 'computeStorage()'");
392 const ElementContext&,
396 throw std::logic_error(
"Not implemented: The local residual "+Dune::className<Implementation>()
397 +
" does not implement the required method 'computeFlux()'");
407 const ElementContext&,
411 throw std::logic_error(
"Not implemented: The local residual "+Dune::className<Implementation>()
412 +
" does not implement the required method 'computeSource()'");
420 const ElementContext& elemCtx,
423 if (!elemCtx.onBoundary())
432 size_t numBoundaryFaces =
boundaryCtx.numBoundaryFaces(0);
433 for (
unsigned faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx,
boundaryCtx.increment()) {
444 size_t numDof = elemCtx.numDof(0);
445 for (
unsigned i=0; i < numDof; i++) {
446 for (
unsigned j = 0; j < numEq; ++ j) {
448 Valgrind::CheckDefined(
residual[i][j]);
464 BoundaryRateVector values;
466 Valgrind::SetUndefined(values);
468 Valgrind::CheckDefined(values);
471 unsigned dofIdx = stencil.boundaryFace(
boundaryFaceIdx).interiorIndex();
478 Valgrind::CheckDefined(values[
eqIdx]);
492 ElementContext& elemCtx)
const
502 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
503 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
504 Scalar extrusionFactor =
505 elemCtx.intensiveQuantities(dofIdx, 0).extrusionFactor();
506 Valgrind::CheckDefined(extrusionFactor);
508 assert(extrusionFactor > 0.0);
510 elemCtx.stencil(0).subControlVolume(dofIdx).volume() * extrusionFactor;
518 if (!extensiveStorageTerm &&
519 !std::is_same<Scalar, Evaluation>::value &&
520 dofIdx != elemCtx.focusDofIndex())
522 asImp_().computeStorage(
tmp2, elemCtx, dofIdx, 0);
527 asImp_().computeStorage(tmp, elemCtx, dofIdx, 0);
530 Valgrind::CheckDefined(tmp);
535 if (elemCtx.enableStorageCache()) {
536 const auto& model = elemCtx.model();
537 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
538 if (model.newtonMethod().numIterations() == 0 &&
539 !elemCtx.haveStashedIntensiveQuantities())
541 if (!elemCtx.problem().recycleFirstIterationStorage()) {
546 elemCtx.updatePrimaryIntensiveQuantities(1);
547 asImp_().computeStorage(
tmp2, elemCtx, dofIdx, 1);
561 Valgrind::CheckDefined(
tmp2);
563 model.updateCachedStorage(globalDofIdx, 1,
tmp2);
569 tmp2 = model.cachedStorage(globalDofIdx, 1);
570 Valgrind::CheckDefined(
tmp2);
577 asImp_().computeStorage(
tmp2, elemCtx, dofIdx, 1);
578 Valgrind::CheckDefined(
tmp2);
583 double dt = elemCtx.simulator().timeStepSize();
591 Valgrind::CheckDefined(
residual[dofIdx]);
594 asImp_().computeSource(
sourceRate, elemCtx, dofIdx, 0);
599 if (!extensiveStorageTerm &&
600 !std::is_same<Scalar, Evaluation>::value &&
601 dofIdx != elemCtx.focusDofIndex())
613 Valgrind::CheckDefined(
residual[dofIdx]);
618 size_t numDof = elemCtx.numDof(0);
619 for (
unsigned i=0; i < numDof; i++) {
620 for (
unsigned j = 0; j < numEq; ++ j) {
622 Valgrind::CheckDefined(
residual[i][j]);
630 Implementation& asImp_()
631 {
return *
static_cast<Implementation*
>(
this); }
633 const Implementation& asImp_()
const
634 {
return *
static_cast<const Implementation*
>(
this); }
636 LocalEvalBlockVector internalResidual_;
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1....
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition fvbaselocalresidual.hh:58
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:491
static void registerParameters()
Register all run-time parameters for the local residual.
Definition fvbaselocalresidual.hh:99
void computeSource(RateVector &, const ElementContext &, unsigned, unsigned) const
Calculate the source term of the equation.
Definition fvbaselocalresidual.hh:406
void evalFluxes(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Add the flux term to a local residual.
Definition fvbaselocalresidual.hh:299
void computeStorage(EqVector &, const ElementContext &, unsigned, unsigned) const
Evaluate the amount all conservation quantities (e.g.
Definition fvbaselocalresidual.hh:375
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:204
void eval(const Problem &problem, const Element &element)
Compute the local residual, i.e.
Definition fvbaselocalresidual.hh:128
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:459
void eval(LocalEvalBlockVector &residual, ElementContext &elemCtx) const
Compute the local residual, i.e.
Definition fvbaselocalresidual.hh:158
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:391
void evalBoundary_(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Evaluate the boundary conditions of an element.
Definition fvbaselocalresidual.hh:419
const LocalEvalBlockVector & residual() const
Return the result of the eval() call using internal storage.
Definition fvbaselocalresidual.hh:106
void eval(ElementContext &elemCtx)
Compute the local residual, i.e.
Definition fvbaselocalresidual.hh:144
const EvalVector & residual(unsigned dofIdx) const
Return the result of the eval() call using internal storage.
Definition fvbaselocalresidual.hh:115
Definition alignedallocator.hh:94
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
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
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
This file provides the infrastructure to retrieve run-time parameters.