GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/porousmediumflow/velocity.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 111 137 81.0%
Functions: 110 454 24.2%
Branches: 141 404 34.9%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup PorousmediumflowModels
10 * \brief Velocity computation for implicit (porous media) models.
11 */
12
13 #ifndef DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
14 #define DUMUX_POROUSMEDIUMFLOW_VELOCITY_HH
15
16 #include <vector>
17 #include <type_traits>
18
19 #include <dune/common/fvector.hh>
20 #include <dune/common/float_cmp.hh>
21 #include <dune/geometry/type.hh>
22 #include <dune/geometry/referenceelements.hh>
23
24 #include <dumux/common/parameters.hh>
25 #include <dumux/discretization/extrusion.hh>
26 #include <dumux/discretization/method.hh>
27 #include <dumux/discretization/elementsolution.hh>
28 #include <dumux/flux/traits.hh>
29
30 namespace Dumux {
31
32 #ifndef DOXYGEN
33 namespace Detail {
34 // helper structs and functions detecting if the model is compositional
35
36 template <typename T, typename ...Ts>
37 using MoleFractionDetector = decltype(std::declval<T>().moleFraction(std::declval<Ts>()...));
38
39 template<class T, typename ...Args>
40 static constexpr bool hasMoleFraction()
41 { return Dune::Std::is_detected<MoleFractionDetector, T, Args...>::value; }
42
43 template <typename T, typename ...Ts>
44 using MassFractionDetector = decltype(std::declval<T>().massFraction(std::declval<Ts>()...));
45
46 template<class T, typename ...Args>
47 static constexpr bool hasMassFraction()
48 { return Dune::Std::is_detected<MassFractionDetector, T, Args...>::value; }
49
50 } // end namespace Detail
51 #endif // DOXYGEN
52
53 /*!
54 * \ingroup PorousmediumflowModels
55 * \brief Velocity computation for implicit (porous media) models.
56 */
57 template<class GridVariables, class FluxVariables>
58 ✗ class PorousMediumFlowVelocity
59 {
60 using GridGeometry = typename GridVariables::GridGeometry;
61 using Extrusion = Extrusion_t<GridGeometry>;
62 using FVElementGeometry = typename GridGeometry::LocalView;
63 using SubControlVolume = typename GridGeometry::SubControlVolume;
64 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
65 using GridView = typename GridGeometry::GridView;
66 using Element = typename GridView::template Codim<0>::Entity;
67 using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
68 using ElementFluxVarsCache = typename GridVariables::GridFluxVariablesCache::LocalView;
69 using VolumeVariables = typename GridVariables::VolumeVariables;
70 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
71 using FluidSystem = typename VolumeVariables::FluidSystem;
72 using Scalar = typename GridVariables::Scalar;
73 using FluxTraits = typename Dumux::FluxTraits<FluxVariables>;
74 using AdvectionType = typename FluxVariables::AdvectionType;
75
76 static constexpr int dim = GridView::dimension;
77 static constexpr int dimWorld = GridView::dimensionworld;
78 static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box;
79 static constexpr bool isDiamond = GridGeometry::discMethod == DiscretizationMethods::fcdiamond;
80 static constexpr bool stationaryVelocityField = FluxTraits::hasStationaryVelocityField();
81
82 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
83
84 using Problem = typename GridVolumeVariables::Problem;
85
86 static constexpr bool modelIsCompositional = Detail::hasMoleFraction<typename GridVolumeVariables::VolumeVariables, int, int>() ||
87 Detail::hasMassFraction<typename GridVolumeVariables::VolumeVariables, int, int>();
88
89 static_assert(VolumeVariables::numFluidPhases() >= 1, "Velocity output only makes sense for models with fluid phases.");
90
91 public:
92 static constexpr int numFluidPhases = VolumeVariables::numFluidPhases();
93 using VelocityVector = std::vector<Dune::FieldVector<Scalar, dimWorld>>;
94
95 /*!
96 * \brief Constructor initializes the static data with the initial solution.
97 *
98 * \param gridVariables The grid variables
99 */
100 95 PorousMediumFlowVelocity(const GridVariables& gridVariables)
101
1/2
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
 98 : problem_(gridVariables.curGridVolVars().problem())
102 , gridGeometry_(gridVariables.gridGeometry())
103
2/4
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
 252 , gridVariables_(gridVariables)
104 {
105 // set the number of scvs the vertices are connected to
106 if constexpr (isBox && dim > 1)
107 {
108 // resize to the number of vertices of the grid
109
4/10
✓ Branch 1 taken 21 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 21 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 21 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 66 cellNum_.assign(gridGeometry_.gridView().size(dim), 0);
110
111
7/19
✓ Branch 1 taken 33 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 33 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 16749 times.
✓ Branch 7 taken 12 times.
✓ Branch 8 taken 2878 times.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2878 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 16794 for (const auto& element : elements(gridGeometry_.gridView()))
112
5/5
✓ Branch 0 taken 55400 times.
✓ Branch 1 taken 28240 times.
✓ Branch 2 taken 55400 times.
✓ Branch 3 taken 25362 times.
✓ Branch 4 taken 2878 times.
 152890 for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
113
2/4
✓ Branch 1 taken 66912 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 66912 times.
✗ Branch 5 not taken.
 133824 ++cellNum_[gridGeometry_.vertexMapper().subIndex(element, vIdx, dim)];
114 }
115 33 }
116
117 //! Calculates the velocities for the scvs in the element.
118 //! We assume the local containers to be bound to the complete stencil.
119 1198080 void calculateVelocity(VelocityVector& velocity,
120 const Element& element,
121 const FVElementGeometry& fvGeometry,
122 const ElementVolumeVariables& elemVolVars,
123 const ElementFluxVarsCache& elemFluxVarsCache,
124 int phaseIdx) const
125 {
126 using Velocity = typename VelocityVector::value_type;
127
128 1198080 const auto geometry = element.geometry();
129
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 1198080 const Dune::GeometryType geomType = geometry.type();
130
131 // the upwind term to be used for the volume flux evaluation
132
6/14
✓ Branch 0 taken 174854 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 174854 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 174854 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 131746 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 131746 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 131746 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 15711540 auto upwindTerm = [phaseIdx](const auto& volVars) { return volVars.mobility(phaseIdx); };
133
134 if constexpr (isBox && dim == 1)
135 {
136 2960 Velocity tmpVelocity(0.0);
137 8880 tmpVelocity = (geometry.corner(1) - geometry.corner(0));
138 5920 tmpVelocity /= tmpVelocity.two_norm();
139
140
4/4
✓ Branch 0 taken 2996 times.
✓ Branch 1 taken 2960 times.
✓ Branch 2 taken 2996 times.
✓ Branch 3 taken 2960 times.
 8916 for (auto&& scvf : scvfs(fvGeometry))
141 {
142
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 2960 times.
 2996 if (scvf.boundary())
143 36 continue;
144
145 // instantiate the flux variables
146 2960 FluxVariables fluxVars;
147 2960 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
148
149 // get the volume flux divided by the area of the
150 // subcontrolvolume face in the reference element
151 2960 Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
152 2960 Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
153 5920 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
154 2960 flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area();
155 2960 tmpVelocity *= flux;
156
157 5920 const int eIdxGlobal = gridGeometry_.elementMapper().index(element);
158 5920 velocity[eIdxGlobal] = tmpVelocity;
159 }
160 return;
161 }
162
163 // get the transposed Jacobian of the element mapping
164 using Dune::referenceElement;
165
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 1195120 const auto refElement = referenceElement(geometry);
166
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 800 times.
 1195120 const auto& localPos = refElement.position(0, 0);
167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 800 times.
 1195120 const auto jacobianT2 = geometry.jacobianTransposed(localPos);
168
169 if constexpr (isBox)
170 {
171 using ScvVelocities = Dune::BlockVector<Velocity>;
172
1/9
✗ Branch 1 not taken.
✓ Branch 2 taken 349212 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 1047636 ScvVelocities scvVelocities(fvGeometry.numScv());
173 349212 scvVelocities = 0;
174
175
4/4
✓ Branch 0 taken 1728572 times.
✓ Branch 1 taken 349212 times.
✓ Branch 2 taken 1728572 times.
✓ Branch 3 taken 349212 times.
 2426996 for (auto&& scvf : scvfs(fvGeometry))
176 {
177
2/2
✓ Branch 0 taken 331724 times.
✓ Branch 1 taken 1396848 times.
 1728572 if (scvf.boundary())
178 331724 continue;
179
180 // local position of integration point
181
2/4
✓ Branch 1 taken 263888 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 263888 times.
✗ Branch 5 not taken.
 2793696 const auto localPosIP = geometry.local(scvf.ipGlobal());
182
183 // Transformation of the global normal vector to normal vector in the reference element
184
1/2
✓ Branch 1 taken 263888 times.
✗ Branch 2 not taken.
 1396848 const auto jacobianT1 = geometry.jacobianTransposed(localPosIP);
185 1396848 const auto globalNormal = scvf.unitOuterNormal();
186 1396848 GlobalPosition localNormal(0);
187 jacobianT1.mv(globalNormal, localNormal);
188 1396848 localNormal /= localNormal.two_norm();
189
190 // instantiate the flux variables
191
1/2
✓ Branch 1 taken 1042048 times.
✗ Branch 2 not taken.
 1396848 FluxVariables fluxVars;
192
1/2
✓ Branch 1 taken 1396848 times.
✗ Branch 2 not taken.
 1396848 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
193
194 // get the volume flux divided by the area of the
195 // subcontrolvolume face in the reference element
196
1/2
✓ Branch 1 taken 1396848 times.
✗ Branch 2 not taken.
 1396848 Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
197
1/2
✓ Branch 1 taken 1396848 times.
✗ Branch 2 not taken.
 1396848 Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
198 2793696 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
199 1396848 flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area();
200
201 // transform the volume flux into a velocity vector
202 1396848 Velocity tmpVelocity = localNormal;
203 1396848 tmpVelocity *= flux;
204
205 4190544 scvVelocities[scvf.insideScvIdx()] += tmpVelocity;
206
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1396848 times.
 2793696 scvVelocities[scvf.outsideScvIdx()] += tmpVelocity;
207 }
208
209 // transform vertex velocities from local to global coordinates
210
4/4
✓ Branch 0 taken 1396848 times.
✓ Branch 1 taken 349212 times.
✓ Branch 2 taken 1396848 times.
✓ Branch 3 taken 349212 times.
 3492120 for (auto&& scv : scvs(fvGeometry))
211 {
212 1396848 int vIdxGlobal = scv.dofIndex();
213
214 // calculate the subcontrolvolume velocity by the Piola transformation
215 1396848 Velocity scvVelocity(0);
216
217 2793696 jacobianT2.mtv(scvVelocities[scv.indexInElement()], scvVelocity);
218
1/2
✓ Branch 1 taken 263888 times.
✗ Branch 2 not taken.
 2529808 scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal];
219 // add up the wetting phase subcontrolvolume velocities for each vertex
220 4190544 velocity[vIdxGlobal] += scvVelocity;
221 }
222 }
223 else if constexpr (isDiamond)
224 ✗ DUNE_THROW(Dune::NotImplemented, "Velocity output with diamond discretization");
225 else
226 {
227 // For the number of scvfs per facet (mpfa) we simply obtain the number of
228 // corners of the first facet as prisms/pyramids are not supported here anyway
229 // -> for prisms/pyramids the number of scvfs would differ from facet to facet
230 static constexpr bool isMpfa = GridGeometry::discMethod == DiscretizationMethods::ccmpfa;
231
5/11
✗ Branch 1 not taken.
✓ Branch 2 taken 2592 times.
✓ Branch 3 taken 27800 times.
✓ Branch 4 taken 2592 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2592 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2592 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 848500 const int numScvfsPerFace = isMpfa ? element.template subEntity<1>(0).geometry().corners() : 1;
232
233
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 840736 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 845908 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 840736 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2532552 if (fvGeometry.numScvf() != element.subEntities(1)*numScvfsPerFace)
234 ✗ DUNE_THROW(Dune::NotImplemented, "Velocity output for non-conforming grids");
235
236
2/6
✓ Branch 0 taken 5172 times.
✓ Branch 1 taken 840736 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 845908 if (!geomType.isCube() && !geomType.isSimplex())
237 ✗ DUNE_THROW(Dune::NotImplemented, "Velocity output for other geometry types than cube and simplex");
238
239 // first we extract the corner indices for each scv for the CIV method
240 // for network grids there might be multiple intersection with the same geometryInInside
241 // we identify those by the indexInInside for now (assumes conforming grids at branching facets)
242 // here we keep track of them
243
1/2
✓ Branch 1 taken 845908 times.
✗ Branch 2 not taken.
 845908 std::vector<bool> handledScvf;
244 if (dim < dimWorld)
245 ✗ handledScvf.resize(element.subEntities(1), false);
246
247 // find the local face indices of the scvfs (for conforming meshes)
248
3/6
✓ Branch 1 taken 845908 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 845908 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 845908 times.
✗ Branch 8 not taken.
 3383632 std::vector<unsigned int> scvfIndexInInside(fvGeometry.numScvf());
249 845908 int localScvfIdx = 0;
250
11/18
✓ Branch 1 taken 5172 times.
✓ Branch 2 taken 840736 times.
✓ Branch 3 taken 3361344 times.
✓ Branch 4 taken 6762 times.
✓ Branch 5 taken 10 times.
✓ Branch 7 taken 5172 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 5172 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 25860 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 9412 times.
✓ Branch 15 taken 11276 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 20 taken 20688 times.
✗ Branch 21 not taken.
 5099708 for (const auto& intersection : intersections(gridGeometry_.gridView(), element))
251 {
252 if (dim < dimWorld)
253 ✗ if (handledScvf[intersection.indexInInside()])
254 continue;
255
256
5/8
✓ Branch 0 taken 3285270 times.
✓ Branch 1 taken 96762 times.
✓ Branch 2 taken 180384 times.
✓ Branch 3 taken 11276 times.
✓ Branch 4 taken 9412 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3488206 if (intersection.neighbor() || intersection.boundary())
257 {
258
2/2
✓ Branch 0 taken 3503600 times.
✓ Branch 1 taken 3382032 times.
 6885632 for (int i = 0; i < numScvfsPerFace; ++i)
259
1/3
✓ Branch 0 taken 31056 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 10510800 scvfIndexInInside[localScvfIdx++] = intersection.indexInInside();
260
261 // for surface and network grids mark that we handled this face
262 if (dim < dimWorld)
263 ✗ handledScvf[intersection.indexInInside()] = true;
264 }
265 }
266
267
4/15
✓ Branch 1 taken 845908 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 845908 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 845908 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 845908 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
 3378460 std::vector<Scalar> scvfFluxes(element.subEntities(1), 0.0);
268 845908 localScvfIdx = 0;
269
4/4
✓ Branch 0 taken 3503600 times.
✓ Branch 1 taken 845908 times.
✓ Branch 2 taken 3503600 times.
✓ Branch 3 taken 845908 times.
 8699016 for (auto&& scvf : scvfs(fvGeometry))
270 {
271
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6764064 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
272
2/2
✓ Branch 0 taken 3316328 times.
✓ Branch 1 taken 187272 times.
 3503600 if (!scvf.boundary())
273 {
274
1/2
✓ Branch 1 taken 3100798 times.
✗ Branch 2 not taken.
 3316328 FluxVariables fluxVars;
275
1/2
✓ Branch 1 taken 3316328 times.
✗ Branch 2 not taken.
 3316328 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
276
1/2
✓ Branch 1 taken 3316328 times.
✗ Branch 2 not taken.
 3316328 scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor();
277 }
278 else
279 {
280
2/3
✓ Branch 0 taken 8106 times.
✓ Branch 1 taken 16454 times.
✗ Branch 2 not taken.
 187272 auto bcTypes = problem_.boundaryTypes(element, scvf);
281
4/4
✓ Branch 0 taken 17503 times.
✓ Branch 1 taken 169769 times.
✓ Branch 2 taken 17503 times.
✓ Branch 3 taken 169769 times.
 374544 if (bcTypes.hasOnlyDirichlet())
282 {
283
1/2
✓ Branch 1 taken 13967 times.
✗ Branch 2 not taken.
 17503 FluxVariables fluxVars;
284
1/2
✓ Branch 1 taken 17503 times.
✗ Branch 2 not taken.
 17503 fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
285
1/2
✓ Branch 1 taken 17503 times.
✗ Branch 2 not taken.
 17503 scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor();
286 }
287 }
288
289 // increment scvf counter
290 3503600 localScvfIdx++;
291 }
292
293 // Correct boundary fluxes in case of Neumann conditions.
294 // In this general setting, it would be very difficult to
295 // calculate correct phase, i.e., volume, fluxes from arbitrary
296 // Neumann conditions. We approximate the Neumann flux by the
297 // flux on the opposite face. For extremely distorted grids this can
298 // lead to unexpected results (but then TPFA also leads to unexpected results).
299 845908 localScvfIdx = 0;
300
6/6
✓ Branch 0 taken 3503600 times.
✓ Branch 1 taken 845908 times.
✓ Branch 2 taken 3503600 times.
✓ Branch 3 taken 845908 times.
✓ Branch 4 taken 17600 times.
✓ Branch 5 taken 862400 times.
 8699016 for (auto&& scvf : scvfs(fvGeometry))
301 {
302
2/2
✓ Branch 0 taken 187272 times.
✓ Branch 1 taken 3316328 times.
 3503600 if (scvf.boundary())
303 {
304
2/3
✓ Branch 0 taken 10306 times.
✓ Branch 1 taken 14254 times.
✗ Branch 2 not taken.
 187272 auto bcTypes = problem_.boundaryTypes(element, scvf);
305
4/4
✓ Branch 0 taken 169769 times.
✓ Branch 1 taken 17503 times.
✓ Branch 2 taken 169769 times.
✓ Branch 3 taken 17503 times.
 374544 if (bcTypes.hasNeumann())
306 {
307 // for stationary velocity fields we can easily compute the correct velocity
308 // this special treatment makes sure that the velocity field is also correct on Neumann boundaries
309 // of tracer problems where the velocity field is given.
310 // (For Dirichlet boundaries no special treatment is necessary.)
311 if (stationaryVelocityField)
312 {
313 21836 const auto flux = AdvectionType::flux(problem_, element, fvGeometry, elemVolVars,
314 scvf, phaseIdx, elemFluxVarsCache);
315
316 39712 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
317 65508 scvfFluxes[scvfIndexInInside[localScvfIdx]] += flux / insideVolVars.extrusionFactor();
318 }
319 else
320 {
321 // check if we have Neumann no flow, we can just use 0
322
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 113605 times.
✗ Branch 2 not taken.
 147933 const auto neumannFlux = problem_.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
323 using NumEqVector = std::decay_t<decltype(neumannFlux)>;
324
2/2
✓ Branch 0 taken 140324 times.
✓ Branch 1 taken 7293 times.
 288970 if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30))
325 421920 scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
326
327 // otherwise, we try some reconstruction
328 // for cubes
329
1/2
✓ Branch 0 taken 2253 times.
✗ Branch 1 not taken.
 7283 else if (dim == 1 || geomType.isCube())
330 {
331
2/2
✓ Branch 0 taken 2223 times.
✓ Branch 1 taken 30 times.
 7293 const auto fIdx = scvfIndexInInside[localScvfIdx];
332
333 if constexpr (!modelIsCompositional)
334 {
335 // We assume that the density at the face equals the one at the cell center and reconstruct a volume flux from the Neumann mass flux.
336 10080 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
337 5040 const auto eqIdx = VolumeVariables::Indices::conti0EqIdx + phaseIdx;
338 16160 scvfFluxes[fIdx] += neumannFlux[eqIdx] / insideVolVars.density(phaseIdx) * scvf.area() * insideVolVars.extrusionFactor();
339 }
340 else
341 {
342 // For compositional models, we generally can't reconstruct the volume flow from the Neumann flux (which is a species flux rather
343 // than a phase flux here). Instead, we use the velocity of the opposing face.
344
2/2
✓ Branch 0 taken 2223 times.
✓ Branch 1 taken 30 times.
 2253 const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
345 6759 scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite];
346 }
347 }
348
349 // for simplices
350 ✗ else if (geomType.isSimplex())
351 ✗ scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
352
353 else
354 ✗ DUNE_THROW(Dune::NotImplemented, "Velocity computation at Neumann boundaries for cell-centered and prism/pyramid");
355 }
356 }
357 }
358
359 // increment scvf counter
360 3503600 localScvfIdx++;
361 }
362
363
364 845108 Velocity refVelocity;
365 // cubes: On the reference element simply average over opposite fluxes
366 // note that this is equal to a corner velocity interpolation method
367
1/2
✓ Branch 0 taken 5172 times.
✗ Branch 1 not taken.
 845108 if (dim == 1 || geomType.isCube())
