28 #ifndef EWOMS_FV_BASE_LOCAL_RESIDUAL_HH 29 #define EWOMS_FV_BASE_LOCAL_RESIDUAL_HH 36 #include <opm/common/Valgrind.hpp> 37 #include <opm/common/Unused.hpp> 38 #include <opm/common/ErrorMacros.hpp> 39 #include <opm/common/Exceptions.hpp> 41 #include <dune/istl/bvector.hh> 42 #include <dune/grid/common/geometry.hh> 44 #include <dune/common/fvector.hh> 46 #include <dune/common/classname.hh> 59 template<
class TypeTag>
63 typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
65 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
66 typedef typename GridView::template Codim<0>::Entity Element;
68 typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
69 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
70 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
71 typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
72 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
73 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
74 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
75 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
76 typedef typename GET_PROP_TYPE(TypeTag, BoundaryContext) BoundaryContext;
78 static constexpr
bool useVolumetricResidual =
GET_PROP_VALUE(TypeTag, UseVolumetricResidual);
81 enum { extensiveStorageTerm =
GET_PROP_VALUE(TypeTag, ExtensiveStorageTerm) };
83 typedef Opm::MathToolbox<Evaluation> Toolbox;
84 typedef Dune::FieldVector<Evaluation, numEq> EvalVector;
91 typedef Dune::BlockVector<EvalVector, Ewoms::aligned_allocator<EvalVector, alignof(EvalVector)> > LocalEvalBlockVector;
110 {
return internalResidual_; }
119 {
return internalResidual_[dofIdx]; }
131 void eval(
const Problem& problem,
const Element& element)
133 ElementContext elemCtx(problem);
134 elemCtx.updateAll(element);
147 void eval(
const ElementContext& elemCtx)
149 size_t numDof = elemCtx.numDof(0);
150 internalResidual_.resize(numDof);
151 asImp_().eval(internalResidual_, elemCtx);
162 const ElementContext& elemCtx)
const 164 assert(
residual.size() == elemCtx.numDof(0));
169 asImp_().evalFluxes(
residual, elemCtx, 0);
172 asImp_().evalVolumeTerms_(
residual, elemCtx);
175 asImp_().evalBoundary_(
residual, elemCtx, 0);
177 if (useVolumetricResidual) {
180 size_t numDof = elemCtx.numDof(0);
181 for (
unsigned dofIdx=0; dofIdx < numDof; ++dofIdx) {
182 if (elemCtx.dofTotalVolume(dofIdx, 0) > 0.0) {
184 Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, 0);
186 assert(std::isfinite(dofVolume));
187 Opm::Valgrind::CheckDefined(dofVolume);
189 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
190 residual[dofIdx][eqIdx] /= dofVolume;
208 const ElementContext& elemCtx,
209 unsigned timeIdx)
const 217 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
218 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
219 storage[dofIdx] = 0.0;
223 elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
224 * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
230 if (dofIdx == elemCtx.focusDofIndex()) {
231 asImp_().computeStorage(storage[dofIdx],
236 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
237 storage[dofIdx][eqIdx] *= alpha;
240 Dune::FieldVector<Scalar, numEq> tmp;
241 asImp_().computeStorage(tmp,
246 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
247 storage[dofIdx][eqIdx] = tmp[eqIdx]*alpha;
254 if (elemCtx.enableStorageCache()) {
255 size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
256 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
257 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
258 const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx, timeIdx);
259 for (
unsigned eqIdx=0; eqIdx < numEq; eqIdx++)
260 storage[dofIdx][eqIdx] = cachedStorage[eqIdx];
266 Dune::FieldVector<Scalar, numEq> tmp;
267 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
268 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
270 asImp_().computeStorage(tmp,
275 elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
276 * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
278 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
279 storage[dofIdx][eqIdx] = tmp[eqIdx];
293 const ElementContext& elemCtx,
294 unsigned timeIdx)
const 298 const auto& stencil = elemCtx.stencil(timeIdx);
300 size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
301 for (
unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
302 const auto& face = stencil.interiorFace(scvfIdx);
303 unsigned i = face.interiorIndex();
304 unsigned j = face.exteriorIndex();
306 Opm::Valgrind::SetUndefined(flux);
307 asImp_().computeFlux(flux, elemCtx, scvfIdx, timeIdx);
308 Opm::Valgrind::CheckDefined(flux);
310 Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor();
311 alpha *= face.area();
312 Opm::Valgrind::CheckDefined(alpha);
313 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
314 flux[eqIdx] *= alpha;
329 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
337 size_t numDof = elemCtx.numDof(timeIdx);
338 for (
unsigned i=0; i < numDof; i++) {
339 for (
unsigned j = 0; j < numEq; ++ j) {
340 assert(std::isfinite(Toolbox::value(
residual[i][j])));
341 Opm::Valgrind::CheckDefined(
residual[i][j]);
361 const ElementContext& elemCtx OPM_UNUSED,
362 unsigned dofIdx OPM_UNUSED,
363 unsigned timeIdx OPM_UNUSED)
const 365 OPM_THROW(std::logic_error,
366 "Not implemented: The local residual " << Dune::className<Implementation>()
367 <<
" does not implement the required method 'computeStorage()'");
378 const ElementContext& elemCtx OPM_UNUSED,
379 unsigned scvfIdx OPM_UNUSED,
380 unsigned timeIdx OPM_UNUSED)
const 382 OPM_THROW(std::logic_error,
383 "Not implemented: The local residual " << Dune::className<Implementation>()
384 <<
" does not implement the required method 'computeFlux()'");
394 const ElementContext& elemCtx OPM_UNUSED,
395 unsigned dofIdx OPM_UNUSED,
396 unsigned timeIdx OPM_UNUSED)
const 398 OPM_THROW(std::logic_error,
399 "Not implemented: The local residual " << Dune::className<Implementation>()
400 <<
" does not implement the required method 'computeSource()'");
408 const ElementContext& elemCtx,
409 unsigned timeIdx)
const 411 if (!elemCtx.onBoundary())
414 BoundaryContext boundaryCtx(elemCtx);
417 size_t numBoundaryFaces = boundaryCtx.numBoundaryFaces(0);
418 for (
unsigned faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx) {
429 size_t numDof = elemCtx.numDof(0);
430 for (
unsigned i=0; i < numDof; i++) {
431 for (
unsigned j = 0; j < numEq; ++ j) {
432 assert(std::isfinite(Toolbox::value(
residual[i][j])));
433 Opm::Valgrind::CheckDefined(
residual[i][j]);
445 const BoundaryContext& boundaryCtx,
446 unsigned boundaryFaceIdx,
447 unsigned timeIdx)
const 449 BoundaryRateVector values;
451 Opm::Valgrind::SetUndefined(values);
452 boundaryCtx.problem().boundary(values, boundaryCtx, boundaryFaceIdx, timeIdx);
453 Opm::Valgrind::CheckDefined(values);
455 const auto& stencil = boundaryCtx.stencil(timeIdx);
456 unsigned dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex();
457 const auto& insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx);
458 for (
unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx)
460 stencil.boundaryFace(boundaryFaceIdx).area()
461 * insideIntQuants.extrusionFactor();
463 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
464 residual[dofIdx][eqIdx] += values[eqIdx];
473 const ElementContext& elemCtx)
const 477 RateVector sourceRate;
483 size_t numPrimaryDof = elemCtx.numPrimaryDof(0);
484 for (
unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
485 Scalar extrusionFactor =
486 elemCtx.intensiveQuantities(dofIdx, 0).extrusionFactor();
488 elemCtx.stencil(0).subControlVolume(dofIdx).volume() * extrusionFactor;
489 Opm::Valgrind::CheckDefined(scvVolume);
494 if (!extensiveStorageTerm &&
495 !std::is_same<Scalar, Evaluation>::value &&
496 dofIdx != elemCtx.focusDofIndex())
498 asImp_().computeStorage(tmp2, elemCtx, dofIdx, 0);
499 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
500 tmp[eqIdx] = tmp2[eqIdx];
503 asImp_().computeStorage(tmp, elemCtx, dofIdx, 0);
504 Opm::Valgrind::CheckDefined(tmp);
506 if (elemCtx.enableStorageCache()) {
507 const auto& model = elemCtx.model();
508 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
509 if (model.newtonMethod().numIterations() == 0 &&
510 !elemCtx.haveStashedIntensiveQuantities())
518 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
519 tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]);
520 Opm::Valgrind::CheckDefined(tmp2);
522 model.updateCachedStorage(globalDofIdx, 1, tmp2);
527 tmp2 = model.cachedStorage(globalDofIdx, 1);
528 Opm::Valgrind::CheckDefined(tmp2);
535 asImp_().computeStorage(tmp2, elemCtx, dofIdx, 1);
536 Opm::Valgrind::CheckDefined(tmp2);
540 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
541 tmp[eqIdx] -= tmp2[eqIdx];
542 tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
544 residual[dofIdx][eqIdx] += tmp[eqIdx];
547 Opm::Valgrind::CheckDefined(
residual[dofIdx]);
550 asImp_().computeSource(sourceRate, elemCtx, dofIdx, 0);
555 if (!extensiveStorageTerm &&
556 !std::is_same<Scalar, Evaluation>::value &&
557 dofIdx != elemCtx.focusDofIndex())
559 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
560 residual[dofIdx][eqIdx] -= Opm::scalarValue(sourceRate[eqIdx])*scvVolume;
563 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
564 sourceRate[eqIdx] *= scvVolume;
565 residual[dofIdx][eqIdx] -= sourceRate[eqIdx];
569 Opm::Valgrind::CheckDefined(
residual[dofIdx]);
574 size_t numDof = elemCtx.numDof(0);
575 for (
unsigned i=0; i < numDof; i++) {
576 for (
unsigned j = 0; j < numEq; ++ j) {
577 assert(Opm::isfinite(
residual[i][j]));
578 Opm::Valgrind::CheckDefined(
residual[i][j]);
586 Implementation& asImp_()
587 {
return *
static_cast<Implementation*
>(
this); }
589 const Implementation& asImp_()
const 590 {
return *
static_cast<const Implementation*
>(
this); }
592 LocalEvalBlockVector internalResidual_;
void eval(LocalEvalBlockVector &residual, const ElementContext &elemCtx) const
Compute the local residual, i.e.
Definition: fvbaselocalresidual.hh:161
This is a stand-alone version of boost::alignment::aligned_allocator from Boost 1...
void eval(const Problem &problem, const Element &element)
Compute the local residual, i.e.
Definition: fvbaselocalresidual.hh:131
Definition: baseauxiliarymodule.hh:37
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
This file provides the infrastructure to retrieve run-time parameters.
Declare the properties used by the infrastructure code of the finite volume discretizations.
void evalFluxes(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Add the flux term to a local residual.
Definition: fvbaselocalresidual.hh:292
const LocalEvalBlockVector & residual() const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:109
void computeFlux(RateVector &flux OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned scvfIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume...
Definition: fvbaselocalresidual.hh:377
void computeStorage(EqVector &storage OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned dofIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the amount all conservation quantities (e.g.
Definition: fvbaselocalresidual.hh:360
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition: fvbaselocalresidual.hh:60
const EvalVector & residual(unsigned dofIdx) const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:118
void evalBoundary_(LocalEvalBlockVector &residual, const ElementContext &elemCtx, unsigned timeIdx) const
Evaluate the boundary conditions of an element.
Definition: fvbaselocalresidual.hh:407
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:207
void evalVolumeTerms_(LocalEvalBlockVector &residual, const 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:472
void computeSource(RateVector &source OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned dofIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Calculate the source term of the equation.
Definition: fvbaselocalresidual.hh:393
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:444
void eval(const ElementContext &elemCtx)
Compute the local residual, i.e.
Definition: fvbaselocalresidual.hh:147
static void registerParameters()
Register all run-time parameters for the local residual.
Definition: fvbaselocalresidual.hh:102