26 #ifndef EWOMS_FV_BASE_LOCAL_RESIDUAL_HH
27 #define EWOMS_FV_BASE_LOCAL_RESIDUAL_HH
29 #include <opm/material/common/Valgrind.hpp>
31 #include <dune/istl/bvector.hh>
32 #include <dune/grid/common/geometry.hh>
34 #include <dune/common/fvector.hh>
36 #include <opm/material/common/ClassName.hpp>
51 template<
class TypeTag>
55 typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
57 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
58 typedef typename GridView::template Codim<0>::Entity Element;
60 typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
61 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
62 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
63 typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
64 typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
65 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
66 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
67 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
70 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
71 typedef typename GET_PROP_TYPE(TypeTag, BoundaryContext) BoundaryContext;
72 typedef typename GET_PROP_TYPE(TypeTag, ConstraintsContext) ConstraintsContext;
74 typedef Opm::MathToolbox<Evaluation> Toolbox;
75 typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
76 typedef Dune::FieldVector<Evaluation, numEq> VectorBlock;
77 typedef Dune::BlockVector<VectorBlock> LocalBlockVector;
101 {
return internalResidual_; }
110 {
return internalResidual_[dofIdx]; }
117 {
return internalStorageTerm_; }
126 {
return internalStorageTerm_[dofIdx]; }
139 void eval(
const Problem &problem,
const Element &element)
141 ElementContext elemCtx(problem);
142 elemCtx.updateAll(element);
156 void eval(
const ElementContext &elemCtx)
158 int numDof = elemCtx.numDof(0);
159 int numPrimaryDof = elemCtx.numPrimaryDof(0);
160 internalResidual_.resize(numDof);
161 internalStorageTerm_.resize(numPrimaryDof);
162 asImp_().eval(internalResidual_, internalStorageTerm_, elemCtx);
174 LocalBlockVector &storage,
175 const ElementContext &elemCtx)
const
177 int numPrimaryDof = elemCtx.numPrimaryDof(0);
178 int numDof = elemCtx.numDof(0);
180 typedef Opm::MathToolbox<Evaluation> Toolbox;
182 residual = Toolbox::createConstant(0.0);
183 storage = Toolbox::createConstant(0.0);
186 asImp_().evalFluxes(residual, elemCtx, 0);
187 assert(static_cast<int>(residual.size()) == numDof);
188 assert(static_cast<int>(storage.size()) == numPrimaryDof);
191 for (
int i=0; i < numDof; i++) {
192 for (
int j = 0; j < numEq; ++ j) {
193 assert(std::isfinite(Toolbox::value(residual[i][j])));
194 Valgrind::CheckDefined(residual[i][j]);
200 asImp_().evalVolumeTerms_(residual, storage, elemCtx);
201 for (
int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
202 for (
int eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
203 storage[dofIdx][eqIdx] /= elemCtx.dofTotalVolume(dofIdx, 0);
205 assert(std::isfinite(Toolbox::value(residual[dofIdx][eqIdx])));
206 Valgrind::CheckDefined(residual[dofIdx][eqIdx]);
211 asImp_().evalBoundary_(residual, elemCtx, 0);
213 for (
int i=0; i < numDof; i++) {
214 for (
int j = 0; j < numEq; ++ j) {
215 assert(std::isfinite(Toolbox::value(residual[i][j])));
216 Valgrind::CheckDefined(residual[i][j]);
222 asImp_().evalConstraints_(residual, storage, elemCtx, 0);
226 for (
int dofIdx=0; dofIdx < numDof; ++dofIdx) {
227 if (elemCtx.dofTotalVolume(dofIdx, 0) > 0) {
229 Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, 0);
231 for (
int eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
232 residual[dofIdx][eqIdx] /= dofVolume;
233 assert(std::isfinite(Toolbox::value(residual[dofIdx][eqIdx])));
234 Valgrind::CheckDefined(residual[dofIdx][eqIdx]);
252 int numPrimaryDof = elemCtx.numPrimaryDof(0);
253 internalStorageTerm_.resize(numPrimaryDof);
271 const ElementContext &elemCtx,
280 int numPrimaryDof = elemCtx.numPrimaryDof(0);
281 for (
int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
282 storage[dofIdx] = Toolbox::createConstant(0.0);
283 asImp_().computeStorage(storage[dofIdx], elemCtx, dofIdx, timeIdx);
285 for (
int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
286 storage[dofIdx][eqIdx] *=
287 elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
288 * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
298 Dune::FieldVector<Scalar, numEq> tmp;
299 int numPrimaryDof = elemCtx.numPrimaryDof(0);
300 for (
int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
302 asImp_().computeStorage(tmp,
307 elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
308 * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
311 for (
int eqIdx = 0; eqIdx < numEq; ++eqIdx)
312 storage[dofIdx][eqIdx] = Toolbox::createConstant(tmp[eqIdx]);
325 const ElementContext &elemCtx,
330 const auto &stencil = elemCtx.stencil(timeIdx);
332 int numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
333 for (
int scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
334 const auto &face = stencil.interiorFace(scvfIdx);
335 int i = face.interiorIndex();
336 int j = face.exteriorIndex();
338 Valgrind::SetUndefined(flux);
339 asImp_().computeFlux(flux, elemCtx, scvfIdx, timeIdx);
341 Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor();
342 alpha *= face.area();
343 for (
int eqIdx = 0; eqIdx < numEq; ++ eqIdx)
344 flux[eqIdx] *= alpha;
345 Valgrind::CheckDefined(flux);
378 const ElementContext &elemCtx,
382 OPM_THROW(std::logic_error,
383 "Not implemented: The local residual " << Opm::className<Implementation>()
384 <<
" does not implement the required method 'computeStorage()'");
395 const ElementContext &elemCtx,
399 OPM_THROW(std::logic_error,
400 "Not implemented: The local residual " << Opm::className<Implementation>()
401 <<
" does not implement the required method 'computeFlux()'");
411 const ElementContext &elemCtx,
415 OPM_THROW(std::logic_error,
416 "Not implemented: The local residual " << Opm::className<Implementation>()
417 <<
" does not implement the required method 'computeSource()'");
425 const ElementContext &elemCtx,
428 if (!elemCtx.onBoundary())
431 BoundaryContext boundaryCtx(elemCtx);
434 int numBoundaryFaces = boundaryCtx.numBoundaryFaces(0);
435 for (
int faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx) {
450 const BoundaryContext &boundaryCtx,
454 BoundaryRateVector values;
456 Valgrind::SetUndefined(values);
457 boundaryCtx.problem().boundary(values,
461 Valgrind::CheckDefined(values);
463 const auto &stencil = boundaryCtx.stencil(timeIdx);
464 int dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex();
465 const auto &insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx);
466 for (
unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx)
468 stencil.boundaryFace(boundaryFaceIdx).area()
469 * insideIntQuants.extrusionFactor();
471 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
472 residual[dofIdx][eqIdx] += values[eqIdx];
480 LocalBlockVector &storage,
481 const ElementContext &elemCtx,
487 const auto &problem = elemCtx.problem();
488 Constraints constraints;
489 ConstraintsContext constraintsCtx(elemCtx);
490 int numPrimaryDof = elemCtx.numPrimaryDof(0);
491 for (
int dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
494 problem.constraints(constraints, elemCtx, dofIdx, timeIdx);
496 if (!constraints.isConstraint())
500 const PrimaryVariables &priVars = elemCtx.primaryVars(dofIdx, timeIdx);
501 for (
int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
502 if (!constraints.isConstraint(eqIdx))
505 int pvIdx = constraints.eqToPvIndex(eqIdx);
507 assert(0 <= pvIdx && pvIdx < numEq);
508 Valgrind::CheckDefined(constraints[pvIdx]);
510 residual[dofIdx][eqIdx] = priVars.makeEvaluation(pvIdx, timeIdx) - constraints[pvIdx];
511 storage[dofIdx][eqIdx] = 0.0;
522 LocalBlockVector &storage,
523 const ElementContext &elemCtx)
const
527 RateVector sourceRate;
529 tmp = Toolbox::createConstant(0.0);
533 int numPrimaryDof = elemCtx.numPrimaryDof(0);
534 for (
int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
535 Scalar extrusionFactor =
536 elemCtx.intensiveQuantities(dofIdx, 0).extrusionFactor();
538 elemCtx.stencil(0).subControlVolume(dofIdx).volume() * extrusionFactor;
539 Valgrind::CheckDefined(scvVolume);
547 asImp_().computeStorage(tmp,
551 Valgrind::CheckDefined(tmp);
553 asImp_().computeStorage(tmp2,
557 Valgrind::CheckDefined(tmp2);
560 for (
int eqIdx = 0; eqIdx < numEq; eqIdx++)
561 tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
563 storage[dofIdx] += tmp;
564 residual[dofIdx] += tmp;
566 Valgrind::CheckDefined(storage[dofIdx]);
567 Valgrind::CheckDefined(residual[dofIdx]);
570 asImp_().computeSource(sourceRate, elemCtx, dofIdx, 0);
571 for (
int eqIdx = 0; eqIdx < numEq; ++eqIdx)
572 sourceRate[eqIdx] *= scvVolume;
573 residual[dofIdx] -= sourceRate;
577 Valgrind::CheckDefined(storage[dofIdx]);
578 Valgrind::CheckDefined(residual[dofIdx]);
584 Implementation &asImp_()
586 assert(static_cast<Implementation*>(
this) != 0);
587 return *
static_cast<Implementation*
>(
this);
590 const Implementation &asImp_()
const
592 assert(static_cast<const Implementation*>(
this) != 0);
593 return *
static_cast<const Implementation*
>(
this);
596 LocalBlockVector internalResidual_;
597 LocalBlockVector internalStorageTerm_;
~FvBaseLocalResidual()
Definition: fvbaselocalresidual.hh:87
void eval(LocalBlockVector &residual, LocalBlockVector &storage, const ElementContext &elemCtx) const
Compute the local residual, i.e. the deviation of the conservation equations from zero...
Definition: fvbaselocalresidual.hh:173
void eval(const ElementContext &elemCtx)
Compute the local residual, i.e. the deviation of the conservation equations from zero and store the ...
Definition: fvbaselocalresidual.hh:156
Declare the properties used by the infrastructure code of the finite volume discretizations.
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
const LocalBlockVector & storageTerm() const
Return the storage term calculated using the last call to eval() using internal storage.
Definition: fvbaselocalresidual.hh:116
void evalConstraints_(LocalBlockVector &residual, LocalBlockVector &storage, const ElementContext &elemCtx, int timeIdx) const
Set the values of the constraint volumes of the current element.
Definition: fvbaselocalresidual.hh:479
void evalVolumeTerms_(LocalBlockVector &residual, LocalBlockVector &storage, 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:521
void computeSource(RateVector &source, const ElementContext &elemCtx, int dofIdx, int timeIdx) const
Calculate the source term of the equation.
Definition: fvbaselocalresidual.hh:410
void computeFlux(RateVector &flux, const ElementContext &elemCtx, int scvfIdx, int timeIdx) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume...
Definition: fvbaselocalresidual.hh:394
void eval(const Problem &problem, const Element &element)
Compute the local residual, i.e. the deviation of the conservation equations from zero and store the ...
Definition: fvbaselocalresidual.hh:139
void computeStorage(EqVector &storage, const ElementContext &elemCtx, int dofIdx, int timeIdx) const
Evaluate the amount all conservation quantities (e.g. phase mass) within a finite sub-control volume...
Definition: fvbaselocalresidual.hh:377
const VectorBlock & storageTerm(int dofIdx) const
Return the storage term calculated using the last call to eval() using internal storage.
Definition: fvbaselocalresidual.hh:125
Definition: baseauxiliarymodule.hh:35
const LocalBlockVector & residual() const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:100
void evalFluxes(LocalBlockVector &residual, const ElementContext &elemCtx, int timeIdx) const
Add the flux term to a local residual.
Definition: fvbaselocalresidual.hh:324
const VectorBlock & residual(int dofIdx) const
Return the result of the eval() call using internal storage.
Definition: fvbaselocalresidual.hh:109
void evalStorage(const ElementContext &elemCtx, int timeIdx)
Calculate the amount of all conservation quantities stored in all element's primary sub-control volum...
Definition: fvbaselocalresidual.hh:250
void evalBoundary_(LocalBlockVector &residual, const ElementContext &elemCtx, int timeIdx) const
Evaluate the boundary conditions of an element.
Definition: fvbaselocalresidual.hh:424
FvBaseLocalResidual()
Definition: fvbaselocalresidual.hh:84
void evalBoundarySegment_(LocalBlockVector &residual, const BoundaryContext &boundaryCtx, int boundaryFaceIdx, int timeIdx) const
Evaluate all boundary conditions for a single sub-control volume face to the local residual...
Definition: fvbaselocalresidual.hh:449
Element-wise caculation of the residual matrix for models based on a finite volume spatial discretiza...
Definition: fvbaselocalresidual.hh:52
static void registerParameters()
Register all run-time parameters for the local residual.
Definition: fvbaselocalresidual.hh:93
void evalStorage(LocalBlockVector &storage, const ElementContext &elemCtx, int timeIdx) const
Calculate the amount of all conservation quantities stored in all element's sub-control volumes for a...
Definition: fvbaselocalresidual.hh:270