368 {
369
2/2
✓ Branch 0 taken 1691016 times.
✓ Branch 1 taken 845908 times.
 2536924 for (int i = 0; i < dim; i++)
370 6763264 refVelocity[i] = 0.5 * (scvfFluxes[2*i + 1] - scvfFluxes[2*i]);
371 }
372 // simplices: Raviart-Thomas-0 interpolation evaluated at the cell center
373 ✗ else if (geomType.isSimplex())
374 {
375 ✗ for (int dimIdx = 0; dimIdx < dim; dimIdx++)
376 {
377 ✗ refVelocity[dimIdx] = -scvfFluxes[dim - 1 - dimIdx];
378 ✗ for (int fIdx = 0; fIdx < dim + 1; fIdx++)
379 ✗ refVelocity[dimIdx] += scvfFluxes[fIdx]/(dim + 1);
380 }
381 }
382 // 3D prism and pyramids
383 else
384 ✗ DUNE_THROW(Dune::NotImplemented, "Velocity computation for cell-centered and prism/pyramid");
385
386 845108 Velocity scvVelocity(0);
387 5972 jacobianT2.mtv(refVelocity, scvVelocity);
388
389
1/2
✓ Branch 1 taken 5172 times.
✗ Branch 2 not taken.
 846708 scvVelocity /= geometry.integrationElement(localPos);
