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 | 97 | PorousMediumFlowVelocity(const GridVariables& gridVariables) | |
101 |
1/2✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
|
99 | : problem_(gridVariables.curGridVolVars().problem()) |
102 | , gridGeometry_(gridVariables.gridGeometry()) | ||
103 |
2/4✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
|
257 | , 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 22 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
68 | cellNum_.assign(gridGeometry_.gridView().size(dim), 0); |
110 | |||
111 |
7/19✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 34 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 21950 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.
|
21996 | for (const auto& element : elements(gridGeometry_.gridView())) |
112 |
5/5✓ Branch 0 taken 76200 times.
✓ Branch 1 taken 33440 times.
✓ Branch 2 taken 76200 times.
✓ Branch 3 taken 30562 times.
✓ Branch 4 taken 2878 times.
|
204890 | for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx) |
113 |
2/4✓ Branch 1 taken 87712 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87712 times.
✗ Branch 5 not taken.
|
175424 | ++cellNum_[gridGeometry_.vertexMapper().subIndex(element, vIdx, dim)]; |
114 | } | ||
115 | 34 | } | |
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 | 1215480 | 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 | 1215480 | const auto geometry = element.geometry(); | |
129 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
1215480 | 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.
|
15987460 | auto upwindTerm = [phaseIdx](const auto& volVars) { return volVars.mobility(phaseIdx); }; |
133 | |||
134 | if constexpr (isBox && dim == 1) | ||
135 | { | ||
136 | 2400 | Velocity tmpVelocity(0.0); | |
137 | 7200 | tmpVelocity = (geometry.corner(1) - geometry.corner(0)); | |
138 | 4800 | tmpVelocity /= tmpVelocity.two_norm(); | |
139 | |||
140 |
4/4✓ Branch 0 taken 2422 times.
✓ Branch 1 taken 2400 times.
✓ Branch 2 taken 2422 times.
✓ Branch 3 taken 2400 times.
|
7222 | for (auto&& scvf : scvfs(fvGeometry)) |
141 | { | ||
142 |
2/2✓ Branch 0 taken 22 times.
✓ Branch 1 taken 2400 times.
|
2422 | if (scvf.boundary()) |
143 | 22 | continue; | |
144 | |||
145 | // instantiate the flux variables | ||
146 | 2400 | FluxVariables fluxVars; | |
147 | 2400 | 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 | 2400 | Scalar localArea = scvfReferenceArea_(geomType, scvf.index()); | |
152 | 2400 | Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea; | |
153 | 4800 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
154 | 2400 | flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area(); | |
155 | 2400 | tmpVelocity *= flux; | |
156 | |||
157 | 4800 | const int eIdxGlobal = gridGeometry_.elementMapper().index(element); | |
158 | 4800 | 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.
|
1213080 | const auto refElement = referenceElement(geometry); |
166 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 800 times.
|
1213080 | const auto& localPos = refElement.position(0, 0); |
167 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 800 times.
|
1213080 | 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 356772 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.
|
1070316 | ScvVelocities scvVelocities(fvGeometry.numScv()); |
173 | 356772 | scvVelocities = 0; | |
174 | |||
175 |
4/4✓ Branch 0 taken 1761964 times.
✓ Branch 1 taken 356772 times.
✓ Branch 2 taken 1761964 times.
✓ Branch 3 taken 356772 times.
|
2475508 | for (auto&& scvf : scvfs(fvGeometry)) |
176 | { | ||
177 |
2/2✓ Branch 0 taken 334876 times.
✓ Branch 1 taken 1427088 times.
|
1761964 | if (scvf.boundary()) |
178 | 334876 | 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.
|
2854176 | 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.
|
1427088 | const auto jacobianT1 = geometry.jacobianTransposed(localPosIP); |
185 | 1427088 | const auto globalNormal = scvf.unitOuterNormal(); | |
186 | 1427088 | GlobalPosition localNormal(0); | |
187 | jacobianT1.mv(globalNormal, localNormal); | ||
188 | 1427088 | localNormal /= localNormal.two_norm(); | |
189 | |||
190 | // instantiate the flux variables | ||
191 |
1/2✓ Branch 1 taken 1085248 times.
✗ Branch 2 not taken.
|
1427088 | FluxVariables fluxVars; |
192 |
1/2✓ Branch 1 taken 1427088 times.
✗ Branch 2 not taken.
|
1427088 | 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 1427088 times.
✗ Branch 2 not taken.
|
1427088 | Scalar localArea = scvfReferenceArea_(geomType, scvf.index()); |
197 |
1/2✓ Branch 1 taken 1427088 times.
✗ Branch 2 not taken.
|
1427088 | Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea; |
198 | 2854176 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
199 | 1427088 | flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area(); | |
200 | |||
201 | // transform the volume flux into a velocity vector | ||
202 | 1427088 | Velocity tmpVelocity = localNormal; | |
203 | 1427088 | tmpVelocity *= flux; | |
204 | |||
205 | 4281264 | scvVelocities[scvf.insideScvIdx()] += tmpVelocity; | |
206 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1427088 times.
|
2854176 | scvVelocities[scvf.outsideScvIdx()] += tmpVelocity; |
207 | } | ||
208 | |||
209 | // transform vertex velocities from local to global coordinates | ||
210 |
4/4✓ Branch 0 taken 1427088 times.
✓ Branch 1 taken 356772 times.
✓ Branch 2 taken 1427088 times.
✓ Branch 3 taken 356772 times.
|
3567720 | for (auto&& scv : scvs(fvGeometry)) |
211 | { | ||
212 | 1427088 | int vIdxGlobal = scv.dofIndex(); | |
213 | |||
214 | // calculate the subcontrolvolume velocity by the Piola transformation | ||
215 | 1427088 | Velocity scvVelocity(0); | |
216 | |||
217 | 2854176 | jacobianT2.mtv(scvVelocities[scv.indexInElement()], scvVelocity); | |
218 |
1/2✓ Branch 1 taken 263888 times.
✗ Branch 2 not taken.
|
2590288 | scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal]; |
219 | // add up the wetting phase subcontrolvolume velocities for each vertex | ||
220 | 4281264 | 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.
|
858900 | const int numScvfsPerFace = isMpfa ? element.template subEntity<1>(0).geometry().corners() : 1; |
232 | |||
233 |
3/8✗ Branch 0 not taken.
✓ Branch 1 taken 851136 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 856308 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 851136 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2563752 | 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 851136 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
856308 | 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 856308 times.
✗ Branch 2 not taken.
|
856308 | 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 856308 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 856308 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 856308 times.
✗ Branch 8 not taken.
|
3425232 | std::vector<unsigned int> scvfIndexInInside(fvGeometry.numScvf()); |
249 | 856308 | int localScvfIdx = 0; | |
250 |
11/18✓ Branch 1 taken 5172 times.
✓ Branch 2 taken 851136 times.
✓ Branch 3 taken 3402944 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.
|
5162108 | 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 3325694 times.
✓ Branch 1 taken 97938 times.
✓ Branch 2 taken 182736 times.
✓ Branch 3 taken 11276 times.
✓ Branch 4 taken 9412 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
3530982 | if (intersection.neighbor() || intersection.boundary()) |
257 | { | ||
258 |
2/2✓ Branch 0 taken 3545200 times.
✓ Branch 1 taken 3423632 times.
|
6968832 | for (int i = 0; i < numScvfsPerFace; ++i) |
259 |
1/3✓ Branch 0 taken 31056 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
10635600 | 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 856308 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 856308 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 856308 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 856308 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.
|
3420060 | std::vector<Scalar> scvfFluxes(element.subEntities(1), 0.0); |
268 | 856308 | localScvfIdx = 0; | |
269 |
4/4✓ Branch 0 taken 3545200 times.
✓ Branch 1 taken 856308 times.
✓ Branch 2 taken 3545200 times.
✓ Branch 3 taken 856308 times.
|
8803016 | 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.
|
6847264 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; |
272 |
2/2✓ Branch 0 taken 3355576 times.
✓ Branch 1 taken 189624 times.
|
3545200 | if (!scvf.boundary()) |
273 | { | ||
274 |
1/2✓ Branch 1 taken 3141526 times.
✗ Branch 2 not taken.
|
3355576 | FluxVariables fluxVars; |
275 |
1/2✓ Branch 1 taken 3355576 times.
✗ Branch 2 not taken.
|
3355576 | fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); |
276 |
1/2✓ Branch 1 taken 3355576 times.
✗ Branch 2 not taken.
|
3355576 | 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.
|
189624 | auto bcTypes = problem_.boundaryTypes(element, scvf); |
281 |
4/4✓ Branch 0 taken 17555 times.
✓ Branch 1 taken 172069 times.
✓ Branch 2 taken 17555 times.
✓ Branch 3 taken 172069 times.
|
379248 | if (bcTypes.hasOnlyDirichlet()) |
282 | { | ||
283 |
1/2✓ Branch 1 taken 14039 times.
✗ Branch 2 not taken.
|
17555 | FluxVariables fluxVars; |
284 |
1/2✓ Branch 1 taken 17555 times.
✗ Branch 2 not taken.
|
17555 | fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); |
285 |
1/2✓ Branch 1 taken 17555 times.
✗ Branch 2 not taken.
|
17555 | scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor(); |
286 | } | ||
287 | } | ||
288 | |||
289 | // increment scvf counter | ||
290 | 3545200 | 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 | 856308 | localScvfIdx = 0; | |
300 |
6/6✓ Branch 0 taken 3545200 times.
✓ Branch 1 taken 856308 times.
✓ Branch 2 taken 3545200 times.
✓ Branch 3 taken 856308 times.
✓ Branch 4 taken 17600 times.
✓ Branch 5 taken 862400 times.
|
8803016 | for (auto&& scvf : scvfs(fvGeometry)) |
301 | { | ||
302 |
2/2✓ Branch 0 taken 189624 times.
✓ Branch 1 taken 3355576 times.
|
3545200 | if (scvf.boundary()) |
303 | { | ||
304 |
2/3✓ Branch 0 taken 10306 times.
✓ Branch 1 taken 14254 times.
✗ Branch 2 not taken.
|
189624 | auto bcTypes = problem_.boundaryTypes(element, scvf); |
305 |
4/4✓ Branch 0 taken 172069 times.
✓ Branch 1 taken 17555 times.
✓ Branch 2 taken 172069 times.
✓ Branch 3 taken 17555 times.
|
379248 | 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 113505 times.
✗ Branch 2 not taken.
|
150233 | const auto neumannFlux = problem_.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); |
323 | using NumEqVector = std::decay_t<decltype(neumannFlux)>; | ||
324 |
2/2✓ Branch 0 taken 142644 times.
✓ Branch 1 taken 7273 times.
|
293570 | if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30)) |
325 | 428880 | scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; | |
326 | |||
327 | // otherwise, we try some reconstruction | ||
328 | // for cubes | ||
329 |
1/2✓ Branch 0 taken 2233 times.
✗ Branch 1 not taken.
|
7263 | else if (dim == 1 || geomType.isCube()) |
330 | { | ||
331 |
2/2✓ Branch 0 taken 2203 times.
✓ Branch 1 taken 30 times.
|
7273 | 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 2203 times.
✓ Branch 1 taken 30 times.
|
2233 | const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1; |
345 | 6699 | 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 | 3545200 | localScvfIdx++; | |
361 | } | ||
362 | |||
363 | |||
364 | 855508 | 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.
|
855508 | if (dim == 1 || geomType.isCube()) |
368 | { | ||
369 |
2/2✓ Branch 0 taken 1711816 times.
✓ Branch 1 taken 856308 times.
|
2568124 | for (int i = 0; i < dim; i++) |
370 | 6846464 | 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 | 855508 | Velocity scvVelocity(0); | |
387 | 5972 | jacobianT2.mtv(refVelocity, scvVelocity); | |
388 | |||
389 |
1/2✓ Branch 1 taken 5172 times.
✗ Branch 2 not taken.
|
857108 | 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.
|
1712616 | int eIdxGlobal = gridGeometry_.elementMapper().index(element); |
392 | |||
393 |
2/4✓ Branch 0 taken 856308 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 856308 times.
✗ Branch 3 not taken.
|
1712616 | velocity[eIdxGlobal] = scvVelocity; |
394 | |||
395 | } // cell-centered | ||
396 | 1213080 | } | |
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 | 1427088 | static Scalar scvfReferenceArea_(Dune::GeometryType geomType, int fIdx) | |
403 | { | ||
404 |
1/2✓ Branch 0 taken 1427088 times.
✗ Branch 1 not taken.
|
1427088 | 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 |