25#ifndef DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
26#define DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
30#include <dune/common/fvector.hh>
31#include <dune/common/float_cmp.hh>
32#include <dune/geometry/referenceelements.hh>
45template<
class Gr
idVariables,
class FluxVariables>
48 using GridGeometry =
typename GridVariables::GridGeometry;
49 using FVElementGeometry =
typename GridGeometry::LocalView;
50 using SubControlVolume =
typename GridGeometry::SubControlVolume;
51 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
52 using GridView =
typename GridGeometry::GridView;
53 using Element =
typename GridView::template Codim<0>::Entity;
54 using GridVolumeVariables =
typename GridVariables::GridVolumeVariables;
55 using ElementFluxVarsCache =
typename GridVariables::GridFluxVariablesCache::LocalView;
56 using VolumeVariables =
typename GridVariables::VolumeVariables;
57 using ElementVolumeVariables =
typename GridVolumeVariables::LocalView;
58 using FluidSystem =
typename VolumeVariables::FluidSystem;
59 using Scalar =
typename GridVariables::Scalar;
61 static constexpr int dim = GridView::dimension;
62 static constexpr int dimWorld = GridView::dimensionworld;
64 static constexpr int dofCodim = isBox ? dim : 0;
66 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
67 using ReferenceElements = Dune::ReferenceElements<typename GridView::ctype, dim>;
69 using Problem =
typename GridVolumeVariables::Problem;
70 using BoundaryTypes =
typename Problem::Traits::BoundaryTypes;
82 : problem_(gridVariables.curGridVolVars().problem())
83 , gridGeometry_(gridVariables.gridGeometry())
84 , gridVariables_(gridVariables)
90 cellNum_.assign(gridGeometry_.gridView().size(dim), 0);
92 for (
const auto& element : elements(gridGeometry_.gridView()))
93 for (
unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
94 ++cellNum_[gridGeometry_.vertexMapper().subIndex(element, vIdx, dim)];
101 const Element& element,
102 const FVElementGeometry& fvGeometry,
103 const ElementVolumeVariables& elemVolVars,
104 const ElementFluxVarsCache& elemFluxVarsCache,
107 using Velocity =
typename VelocityVector::value_type;
109 const auto geometry = element.geometry();
110 const Dune::GeometryType geomType = geometry.type();
113 auto upwindTerm = [phaseIdx](
const auto& volVars) {
return volVars.mobility(phaseIdx); };
115 if(isBox && dim == 1)
118 tmpVelocity = (geometry.corner(1) - geometry.corner(0));
119 tmpVelocity /= tmpVelocity.two_norm();
121 for (
auto&& scvf : scvfs(fvGeometry))
127 FluxVariables fluxVars;
128 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
132 Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
133 Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
134 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
135 flux /= insideVolVars.extrusionFactor();
138 const int eIdxGlobal = gridGeometry_.elementMapper().index(element);
139 velocity[eIdxGlobal] = tmpVelocity;
145 const auto referenceElement = ReferenceElements::general(geomType);
146 const auto& localPos = referenceElement.position(0, 0);
147 const auto jacobianT2 = geometry.jacobianTransposed(localPos);
151 using ScvVelocities = Dune::BlockVector<Velocity>;
152 ScvVelocities scvVelocities(fvGeometry.numScv());
155 for (
auto&& scvf : scvfs(fvGeometry))
161 const auto localPosIP = geometry.local(scvf.ipGlobal());
164 const auto jacobianT1 = geometry.jacobianTransposed(localPosIP);
165 const auto globalNormal = scvf.unitOuterNormal();
166 GlobalPosition localNormal(0);
167 jacobianT1.mv(globalNormal, localNormal);
168 localNormal /= localNormal.two_norm();
171 FluxVariables fluxVars;
172 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
176 Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
177 Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
178 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
179 flux /= insideVolVars.extrusionFactor();
185 scvVelocities[scvf.insideScvIdx()] += tmpVelocity;
186 scvVelocities[scvf.outsideScvIdx()] += tmpVelocity;
190 for (
auto&& scv : scvs(fvGeometry))
192 int vIdxGlobal = scv.dofIndex();
197 jacobianT2.mtv(scvVelocities[scv.indexInElement()], scvVelocity);
198 scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal];
200 velocity[vIdxGlobal] += scvVelocity;
209 const int numScvfsPerFace = isMpfa ? element.template subEntity<1>(0).geometry().corners() : 1;
211 if (fvGeometry.numScvf() != element.subEntities(1)*numScvfsPerFace)
212 DUNE_THROW(Dune::NotImplemented,
"Velocity output for non-conforming grids");
214 if (!geomType.isCube() && !geomType.isSimplex())
215 DUNE_THROW(Dune::NotImplemented,
"Velocity output for other geometry types than cube and simplex");
221 std::vector<bool> handledScvf;
223 handledScvf.resize(element.subEntities(1),
false);
226 std::vector<unsigned int> scvfIndexInInside(fvGeometry.numScvf());
227 int localScvfIdx = 0;
228 for (
const auto& intersection :
intersections(gridGeometry_.gridView(), element))
231 if (handledScvf[intersection.indexInInside()])
234 if (intersection.neighbor() || intersection.boundary())
236 for (
int i = 0; i < numScvfsPerFace; ++i)
237 scvfIndexInInside[localScvfIdx++] = intersection.indexInInside();
241 handledScvf[intersection.indexInInside()] =
true;
245 std::vector<Scalar> scvfFluxes(element.subEntities(1), 0.0);
247 for (
auto&& scvf : scvfs(fvGeometry))
249 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
250 if (!scvf.boundary())
252 FluxVariables fluxVars;
253 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
254 scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor();
258 auto bcTypes = problemBoundaryTypes_(element, scvf);
259 if (bcTypes.hasOnlyDirichlet())
261 FluxVariables fluxVars;
262 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
263 scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor();
278 for (
auto&& scvf : scvfs(fvGeometry))
282 auto bcTypes = problemBoundaryTypes_(element, scvf);
283 if (bcTypes.hasNeumann())
286 const auto neumannFlux = Deprecated::neumann(problem_, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
287 using NumEqVector = std::decay_t<
decltype(neumannFlux)>;
288 if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30))
289 scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
291 else if (dim == 1 || geomType.isCube())
293 const auto fIdx = scvfIndexInInside[localScvfIdx];
294 const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
295 scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
298 else if (geomType.isSimplex())
299 scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
310 if (dim == 1 || geomType.isCube())
312 for (
int i = 0; i < dim; i++)
313 refVelocity[i] = 0.5 * (scvfFluxes[2*i + 1] - scvfFluxes[2*i]);
316 else if (geomType.isSimplex())
318 for (
int dimIdx = 0; dimIdx < dim; dimIdx++)
320 refVelocity[dimIdx] = -scvfFluxes[dim - 1 - dimIdx];
321 for (
int fIdx = 0; fIdx < dim + 1; fIdx++)
322 refVelocity[dimIdx] += scvfFluxes[fIdx]/(dim + 1);
327 DUNE_THROW(Dune::NotImplemented,
"velocity output for cell-centered and prism/pyramid");
330 jacobianT2.mtv(refVelocity, scvVelocity);
332 scvVelocity /= geometry.integrationElement(localPos);
334 int eIdxGlobal = gridGeometry_.elementMapper().index(element);
336 velocity[eIdxGlobal] = scvVelocity;
345 static Scalar scvfReferenceArea_(Dune::GeometryType geomType,
int fIdx)
347 if (dim == 1 || geomType.isCube())
349 return 1.0/(1 << (dim-1));
351 else if (geomType.isTriangle())
353 static const Scalar faceToArea[] = {0.372677996249965,
356 return faceToArea[fIdx];
358 else if (geomType.isTetrahedron())
360 static const Scalar faceToArea[] = {0.102062072615966,
366 return faceToArea[fIdx];
368 else if (geomType.isPyramid())
370 static const Scalar faceToArea[] = {0.130437298687488,
378 return faceToArea[fIdx];
380 else if (geomType.isPrism())
382 static const Scalar faceToArea[] = {0.166666666666667,
391 return faceToArea[fIdx];
394 DUNE_THROW(Dune::NotImplemented,
"scvf area for unknown GeometryType");
413 template <
bool enable = isBox,
typename std::enable_if_t<!enable,
int> = 0>
414 BoundaryTypes problemBoundaryTypes_(
const Element& element,
const SubControlVolumeFace& scvf)
const
415 {
return problem_.boundaryTypes(element, scvf); }
418 template <
bool enable = isBox,
typename std::enable_if_t<enable,
int> = 0>
419 BoundaryTypes problemBoundaryTypes_(
const Element& element,
const SubControlVolumeFace& scvf)
const
420 {
return BoundaryTypes(); }
422 const Problem& problem_;
423 const GridGeometry& gridGeometry_;
424 const GridVariables& gridVariables_;
426 std::vector<int> cellNum_;
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Element solution classes and factory functions.
The available discretization methods in Dumux.
Dune::IteratorRange< typename MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper >::Intersections::const_iterator > intersections(const MultiDomainGlue< DomainGridView, TargetGridView, DomainMapper, TargetMapper > &glue)
Range generator to iterate with range-based for loops over all intersections as follows: for (const a...
Definition: glue.hh:62
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
Property tag Velocity
The type velocity reconstruction.
Definition: porousmediumflow/sequential/properties.hh:67
Velocity computation for implicit (porous media) models.
Definition: velocity.hh:47
void calculateVelocity(VelocityVector &velocity, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const ElementFluxVarsCache &elemFluxVarsCache, int phaseIdx) const
Definition: velocity.hh:100
PorousMediumFlowVelocity(const GridVariables &gridVariables)
Constructor initializes the static data with the initial solution.
Definition: velocity.hh:81
std::vector< Dune::FieldVector< Scalar, dimWorld > > VelocityVector
Definition: velocity.hh:74
static constexpr int numFluidPhases
Definition: velocity.hh:73