| 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-FileCopyrightText: 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/common/reservedvector.hh> | ||
| 22 | #include <dune/geometry/type.hh> | ||
| 23 | #include <dune/geometry/referenceelements.hh> | ||
| 24 | |||
| 25 | #include <dumux/common/parameters.hh> | ||
| 26 | #include <dumux/discretization/extrusion.hh> | ||
| 27 | #include <dumux/discretization/method.hh> | ||
| 28 | #include <dumux/discretization/elementsolution.hh> | ||
| 29 | #include <dumux/flux/traits.hh> | ||
| 30 | |||
| 31 | namespace Dumux { | ||
| 32 | |||
| 33 | #ifndef DOXYGEN | ||
| 34 | namespace Detail { | ||
| 35 | // helper structs and functions detecting if the model is compositional | ||
| 36 | |||
| 37 | template <typename T, typename ...Ts> | ||
| 38 | using MoleFractionDetector = decltype(std::declval<T>().moleFraction(std::declval<Ts>()...)); | ||
| 39 | |||
| 40 | template<class T, typename ...Args> | ||
| 41 | static constexpr bool hasMoleFraction() | ||
| 42 | { return Dune::Std::is_detected<MoleFractionDetector, T, Args...>::value; } | ||
| 43 | |||
| 44 | template <typename T, typename ...Ts> | ||
| 45 | using MassFractionDetector = decltype(std::declval<T>().massFraction(std::declval<Ts>()...)); | ||
| 46 | |||
| 47 | template<class T, typename ...Args> | ||
| 48 | static constexpr bool hasMassFraction() | ||
| 49 | { return Dune::Std::is_detected<MassFractionDetector, T, Args...>::value; } | ||
| 50 | |||
| 51 | } // end namespace Detail | ||
| 52 | #endif // DOXYGEN | ||
| 53 | |||
| 54 | /*! | ||
| 55 | * \ingroup PorousmediumflowModels | ||
| 56 | * \brief Velocity computation for implicit (porous media) models. | ||
| 57 | */ | ||
| 58 | template<class GridVariables, class FluxVariables> | ||
| 59 | 99 | class PorousMediumFlowVelocity | |
| 60 | { | ||
| 61 | using GridGeometry = typename GridVariables::GridGeometry; | ||
| 62 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 63 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 64 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 65 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 66 | using GridView = typename GridGeometry::GridView; | ||
| 67 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 68 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 69 | using ElementFluxVarsCache = typename GridVariables::GridFluxVariablesCache::LocalView; | ||
| 70 | using VolumeVariables = typename GridVariables::VolumeVariables; | ||
| 71 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 72 | using FluidSystem = typename VolumeVariables::FluidSystem; | ||
| 73 | using Scalar = typename GridVariables::Scalar; | ||
| 74 | using FluxTraits = typename Dumux::FluxTraits<FluxVariables>; | ||
| 75 | using AdvectionType = typename FluxVariables::AdvectionType; | ||
| 76 | |||
| 77 | static constexpr int dim = GridView::dimension; | ||
| 78 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 79 | static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box; | ||
| 80 | static constexpr bool isDiamond = GridGeometry::discMethod == DiscretizationMethods::fcdiamond; | ||
| 81 | static constexpr bool stationaryVelocityField = FluxTraits::hasStationaryVelocityField(); | ||
| 82 | |||
| 83 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 84 | |||
| 85 | using Problem = typename GridVolumeVariables::Problem; | ||
| 86 | |||
| 87 | static constexpr bool modelIsCompositional = Detail::hasMoleFraction<typename GridVolumeVariables::VolumeVariables, int, int>() || | ||
| 88 | Detail::hasMassFraction<typename GridVolumeVariables::VolumeVariables, int, int>(); | ||
| 89 | |||
| 90 | static_assert(VolumeVariables::numFluidPhases() >= 1, "Velocity output only makes sense for models with fluid phases."); | ||
| 91 | |||
| 92 | using Velocity = Dune::FieldVector<Scalar, dimWorld>; | ||
| 93 | using ReferenceElementVelocity = Dune::FieldVector<Scalar, dim>; | ||
| 94 | public: | ||
| 95 | static constexpr int numFluidPhases = VolumeVariables::numFluidPhases(); | ||
| 96 | using VelocityVector = std::vector<Velocity>; | ||
| 97 | |||
| 98 | /*! | ||
| 99 | * \brief Constructor initializes the static data with the initial solution. | ||
| 100 | * | ||
| 101 | * \param gridVariables The grid variables | ||
| 102 | */ | ||
| 103 |
1/2✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
|
99 | PorousMediumFlowVelocity(const GridVariables& gridVariables) |
| 104 |
1/2✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
|
99 | : problem_(gridVariables.curGridVolVars().problem()) |
| 105 |
1/2✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
|
99 | , gridGeometry_(gridVariables.gridGeometry()) |
| 106 |
1/2✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
|
99 | , gridVariables_(gridVariables) |
| 107 | { | ||
| 108 | // set the number of scvs the vertices are connected to | ||
| 109 | if constexpr (isBox && dim > 1) | ||
| 110 | { | ||
| 111 | // resize to the number of vertices of the grid | ||
| 112 |
2/5✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
|
35 | cellNum_.assign(gridGeometry_.gridView().size(dim), 0); |
| 113 | |||
| 114 |
4/12✓ Branch 1 taken 35 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19850 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 6 taken 3278 times.
✗ Branch 3 not taken.
|
39691 | for (const auto& element : elements(gridGeometry_.gridView())) |
| 115 |
4/5✓ Branch 1 taken 32940 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13112 times.
✓ Branch 4 taken 3278 times.
✓ Branch 0 taken 66200 times.
|
99140 | for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx) |
| 116 |
1/2✓ Branch 1 taken 79312 times.
✗ Branch 2 not taken.
|
79312 | ++cellNum_[gridGeometry_.vertexMapper().subIndex(element, vIdx, dim)]; |
| 117 | } | ||
| 118 | 35 | } | |
| 119 | |||
| 120 | //! Calculates the velocities for the scvs in the element. | ||
| 121 | //! We assume the local containers to be bound to the complete stencil. | ||
| 122 | 1227640 | void calculateVelocity(VelocityVector& velocity, | |
| 123 | const Element& element, | ||
| 124 | const FVElementGeometry& fvGeometry, | ||
| 125 | const ElementVolumeVariables& elemVolVars, | ||
| 126 | const ElementFluxVarsCache& elemFluxVarsCache, | ||
| 127 | int phaseIdx) const | ||
| 128 | { | ||
| 129 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
1227640 | const auto geometry = element.geometry(); |
| 130 | 1227640 | const Dune::GeometryType geomType = geometry.type(); | |
| 131 | |||
| 132 | // the upwind term to be used for the volume flux evaluation | ||
| 133 |
4/10✓ Branch 0 taken 174854 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 174854 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 131746 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 131746 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
6687059 | auto upwindTerm = [phaseIdx](const auto& volVars) { return volVars.mobility(phaseIdx); }; |
| 134 | |||
| 135 | if constexpr (isBox && dim == 1) | ||
| 136 | { | ||
| 137 | 2960 | Velocity tmpVelocity(0.0); | |
| 138 | 2960 | tmpVelocity = (geometry.corner(1) - geometry.corner(0)); | |
| 139 | 2960 | tmpVelocity /= tmpVelocity.two_norm(); | |
| 140 | |||
| 141 |
4/4✓ Branch 0 taken 36 times.
✓ Branch 1 taken 2960 times.
✓ Branch 2 taken 2996 times.
✓ Branch 3 taken 2960 times.
|
5956 | for (auto&& scvf : scvfs(fvGeometry)) |
| 142 | { | ||
| 143 |
2/2✓ Branch 0 taken 36 times.
✓ Branch 1 taken 2960 times.
|
2996 | if (scvf.boundary()) |
| 144 | 36 | continue; | |
| 145 | |||
| 146 | // instantiate the flux variables | ||
| 147 | 2960 | FluxVariables fluxVars; | |
| 148 | 2960 | fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); | |
| 149 | |||
| 150 | // get the volume flux divided by the area of the | ||
| 151 | // subcontrolvolume face in the reference element | ||
| 152 | 2960 | Scalar localArea = scvfReferenceArea_(geomType, scvf.index()); | |
| 153 | 2960 | Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea; | |
| 154 | 2960 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 155 | 2960 | flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area(); | |
| 156 | 2960 | tmpVelocity *= flux; | |
| 157 | |||
| 158 | 2960 | const int eIdxGlobal = gridGeometry_.elementMapper().index(element); | |
| 159 | 2960 | velocity[eIdxGlobal] = tmpVelocity; | |
| 160 | } | ||
| 161 | return; | ||
| 162 | } | ||
| 163 | |||
| 164 | // get the transposed Jacobian of the element mapping | ||
| 165 | using Dune::referenceElement; | ||
| 166 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 800 times.
|
1224680 | const auto localPos = referenceElement(geometry).position(0, 0); |
| 167 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
1224680 | const auto jacobianT2 = geometry.jacobianTransposed(localPos); |
| 168 | |||
| 169 | if constexpr (isBox) | ||
| 170 | { | ||
| 171 | using ScvVelocities = Dune::BlockVector<Velocity>; | ||
| 172 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
365772 | ScvVelocities scvVelocities(fvGeometry.numScv()); |
| 173 | 731544 | scvVelocities = 0; | |
| 174 | |||
| 175 |
4/4✓ Branch 0 taken 338644 times.
✓ Branch 1 taken 1463088 times.
✓ Branch 2 taken 1801732 times.
✓ Branch 3 taken 365772 times.
|
2167504 | for (auto&& scvf : scvfs(fvGeometry)) |
| 176 | { | ||
| 177 |
2/2✓ Branch 0 taken 338644 times.
✓ Branch 1 taken 1463088 times.
|
1801732 | if (scvf.boundary()) |
| 178 | 338644 | continue; | |
| 179 | |||
| 180 | // local position of integration point | ||
| 181 |
2/4✓ Branch 1 taken 311888 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 311888 times.
✗ Branch 5 not taken.
|
1463088 | const auto localPosIP = geometry.local(scvf.ipGlobal()); |
| 182 | |||
| 183 | // Transformation of the global normal vector to normal vector in the reference element | ||
| 184 | 1463088 | const auto jacobianT1 = geometry.jacobianTransposed(localPosIP); | |
| 185 | 1463088 | const auto globalNormal = scvf.unitOuterNormal(); | |
| 186 | 1463088 | GlobalPosition localNormal(0); | |
| 187 |
2/2✓ Branch 0 taken 2926176 times.
✓ Branch 1 taken 1463088 times.
|
4389264 | jacobianT1.mv(globalNormal, localNormal); |
| 188 | 2926176 | localNormal /= localNormal.two_norm(); | |
| 189 | |||
| 190 | // instantiate the flux variables | ||
| 191 | 1463088 | FluxVariables fluxVars; | |
| 192 |
1/2✓ Branch 1 taken 1463088 times.
✗ Branch 2 not taken.
|
1463088 | 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 1463088 times.
✗ Branch 2 not taken.
|
1463088 | Scalar localArea = scvfReferenceArea_(geomType, scvf.index()); |
| 197 |
1/2✓ Branch 1 taken 1463088 times.
✗ Branch 2 not taken.
|
1463088 | Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea; |
| 198 | 1463088 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 199 | 1463088 | flux /= insideVolVars.extrusionFactor() * Extrusion::area(fvGeometry, scvf) / scvf.area(); | |
| 200 | |||
| 201 | // transform the volume flux into a velocity vector | ||
| 202 | 1463088 | Velocity tmpVelocity = localNormal; | |
| 203 | 1463088 | tmpVelocity *= flux; | |
| 204 | |||
| 205 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1463088 times.
|
2926176 | scvVelocities[scvf.insideScvIdx()] += tmpVelocity; |
| 206 | 2926176 | scvVelocities[scvf.outsideScvIdx()] += tmpVelocity; | |
| 207 | } | ||
| 208 | |||
| 209 | // transform vertex velocities from local to global coordinates | ||
| 210 |
2/2✓ Branch 0 taken 1463088 times.
✓ Branch 1 taken 365772 times.
|
1828860 | for (auto&& scv : scvs(fvGeometry)) |
| 211 | { | ||
| 212 | 1463088 | int vIdxGlobal = scv.dofIndex(); | |
| 213 | |||
| 214 | // calculate the subcontrolvolume velocity by the Piola transformation | ||
| 215 | 1463088 | Velocity scvVelocity(0); | |
| 216 | |||
| 217 |
2/3✓ Branch 0 taken 2302400 times.
✓ Branch 1 taken 1463088 times.
✗ Branch 2 not taken.
|
5228576 | jacobianT2.mtv(scvVelocities[scv.indexInElement()], scvVelocity); |
| 218 | 1463088 | scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal]; | |
| 219 | // add up the wetting phase subcontrolvolume velocities for each vertex | ||
| 220 | 2926176 | velocity[vIdxGlobal] += scvVelocity; | |
| 221 | } | ||
| 222 | 365772 | } | |
| 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 |
4/8✓ Branch 1 taken 828528 times.
✓ Branch 2 taken 27800 times.
✓ Branch 4 taken 2400 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2400 times.
✗ Branch 7 not taken.
✗ Branch 0 not taken.
✗ Branch 3 not taken.
|
866300 | const int numScvfsPerFace = isMpfa ? element.template subEntity<1>(0).geometry().corners() : 1; |
| 232 | |||
| 233 |
2/5✓ Branch 1 taken 853736 times.
✓ Branch 2 taken 5172 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 0 not taken.
|
858908 | if (fvGeometry.numScvf() != element.subEntities(1)*numScvfsPerFace) |
| 234 |
1/22✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 33 taken 5172 times.
✗ Branch 34 not taken.
|
5172 | DUNE_THROW(Dune::NotImplemented, "Velocity output for non-conforming grids"); |
| 235 | |||
| 236 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
853736 | if (!geomType.isCube() && !geomType.isSimplex()) |
| 237 |
0/20✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
|
853736 | 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 858908 times.
✗ Branch 2 not taken.
|
858908 | 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 |
1/2✓ Branch 1 taken 858908 times.
✗ Branch 2 not taken.
|
858908 | std::vector<unsigned int> scvfIndexInInside(fvGeometry.numScvf()); |
| 249 | 858908 | int localScvfIdx = 0; | |
| 250 |
9/14✓ Branch 1 taken 858908 times.
✓ Branch 2 taken 3413344 times.
✓ Branch 7 taken 5172 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 25860 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 9412 times.
✓ Branch 12 taken 11276 times.
✓ Branch 4 taken 3330092 times.
✓ Branch 5 taken 86834 times.
✗ Branch 6 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 3 taken 1590 times.
|
4329144 | for (const auto& intersection : intersections(gridGeometry_.gridView(), element)) |
| 251 | { | ||
| 252 | if (dim < dimWorld) | ||
| 253 | ✗ | if (handledScvf[intersection.indexInInside()]) | |
| 254 | ✗ | continue; | |
| 255 | |||
| 256 |
4/5✓ Branch 0 taken 173688 times.
✓ Branch 1 taken 20688 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9412 times.
✓ Branch 4 taken 11276 times.
|
215064 | if (intersection.neighbor() || intersection.boundary()) |
| 257 | { | ||
| 258 |
2/2✓ Branch 0 taken 3555600 times.
✓ Branch 1 taken 3434032 times.
|
6989632 | for (int i = 0; i < numScvfsPerFace; ++i) |
| 259 | 3555600 | 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 |
1/2✓ Branch 1 taken 5172 times.
✗ Branch 2 not taken.
|
858908 | Dune::ReservedVector<Scalar, 2*dim> scvfFluxes(element.subEntities(1), 0.0); |
| 268 | 858908 | localScvfIdx = 0; | |
| 269 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 3555600 times.
✓ Branch 2 taken 858908 times.
✗ Branch 3 not taken.
|
4414508 | for (auto&& scvf : scvfs(fvGeometry)) |
| 270 | { | ||
| 271 |
2/3✓ Branch 1 taken 3365612 times.
✓ Branch 2 taken 189988 times.
✗ Branch 0 not taken.
|
3555600 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; |
| 272 |
2/2✓ Branch 0 taken 3365612 times.
✓ Branch 1 taken 189988 times.
|
3555600 | if (!scvf.boundary()) |
| 273 | { | ||
| 274 | 3365612 | FluxVariables fluxVars; | |
| 275 |
1/2✓ Branch 1 taken 3365612 times.
✗ Branch 2 not taken.
|
3365612 | fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); |
| 276 |
1/2✓ Branch 1 taken 3365612 times.
✗ Branch 2 not taken.
|
3365612 | scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor(); |
| 277 | } | ||
| 278 | else | ||
| 279 | { | ||
| 280 |
4/4✓ Branch 0 taken 11513 times.
✓ Branch 1 taken 178475 times.
✓ Branch 2 taken 3644 times.
✓ Branch 3 taken 1436 times.
|
191708 | auto bcTypes = problem_.boundaryTypes(element, scvf); |
| 281 |
2/2✓ Branch 0 taken 17519 times.
✓ Branch 1 taken 172469 times.
|
189988 | if (bcTypes.hasOnlyDirichlet()) |
| 282 | { | ||
| 283 | 17519 | FluxVariables fluxVars; | |
| 284 |
1/2✓ Branch 1 taken 17519 times.
✗ Branch 2 not taken.
|
17519 | fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); |
| 285 |
1/2✓ Branch 1 taken 17519 times.
✗ Branch 2 not taken.
|
17519 | scvfFluxes[scvfIndexInInside[localScvfIdx]] += fluxVars.advectiveFlux(phaseIdx, upwindTerm)/insideVolVars.extrusionFactor(); |
| 286 | } | ||
| 287 | } | ||
| 288 | |||
| 289 | // increment scvf counter | ||
| 290 | 3555600 | 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 | 858908 | localScvfIdx = 0; | |
| 300 |
4/4✓ Branch 0 taken 189988 times.
✓ Branch 1 taken 3365612 times.
✓ Branch 2 taken 3555600 times.
✓ Branch 3 taken 858908 times.
|
4414508 | for (auto&& scvf : scvfs(fvGeometry)) |
| 301 | { | ||
| 302 |
2/2✓ Branch 0 taken 189988 times.
✓ Branch 1 taken 3365612 times.
|
3555600 | if (scvf.boundary()) |
| 303 | { | ||
| 304 |
4/4✓ Branch 0 taken 165831 times.
✓ Branch 1 taken 24157 times.
✓ Branch 2 taken 1436 times.
✓ Branch 3 taken 3644 times.
|
191708 | auto bcTypes = problem_.boundaryTypes(element, scvf); |
| 305 |
2/2✓ Branch 0 taken 172469 times.
✓ Branch 1 taken 17519 times.
|
189988 | 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 | 21836 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 317 | 21836 | 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 115105 times.
✗ Branch 2 not taken.
|
150633 | const auto neumannFlux = problem_.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); |
| 323 | using NumEqVector = std::decay_t<decltype(neumannFlux)>; | ||
| 324 |
4/4✓ Branch 0 taken 509469 times.
✓ Branch 1 taken 143244 times.
✓ Branch 2 taken 138864 times.
✓ Branch 3 taken 3273 times.
|
798439 | if (Dune::FloatCmp::eq<NumEqVector, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, NumEqVector(0.0), 1e-30)) |
| 325 | 142980 | scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; | |
| 326 | |||
| 327 | // otherwise, we try some reconstruction | ||
| 328 | // for cubes | ||
| 329 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
7553 | else if (dim == 1 || geomType.isCube()) |
| 330 | { | ||
| 331 |
2/2✓ Branch 0 taken 2203 times.
✓ Branch 1 taken 30 times.
|
7653 | 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 | 5420 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 337 | 5420 | const auto eqIdx = VolumeVariables::Indices::conti0EqIdx + phaseIdx; | |
| 338 | 5420 | scvfFluxes[fIdx] += neumannFlux[eqIdx] / insideVolVars.density(phaseIdx) * scvf.area(); | |
| 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 | 2233 | scvfFluxes[fIdx] = -scvfFluxes[fIdxOpposite]; | |
| 346 | } | ||
| 347 | } | ||
| 348 | |||
| 349 | // for simplices | ||
| 350 | ✗ | else if (geomType.isSimplex()) | |
| 351 | ✗ | scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; | |
| 352 | |||
| 353 | else | ||
| 354 |
0/20✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
|
150633 | DUNE_THROW(Dune::NotImplemented, "Velocity computation at Neumann boundaries for cell-centered and prism/pyramid"); |
| 355 | } | ||
| 356 | } | ||
| 357 | } | ||
| 358 | |||
| 359 | // increment scvf counter | ||
| 360 | 3555600 | localScvfIdx++; | |
| 361 | } | ||
| 362 | |||
| 363 | |||
| 364 |
2/2✓ Branch 0 taken 5172 times.
✓ Branch 1 taken 800 times.
|
858908 | ReferenceElementVelocity 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/4✗ Branch 0 not taken.
✓ Branch 1 taken 858108 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
858108 | if (dim == 1 || geomType.isCube()) |
| 368 | { | ||
| 369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 858908 times.
|
858908 | assert(scvfFluxes.size() == 2*dim); |
| 370 |
2/2✓ Branch 0 taken 1717016 times.
✓ Branch 1 taken 858908 times.
|
2575924 | for (int i = 0; i < dim; ++i) |
| 371 | 1717016 | refVelocity[i] = 0.5 * (scvfFluxes[2*i + 1] - scvfFluxes[2*i]); | |
| 372 | } | ||
| 373 | |||
| 374 | // simplices: Raviart-Thomas-0 interpolation evaluated at the cell center | ||
| 375 | ✗ | else if (geomType.isSimplex()) | |
| 376 | { | ||
| 377 | ✗ | assert(scvfFluxes.size() == dim+1); | |
| 378 | ✗ | for (int i = 0; i < dim; ++i) | |
| 379 | { | ||
| 380 | ✗ | refVelocity[i] = -scvfFluxes[dim - 1 - i]; | |
| 381 | ✗ | for (int fIdx = 0; fIdx < dim + 1; ++fIdx) | |
| 382 | ✗ | refVelocity[i] += scvfFluxes[fIdx]/(dim + 1); | |
| 383 | } | ||
| 384 | } | ||
| 385 | |||
| 386 | // 3D prism and pyramids | ||
| 387 | else | ||
| 388 |
0/20✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
|
852936 | DUNE_THROW(Dune::NotImplemented, "Velocity computation for cell-centered and prism/pyramid"); |
| 389 | |||
| 390 | 858108 | Velocity scvVelocity(0); | |
| 391 |
2/3✓ Branch 0 taken 1705872 times.
✓ Branch 1 taken 858108 times.
✗ Branch 2 not taken.
|
2564780 | jacobianT2.mtv(refVelocity, scvVelocity); |
| 392 | |||
| 393 |
2/2✓ Branch 0 taken 1716216 times.
✓ Branch 1 taken 858108 times.
|
2575124 | scvVelocity /= geometry.integrationElement(localPos); |
| 394 | |||
| 395 |
1/2✓ Branch 1 taken 5172 times.
✗ Branch 2 not taken.
|
858908 | int eIdxGlobal = gridGeometry_.elementMapper().index(element); |
| 396 | |||
| 397 | 858908 | velocity[eIdxGlobal] = scvVelocity; | |
| 398 | |||
| 399 | 858908 | } // cell-centered | |
| 400 | 1224680 | } | |
| 401 | |||
| 402 | private: | ||
| 403 | // The area of a subcontrolvolume face in a reference element. | ||
| 404 | // The 3d non-cube values have been calculated with quadrilateralArea3D | ||
| 405 | // of boxfvelementgeometry.hh. | ||
| 406 |
1/2✓ Branch 0 taken 1463088 times.
✗ Branch 1 not taken.
|
1463088 | static Scalar scvfReferenceArea_(Dune::GeometryType geomType, int fIdx) |
| 407 | { | ||
| 408 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
1463088 | if (dim == 1 || geomType.isCube()) |
| 409 | { | ||
| 410 | return 1.0/(1 << (dim-1)); | ||
| 411 | } | ||
| 412 | ✗ | else if (geomType.isTriangle()) | |
| 413 | { | ||
| 414 | static const Scalar faceToArea[] = {0.372677996249965, | ||
| 415 | 0.372677996249965, | ||
| 416 | 0.235702260395516}; | ||
| 417 | ✗ | return faceToArea[fIdx]; | |
| 418 | } | ||
| 419 | ✗ | else if (geomType.isTetrahedron()) | |
| 420 | { | ||
| 421 | static const Scalar faceToArea[] = {0.102062072615966, | ||
| 422 | 0.102062072615966, | ||
| 423 | 0.058925565098879, | ||
| 424 | 0.102062072615966, | ||
| 425 | 0.058925565098879, | ||
| 426 | 0.058925565098879}; | ||
| 427 | ✗ | return faceToArea[fIdx]; | |
| 428 | } | ||
| 429 | ✗ | else if (geomType.isPyramid()) | |
| 430 | { | ||
| 431 | static const Scalar faceToArea[] = {0.130437298687488, | ||
| 432 | 0.130437298687488, | ||
| 433 | 0.130437298687488, | ||
| 434 | 0.130437298687488, | ||
| 435 | 0.150923085635624, | ||
| 436 | 0.1092906420717, | ||
| 437 | 0.1092906420717, | ||
| 438 | 0.0781735959970571}; | ||
| 439 | ✗ | return faceToArea[fIdx]; | |
| 440 | } | ||
| 441 | ✗ | else if (geomType.isPrism()) | |
| 442 | { | ||
| 443 | static const Scalar faceToArea[] = {0.166666666666667, | ||
| 444 | 0.166666666666667, | ||
| 445 | 0.166666666666667, | ||
| 446 | 0.186338998124982, | ||
| 447 | 0.186338998124982, | ||
| 448 | 0.117851130197758, | ||
| 449 | 0.186338998124982, | ||
| 450 | 0.186338998124982, | ||
| 451 | 0.117851130197758}; | ||
| 452 | ✗ | return faceToArea[fIdx]; | |
| 453 | } | ||
| 454 | else { | ||
| 455 | ✗ | DUNE_THROW(Dune::NotImplemented, "scvf area for unknown GeometryType"); | |
| 456 | } | ||
| 457 | } | ||
| 458 | |||
| 459 | private: | ||
| 460 | const Problem& problem_; | ||
| 461 | const GridGeometry& gridGeometry_; | ||
| 462 | const GridVariables& gridVariables_; | ||
| 463 | |||
| 464 | std::vector<int> cellNum_; | ||
| 465 | }; | ||
| 466 | |||
| 467 | } // end namespace Dumux | ||
| 468 | |||
| 469 | #endif | ||
| 470 |