27 #ifndef EWOMS_STOKES_INTENSIVE_QUANTITIES_HH
28 #define EWOMS_STOKES_INTENSIVE_QUANTITIES_HH
33 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
34 #include <dune/geometry/quadraturerules.hh>
36 #include <dune/common/fvector.hh>
48 template <
class TypeTag>
53 typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
54 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
55 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
56 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
57 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
58 typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
59 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
61 enum { numComponents = FluidSystem::numComponents };
62 enum { dim = GridView::dimension };
63 enum { dimWorld = GridView::dimensionworld };
64 enum { pressureIdx = Indices::pressureIdx };
65 enum { moleFrac1Idx = Indices::moleFrac1Idx };
69 typedef typename GridView::ctype CoordScalar;
70 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
71 typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
78 void update(
const ElementContext &elemCtx,
int dofIdx,
int timeIdx)
80 ParentType::update(elemCtx, dofIdx, timeIdx);
82 EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
84 const auto &priVars = elemCtx.primaryVars(dofIdx, timeIdx);
85 fluidState_.setPressure(phaseIdx, priVars[pressureIdx]);
86 Valgrind::CheckDefined(fluidState_.pressure(phaseIdx));
92 fluidState_.setSaturation(phaseIdx, 1.0);
96 for (
int compIdx = 1; compIdx < numComponents; ++compIdx) {
97 fluidState_.setMoleFraction(phaseIdx, compIdx,
98 priVars[moleFrac1Idx + compIdx - 1]);
99 sumx += priVars[moleFrac1Idx + compIdx - 1];
101 fluidState_.setMoleFraction(phaseIdx, 0, 1 - sumx);
104 typename FluidSystem::ParameterCache paramCache;
105 paramCache.updateAll(fluidState_);
107 fluidState_.setDensity(phaseIdx, FluidSystem::density(fluidState_, paramCache, phaseIdx));
108 fluidState_.setViscosity(phaseIdx, FluidSystem::viscosity(fluidState_, paramCache, phaseIdx));
111 EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
114 for (
int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
115 velocityCenter_[dimIdx] = priVars[Indices::velocity0Idx + dimIdx];
119 gravity_ = elemCtx.problem().gravity();
130 for (
int i = 0; i < elemCtx.numDof(0); ++i) {
131 const auto &feGrad = elemCtx.stencil(timeIdx).subControlVolume(dofIdx).gradCenter[i];
132 Valgrind::CheckDefined(feGrad);
133 DimVector tmp(feGrad);
134 tmp *= elemCtx.intensiveQuantities(i, timeIdx).fluidState().pressure(phaseIdx);
135 Valgrind::CheckDefined(tmp);
137 pressureGrad_ += tmp;
142 const auto &stencil = elemCtx.stencil(timeIdx);
143 const auto &scvLocalGeom = stencil.subControlVolume(dofIdx).localGeometry();
145 Dune::GeometryType geomType = scvLocalGeom.type();
146 static const int quadratureOrder = 2;
147 const auto &rule = Dune::QuadratureRules<Scalar, dimWorld>::rule(geomType, quadratureOrder);
151 for (
auto it = rule.begin(); it != rule.end(); ++it) {
152 const auto &posScvLocal = it->position();
153 const auto &posElemLocal = scvLocalGeom.global(posScvLocal);
155 DimVector velocityAtPos = velocityAtPos_(elemCtx, timeIdx, posElemLocal);
156 Scalar weight = it->weight();
160 velocity_.axpy(weight * detjac, velocityAtPos);
173 {
return fluidState_; }
190 {
return velocity_; }
196 {
return velocityCenter_; }
202 {
return pressureGrad_; }
212 DimVector velocityAtPos_(
const ElementContext &elemCtx,
214 const LocalPosition &localPos)
const
217 elemCtx.gradientCalculator().localFiniteElementCache();
218 const auto &localFiniteElement =
219 feCache.get(elemCtx.element().type());
221 typedef Dune::FieldVector<Scalar, 1> ShapeValue;
222 std::vector<ShapeValue> shapeValue;
224 localFiniteElement.localBasis().evaluateFunction(localPos, shapeValue);
226 DimVector result(0.0);
227 for (
int dofIdx = 0; dofIdx < elemCtx.numDof(0); dofIdx++) {
228 result.axpy(shapeValue[dofIdx][0], elemCtx.intensiveQuantities(dofIdx, timeIdx).velocityCenter());
235 DimVector velocityCenter_;
237 DimVector pressureGrad_;
238 FluidState fluidState_;
const DimVector & gravity() const
Returns the gravitational acceleration vector in the sub-control volume.
Definition: stokesintensivequantities.hh:208
Provides the volumetric quantities required for the energy equation.
Definition: energymodule.hh:530
Contains the intensive quantities of the Stokes model.
Definition: stokesintensivequantities.hh:49
Declares the properties required by the Stokes model.
Scalar porosity() const
Returns the porosity of the medium.
Definition: stokesintensivequantities.hh:183
const FluidState & fluidState() const
Returns the thermodynamic state of the fluid for the control-volume.
Definition: stokesintensivequantities.hh:172
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:485
void update(const ElementContext &elemCtx, int dofIdx, int timeIdx)
Definition: stokesintensivequantities.hh:78
Definition: baseauxiliarymodule.hh:35
const DimVector & velocity() const
Returns the average velocity in the sub-control volume.
Definition: stokesintensivequantities.hh:189
const DimVector & pressureGradient() const
Returns the pressure gradient in the sub-control volume.
Definition: stokesintensivequantities.hh:201
void updateScvGradients(const ElementContext &elemCtx, int dofIdx, int timeIdx)
Definition: stokesintensivequantities.hh:125
const DimVector & velocityCenter() const
Returns the velocity at the center in the sub-control volume.
Definition: stokesintensivequantities.hh:195
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...