390
391
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 5172 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5172 times.
✗ Branch 5 not taken.
 1691816 int eIdxGlobal = gridGeometry_.elementMapper().index(element);
392
393
2/4
✓ Branch 0 taken 845908 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 845908 times.
✗ Branch 3 not taken.
 1691816 velocity[eIdxGlobal] = scvVelocity;
394
395 } // cell-centered
396 1195120 }
397
398 private:
399 // The area of a subcontrolvolume face in a reference element.
400 // The 3d non-cube values have been calculated with quadrilateralArea3D
401 // of boxfvelementgeometry.hh.
402 1396848 static Scalar scvfReferenceArea_(Dune::GeometryType geomType, int fIdx)
403 {
404
1/2
✓ Branch 0 taken 1396848 times.
✗ Branch 1 not taken.
 1396848 if (dim == 1 || geomType.isCube())
405 {
406 ✗ return 1.0/(1 << (dim-1));
407 }
408 ✗ else if (geomType.isTriangle())
409 {
410 static const Scalar faceToArea[] = {0.372677996249965,
411 0.372677996249965,
412 0.235702260395516};
413 ✗ return faceToArea[fIdx];
414 }
415 ✗ else if (geomType.isTetrahedron())
416 {
417 static const Scalar faceToArea[] = {0.102062072615966,
418 0.102062072615966,
419 0.058925565098879,
420 0.102062072615966,
421 0.058925565098879,
422 0.058925565098879};
423 ✗ return faceToArea[fIdx];
424 }
425 ✗ else if (geomType.isPyramid())
426 {
427 static const Scalar faceToArea[] = {0.130437298687488,
428 0.130437298687488,
429 0.130437298687488,
430 0.130437298687488,
431 0.150923085635624,
432 0.1092906420717,
433 0.1092906420717,
434 0.0781735959970571};
435 ✗ return faceToArea[fIdx];
436 }
437 ✗ else if (geomType.isPrism())
438 {
439 static const Scalar faceToArea[] = {0.166666666666667,
440 0.166666666666667,
441 0.166666666666667,
442 0.186338998124982,
443 0.186338998124982,
444 0.117851130197758,
445 0.186338998124982,
446 0.186338998124982,
447 0.117851130197758};
448 ✗ return faceToArea[fIdx];
449 }
450 else {
451 ✗ DUNE_THROW(Dune::NotImplemented, "scvf area for unknown GeometryType");
452 }
453 }
454
455 private:
456 const Problem& problem_;
457 const GridGeometry& gridGeometry_;
458 const GridVariables& gridVariables_;
459
460 std::vector<int> cellNum_;
461 };
462
463 } // end namespace Dumux
464
465 #endif
466