| 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 TwoPModel | ||
| 10 | * \brief Performs the transfer of data on a grid from before to after adaptation. | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_TWOP_GRIDDATA_TRANSFER_HH | ||
| 14 | #define DUMUX_TWOP_GRIDDATA_TRANSFER_HH | ||
| 15 | |||
| 16 | #include <memory> | ||
| 17 | |||
| 18 | #include <dune/grid/common/partitionset.hh> | ||
| 19 | #include <dune/grid/utility/persistentcontainer.hh> | ||
| 20 | |||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | |||
| 23 | #include <dumux/discretization/method.hh> | ||
| 24 | #include <dumux/discretization/extrusion.hh> | ||
| 25 | #include <dumux/discretization/elementsolution.hh> | ||
| 26 | #include <dumux/geometry/volume.hh> | ||
| 27 | #include <dumux/porousmediumflow/2p/formulation.hh> | ||
| 28 | #include <dumux/adaptive/griddatatransfer.hh> | ||
| 29 | |||
| 30 | namespace Dumux { | ||
| 31 | |||
| 32 | /*! | ||
| 33 | * \ingroup TwoPModel | ||
| 34 | * \brief Class performing the transfer of data on a grid from before to after adaptation. | ||
| 35 | */ | ||
| 36 | template<class TypeTag> | ||
| 37 | class TwoPGridDataTransfer : public GridDataTransfer<GetPropType<TypeTag, Properties::Grid>> | ||
| 38 | { | ||
| 39 | using Grid = GetPropType<TypeTag, Properties::Grid>; | ||
| 40 | using ParentType = GridDataTransfer<Grid>; | ||
| 41 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 42 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 43 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 44 | using GridView = typename GridGeometry::GridView; | ||
| 45 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 46 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 47 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 48 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 49 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 50 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 51 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 52 | using Element = typename Grid::template Codim<0>::Entity; | ||
| 53 | using ElementSolution = std::decay_t<decltype(elementSolution(std::declval<Element>(), | ||
| 54 | std::declval<SolutionVector>(), | ||
| 55 | std::declval<GridGeometry>()))>; | ||
| 56 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 57 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 58 | using Indices = typename ModelTraits::Indices; | ||
| 59 | |||
| 60 | struct AdaptedValues | ||
| 61 | { | ||
| 62 | 65 | AdaptedValues() : associatedMass(0.0) {} | |
| 63 | ElementSolution u; | ||
| 64 | int count = 0; | ||
| 65 | PrimaryVariables associatedMass; | ||
| 66 | bool wasLeaf = false; | ||
| 67 | }; | ||
| 68 | |||
| 69 | using PersistentContainer = Dune::PersistentContainer<Grid, AdaptedValues>; | ||
| 70 | |||
| 71 | static constexpr int dim = Grid::dimension; | ||
| 72 | static constexpr int dimWorld = Grid::dimensionworld; | ||
| 73 | static constexpr bool isBox = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethods::box; | ||
| 74 | |||
| 75 | // saturation primary variable index | ||
| 76 | enum { saturationIdx = Indices::saturationIdx }; | ||
| 77 | |||
| 78 | // phase indices | ||
| 79 | enum | ||
| 80 | { | ||
| 81 | phase0Idx = FluidSystem::phase0Idx, | ||
| 82 | phase1Idx = FluidSystem::phase1Idx, | ||
| 83 | }; | ||
| 84 | |||
| 85 | // formulations | ||
| 86 | static constexpr auto p0s1 = TwoPFormulation::p0s1; | ||
| 87 | static constexpr auto p1s0 = TwoPFormulation::p1s0; | ||
| 88 | |||
| 89 | // the formulation that is actually used | ||
| 90 | static constexpr auto formulation = ModelTraits::priVarFormulation(); | ||
| 91 | |||
| 92 | // This won't work (mass conservative) for compressible fluids | ||
| 93 | static_assert(!FluidSystem::isCompressible(phase0Idx) | ||
| 94 | && !FluidSystem::isCompressible(phase1Idx), | ||
| 95 | "This adaption helper is only mass conservative for incompressible fluids!"); | ||
| 96 | |||
| 97 | // check if the used formulation is implemented here | ||
| 98 | static_assert(formulation == p0s1 || formulation == p1s0, "Chosen formulation not known to the TwoPGridDataTransfer"); | ||
| 99 | |||
| 100 | public: | ||
| 101 | /*! | ||
| 102 | * \brief Constructor | ||
| 103 | * | ||
| 104 | * \param problem The DuMuX problem to be solved | ||
| 105 | * \param gridGeometry The finite volume grid geometry | ||
| 106 | * \param gridVariables The secondary variables on the grid | ||
| 107 | * \param sol The solution (primary variables) on the grid | ||
| 108 | */ | ||
| 109 | 5 | TwoPGridDataTransfer(std::shared_ptr<const Problem> problem, | |
| 110 | std::shared_ptr<GridGeometry> gridGeometry, | ||
| 111 | std::shared_ptr<const GridVariables> gridVariables, | ||
| 112 | SolutionVector& sol) | ||
| 113 | : ParentType() | ||
| 114 | , problem_(problem) | ||
| 115 | , gridGeometry_(gridGeometry) | ||
| 116 | , gridVariables_(gridVariables) | ||
| 117 | , sol_(sol) | ||
| 118 |
4/14✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 5 times.
✓ Branch 10 taken 5 times.
✗ 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.
|
20 | , adaptionMap_(gridGeometry->gridView().grid(), 0) |
| 119 | 5 | {} | |
| 120 | |||
| 121 | /*! | ||
| 122 | * \brief Stores primary variables and additional data | ||
| 123 | * | ||
| 124 | * To reconstruct the solution in father elements, problem properties might | ||
| 125 | * need to be accessed. From upper level on downwards, the old solution is stored | ||
| 126 | * into a container object, before the grid is adapted. Father elements hold averaged | ||
| 127 | * information from the son cells for the case of the sons being coarsened. | ||
| 128 | */ | ||
| 129 | 15 | void store(const Grid& grid) override | |
| 130 | { | ||
| 131 | 30 | adaptionMap_.resize(); | |
| 132 | |||
| 133 |
3/3✓ Branch 0 taken 35 times.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 2 times.
|
55 | for (auto level = grid.maxLevel(); level >= 0; level--) |
| 134 | { | ||
| 135 |
2/4✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
|
85 | auto fvGeometry = localView(*gridGeometry_); |
| 136 |
13/18✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 15530 times.
✓ Branch 7 taken 40 times.
✓ Branch 8 taken 15530 times.
✓ Branch 9 taken 5 times.
✓ Branch 10 taken 35 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 112414 times.
✓ Branch 13 taken 35 times.
✓ Branch 14 taken 96884 times.
✓ Branch 15 taken 35 times.
✓ Branch 17 taken 96884 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 96884 times.
✗ Branch 21 not taken.
|
112529 | for (const auto& element : elements(grid.levelGridView(level))) |
| 137 | { | ||
| 138 | // get map entry | ||
| 139 | 112414 | auto& adaptedValues = adaptionMap_[element]; | |
| 140 | |||
| 141 | // put values in the map for leaf elements | ||
| 142 |
5/5✓ Branch 0 taken 12372 times.
✓ Branch 1 taken 100042 times.
✓ Branch 2 taken 12372 times.
✓ Branch 3 taken 80813 times.
✓ Branch 4 taken 19229 times.
|
127944 | if (element.isLeaf()) |
| 143 | { | ||
| 144 |
1/2✓ Branch 1 taken 90027 times.
✗ Branch 2 not taken.
|
90027 | fvGeometry.bindElement(element); |
| 145 | |||
| 146 | // store current element solution | ||
| 147 |
2/4✓ Branch 1 taken 90027 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 90027 times.
✗ Branch 5 not taken.
|
180054 | adaptedValues.u = ElementSolution(element, sol_, *gridGeometry_); |
| 148 | |||
| 149 | // compute mass in the scvs | ||
| 150 |
4/4✓ Branch 0 taken 127125 times.
✓ Branch 1 taken 90027 times.
✓ Branch 2 taken 49470 times.
✓ Branch 3 taken 12372 times.
|
356649 | for (const auto& scv : scvs(fvGeometry)) |
| 151 | { | ||
| 152 |
1/2✓ Branch 1 taken 127125 times.
✗ Branch 2 not taken.
|
127125 | VolumeVariables volVars; |
| 153 |
2/4✓ Branch 1 taken 127125 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 127125 times.
✗ Branch 5 not taken.
|
254250 | volVars.update(adaptedValues.u, *problem_, element, scv); |
| 154 | |||
| 155 | 254250 | const auto poreVolume = Extrusion::volume(fvGeometry, scv)*volVars.porosity(); | |
| 156 | 381375 | adaptedValues.associatedMass[phase1Idx] += poreVolume * volVars.density(phase1Idx) * volVars.saturation(phase1Idx); | |
| 157 | 508500 | adaptedValues.associatedMass[phase0Idx] += poreVolume * volVars.density(phase0Idx) * volVars.saturation(phase0Idx); | |
| 158 | } | ||
| 159 | |||
| 160 | // leaf elements always start with count = 1 | ||
| 161 | 90027 | adaptedValues.count = 1; | |
| 162 | 90027 | adaptedValues.wasLeaf = true; | |
| 163 | } | ||
| 164 | // Average in father elements | ||
| 165 |
5/5✓ Branch 0 taken 12458 times.
✓ Branch 1 taken 99956 times.
✓ Branch 2 taken 12458 times.
✓ Branch 3 taken 79988 times.
✓ Branch 4 taken 19968 times.
|
127944 | if (element.level() > 0) |
| 166 | { | ||
| 167 |
1/2✓ Branch 1 taken 89374 times.
✗ Branch 2 not taken.
|
89374 | auto& adaptedValuesFather = adaptionMap_[element.father()]; |
| 168 | // For some grids the father element is identical to the son element. | ||
| 169 | // In that case averaging is not necessary. | ||
| 170 |
2/2✓ Branch 0 taken 89322 times.
✓ Branch 1 taken 52 times.
|
89374 | if(&adaptedValues != &adaptedValuesFather) |
| 171 | 76916 | storeAdaptionValues(adaptedValues, adaptedValuesFather); | |
| 172 | } | ||
| 173 | |||
| 174 | // The vertices of the non-leaf elements exist on the leaf as well | ||
| 175 | // This element solution constructor uses the vertex mapper to obtain | ||
| 176 | // the privars at the vertices, thus, this works for non-leaf elements! | ||
| 177 |
4/4✓ Branch 0 taken 3158 times.
✓ Branch 1 taken 12372 times.
✓ Branch 2 taken 3158 times.
✓ Branch 3 taken 12372 times.
|
31060 | if(isBox && !element.isLeaf()) |
| 178 |
2/4✓ Branch 1 taken 3158 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3158 times.
✗ Branch 5 not taken.
|
6316 | adaptedValues.u = ElementSolution(element, sol_, *gridGeometry_); |
| 179 | } | ||
| 180 | } | ||
| 181 | 15 | } | |
| 182 | |||
| 183 | /*! | ||
| 184 | * \brief Reconstruct missing primary variables (where elements are created/deleted) | ||
| 185 | * | ||
| 186 | * To reconstruct the solution in father elements, problem properties might | ||
| 187 | * need to be accessed. | ||
| 188 | * Starting from the lowest level, the old solution is mapped on the new grid: | ||
| 189 | * Where coarsened, new cells get information from old father element. | ||
| 190 | * Where refined, a new solution is reconstructed from the old father cell, | ||
| 191 | * and then a new son is created. That is then stored into the general data | ||
| 192 | * structure (AdaptedValues). | ||
| 193 | */ | ||
| 194 | 15 | void reconstruct(const Grid& grid) override | |
| 195 | { | ||
| 196 | 45 | gridGeometry_->update(grid.leafGridView()); | |
| 197 | 15 | reconstruct_(); | |
| 198 | 15 | } | |
| 199 | |||
| 200 | private: | ||
| 201 | |||
| 202 | 15 | void reconstruct_() | |
| 203 | { | ||
| 204 | // resize stuff (grid might have changed) | ||
| 205 | 30 | adaptionMap_.resize(); | |
| 206 | 30 | sol_.resize(gridGeometry_->numDofs()); | |
| 207 | |||
| 208 | // vectors storing the mass associated with each vertex, when using the box method | ||
| 209 |
1/5✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
15 | std::vector<Scalar> massCoeff; |
| 210 |
2/5✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
|
30 | std::vector<Scalar> associatedMass; |
| 211 | |||
| 212 | if(isBox) | ||
| 213 | { | ||
| 214 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
|
6 | massCoeff.resize(gridGeometry_->numDofs(), 0.0); |
| 215 |
3/8✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
6 | associatedMass.resize(gridGeometry_->numDofs(), 0.0); |
| 216 | } | ||
| 217 | |||
| 218 | // iterate over leaf and reconstruct the solution | ||
| 219 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 15 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
45 | auto fvGeometry = localView(*gridGeometry_); |
| 220 |
17/24✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 15 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 15 times.
✓ Branch 7 taken 15 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 15 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 11283 times.
✓ Branch 13 taken 15 times.
✓ Branch 14 taken 11283 times.
✓ Branch 15 taken 2 times.
✓ Branch 16 taken 13 times.
✓ Branch 17 taken 11283 times.
✓ Branch 18 taken 75234 times.
✓ Branch 19 taken 13 times.
✓ Branch 20 taken 86517 times.
✓ Branch 21 taken 13 times.
✓ Branch 23 taken 75234 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 75234 times.
✗ Branch 27 not taken.
|
86562 | for (const auto& element : elements(gridGeometry_->gridView().grid().leafGridView(), Dune::Partitions::interior)) |
| 221 | { | ||
| 222 |
4/5✗ Branch 0 not taken.
✓ Branch 1 taken 86517 times.
✓ Branch 2 taken 58794 times.
✓ Branch 3 taken 23711 times.
✓ Branch 4 taken 4012 times.
|
86517 | if (!element.isNew()) |
| 223 | { | ||
| 224 | 66065 | const auto& adaptedValues = adaptionMap_[element]; | |
| 225 |
1/2✓ Branch 1 taken 66065 times.
✗ Branch 2 not taken.
|
66065 | fvGeometry.bindElement(element); |
| 226 | |||
| 227 | // obtain element solution from map (divide by count!) | ||
| 228 | 66065 | auto elemSol = adaptedValues.u; | |
| 229 | if (!isBox) | ||
| 230 | 58794 | elemSol[0] /= adaptedValues.count; | |
| 231 | |||
| 232 |
2/4✓ Branch 1 taken 66065 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 66065 times.
✗ Branch 5 not taken.
|
73336 | const auto elementVolume = volume(element.geometry(), Extrusion{}); |
| 233 |
4/4✓ Branch 0 taken 87868 times.
✓ Branch 1 taken 66065 times.
✓ Branch 2 taken 29074 times.
✓ Branch 3 taken 7271 times.
|
249072 | for (const auto& scv : scvs(fvGeometry)) |
| 234 | { | ||
| 235 |
1/2✓ Branch 1 taken 87868 times.
✗ Branch 2 not taken.
|
87868 | VolumeVariables volVars; |
| 236 |
2/4✓ Branch 1 taken 87868 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87868 times.
✗ Branch 5 not taken.
|
175736 | volVars.update(elemSol, *problem_, element, scv); |
| 237 | |||
| 238 | // write solution at dof in current solution vector | ||
| 239 |
4/4✓ Branch 0 taken 4917 times.
✓ Branch 1 taken 53877 times.
✓ Branch 2 taken 4917 times.
✓ Branch 3 taken 53877 times.
|
175736 | sol_[scv.dofIndex()] = elemSol[scv.localDofIndex()]; |
| 240 | |||
| 241 |
2/2✓ Branch 0 taken 4917 times.
✓ Branch 1 taken 53877 times.
|
87868 | const auto dofIdxGlobal = scv.dofIndex(); |
| 242 | // For cc schemes, overwrite the saturation by a mass conservative one here | ||
| 243 | if (!isBox) | ||
| 244 | { | ||
| 245 | // only recalculate the saturations if element hasn't been leaf before adaptation | ||
| 246 |
2/2✓ Branch 0 taken 4917 times.
✓ Branch 1 taken 53877 times.
|
58794 | if (!adaptedValues.wasLeaf) |
| 247 | { | ||
| 248 | if (formulation == p0s1) | ||
| 249 | { | ||
| 250 | 9834 | sol_[dofIdxGlobal][saturationIdx] = adaptedValues.associatedMass[phase1Idx]; | |
| 251 | 19668 | sol_[dofIdxGlobal][saturationIdx] /= elementVolume * volVars.density(phase1Idx) * volVars.porosity(); | |
| 252 | } | ||
| 253 | else if (formulation == p1s0) | ||
| 254 | { | ||
| 255 | sol_[dofIdxGlobal][saturationIdx] = adaptedValues.associatedMass[phase0Idx]; | ||
| 256 | sol_[dofIdxGlobal][saturationIdx] /= elementVolume * volVars.density(phase0Idx) * volVars.porosity(); | ||
| 257 | } | ||
| 258 | } | ||
| 259 | } | ||
| 260 | |||
| 261 | // For the box scheme, add mass & mass coefficient to container (saturations are recalculated at the end) | ||
| 262 | else | ||
| 263 | { | ||
| 264 | 29074 | const auto scvVolume = Extrusion::volume(fvGeometry, scv); | |
| 265 | if (formulation == p0s1) | ||
| 266 | { | ||
| 267 | 87222 | massCoeff[dofIdxGlobal] += scvVolume * volVars.density(phase1Idx) * volVars.porosity(); | |
| 268 | 87222 | associatedMass[dofIdxGlobal] += scvVolume / elementVolume * adaptedValues.associatedMass[phase1Idx]; | |
| 269 | } | ||
| 270 | else if (formulation == p1s0) | ||
| 271 | { | ||
| 272 | massCoeff[dofIdxGlobal] += scvVolume * volVars.density(phase0Idx) * volVars.porosity(); | ||
| 273 | associatedMass[dofIdxGlobal] += scvVolume / elementVolume * adaptedValues.associatedMass[phase0Idx]; | ||
| 274 | } | ||
| 275 | } | ||
| 276 | } | ||
| 277 | } | ||
| 278 | else | ||
| 279 | { | ||
| 280 | // value is not in map, interpolate from father element | ||
| 281 |
3/5✗ Branch 0 not taken.
✓ Branch 1 taken 20452 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4012 times.
✓ Branch 4 taken 16440 times.
|
24464 | assert(element.hasFather() && "new element does not have a father element!"); |
| 282 | |||
| 283 | // find the ancestor element that existed on the old grid already | ||
| 284 |
1/2✓ Branch 1 taken 20452 times.
✗ Branch 2 not taken.
|
36892 | auto fatherElement = element.father(); |
| 285 |
5/9✗ Branch 0 not taken.
✓ Branch 1 taken 20796 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 16784 times.
✓ Branch 4 taken 4012 times.
✓ Branch 5 taken 344 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 344 times.
✗ Branch 8 not taken.
|
20796 | while(fatherElement.isNew() && fatherElement.level() > 0) |
| 286 |
1/2✓ Branch 1 taken 344 times.
✗ Branch 2 not taken.
|
344 | fatherElement = fatherElement.father(); |
| 287 | |||
| 288 | if(!isBox) | ||
| 289 | { | ||
| 290 | 16440 | const auto& adaptedValuesFather = adaptionMap_[fatherElement]; | |
| 291 | |||
| 292 | // obtain the mass contained in father | ||
| 293 | 16440 | Scalar massFather = 0.0; | |
| 294 | if (formulation == p0s1) | ||
| 295 | 16440 | massFather = adaptedValuesFather.associatedMass[phase1Idx]; | |
| 296 | else if (formulation == p1s0) | ||
| 297 | massFather = adaptedValuesFather.associatedMass[phase0Idx]; | ||
| 298 | |||
| 299 | // obtain the element solution through the father | ||
| 300 | 16440 | auto elemSolSon = adaptedValuesFather.u; | |
| 301 | 16440 | elemSolSon[0] /= adaptedValuesFather.count; | |
| 302 | |||
| 303 |
1/2✓ Branch 1 taken 16440 times.
✗ Branch 2 not taken.
|
16440 | fvGeometry.bindElement(element); |
| 304 | |||
| 305 |
2/2✓ Branch 0 taken 16440 times.
✓ Branch 1 taken 16440 times.
|
32880 | for (const auto& scv : scvs(fvGeometry)) |
| 306 | { | ||
| 307 |
1/2✓ Branch 1 taken 16440 times.
✗ Branch 2 not taken.
|
16440 | VolumeVariables volVars; |
| 308 |
2/4✓ Branch 1 taken 16440 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16440 times.
✗ Branch 5 not taken.
|
32880 | volVars.update(elemSolSon, *problem_, element, scv); |
| 309 | |||
| 310 | // store constructed values of son in the current solution | ||
| 311 |
2/4✓ Branch 1 taken 16440 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16440 times.
✗ Branch 5 not taken.
|
32880 | sol_[scv.dofIndex()] = elemSolSon[0]; |
| 312 | |||
| 313 | // overwrite the saturation by a mass conservative one here | ||
| 314 | 16440 | Scalar massCoeffSon = 0.0; | |
| 315 | if (formulation == p0s1) | ||
| 316 |
3/6✓ Branch 1 taken 16440 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16440 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16440 times.
✗ Branch 8 not taken.
|
49320 | massCoeffSon = Extrusion::volume(fvGeometry, scv) * volVars.density(phase1Idx) * volVars.porosity(); |
| 317 | else if (formulation == p1s0) | ||
| 318 | massCoeffSon = Extrusion::volume(fvGeometry, scv) * volVars.density(phase0Idx) * volVars.porosity(); | ||
| 319 | 32880 | sol_[scv.dofIndex()][saturationIdx] = | |
| 320 |
3/6✓ Branch 1 taken 16440 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16440 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16440 times.
✗ Branch 8 not taken.
|
32880 | ( Extrusion::volume(fvGeometry, scv)/volume(fatherElement.geometry(), Extrusion{})*massFather )/massCoeffSon; |
| 321 | } | ||
| 322 | } | ||
| 323 | else | ||
| 324 | { | ||
| 325 | 4012 | auto& adaptedValuesFather = adaptionMap_[fatherElement]; | |
| 326 | |||
| 327 |
1/2✓ Branch 1 taken 4012 times.
✗ Branch 2 not taken.
|
4012 | fvGeometry.bindElement(element); |
| 328 | |||
| 329 | // interpolate solution in the father to the vertices of the new son | ||
| 330 |
2/4✓ Branch 1 taken 4012 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4012 times.
✗ Branch 5 not taken.
|
8024 | ElementSolution elemSolSon(element, sol_, *gridGeometry_); |
| 331 | 4012 | const auto fatherGeometry = fatherElement.geometry(); | |
| 332 |
4/4✓ Branch 0 taken 15875 times.
✓ Branch 1 taken 4012 times.
✓ Branch 2 taken 15875 times.
✓ Branch 3 taken 4012 times.
|
39774 | for (const auto& scv : scvs(fvGeometry)) |
| 333 |
1/2✓ Branch 1 taken 15875 times.
✗ Branch 2 not taken.
|
15875 | elemSolSon[scv.localDofIndex()] = evalSolution(fatherElement, |
| 334 | fatherGeometry, | ||
| 335 |
1/2✓ Branch 1 taken 15875 times.
✗ Branch 2 not taken.
|
15875 | adaptedValuesFather.u, |
| 336 | scv.dofPosition()); | ||
| 337 | |||
| 338 | // compute mass & mass coefficients for the scvs (saturations are recalculated at the end) | ||
| 339 |
1/2✓ Branch 1 taken 4012 times.
✗ Branch 2 not taken.
|
4012 | const auto fatherElementVolume = volume(fatherGeometry, Extrusion{}); |
| 340 |
4/4✓ Branch 0 taken 15875 times.
✓ Branch 1 taken 4012 times.
✓ Branch 2 taken 15875 times.
✓ Branch 3 taken 4012 times.
|
39774 | for (const auto& scv : scvs(fvGeometry)) |
| 341 | { | ||
| 342 |
1/2✓ Branch 1 taken 15875 times.
✗ Branch 2 not taken.
|
15875 | VolumeVariables volVars; |
| 343 |
2/4✓ Branch 1 taken 15875 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15875 times.
✗ Branch 5 not taken.
|
31750 | volVars.update(elemSolSon, *problem_, element, scv); |
| 344 | |||
| 345 | 15875 | const auto dofIdxGlobal = scv.dofIndex(); | |
| 346 | 15875 | const auto scvVolume = Extrusion::volume(fvGeometry, scv); | |
| 347 | if (formulation == p0s1) | ||
| 348 | { | ||
| 349 | 47625 | massCoeff[dofIdxGlobal] += scvVolume * volVars.density(phase1Idx) * volVars.porosity(); | |
| 350 | 31750 | associatedMass[dofIdxGlobal] += scvVolume / fatherElementVolume * adaptedValuesFather.associatedMass[phase1Idx]; | |
| 351 | } | ||
| 352 | else if (formulation == p1s0) | ||
| 353 | { | ||
| 354 | massCoeff[dofIdxGlobal] += scvVolume * volVars.density(phase0Idx) * volVars.porosity(); | ||
| 355 | associatedMass[dofIdxGlobal] += scvVolume / fatherElementVolume * adaptedValuesFather.associatedMass[phase0Idx]; | ||
| 356 | } | ||
| 357 | |||
| 358 | // store constructed (pressure) values of son in the current solution (saturation comes later) | ||
| 359 | 47625 | sol_[dofIdxGlobal] = elemSolSon[scv.localDofIndex()]; | |
| 360 | } | ||
| 361 | } | ||
| 362 | } | ||
| 363 | } | ||
| 364 | |||
| 365 | if(isBox) | ||
| 366 | { | ||
| 367 |
6/6✓ Branch 0 taken 11457 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 11457 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 11457 times.
✓ Branch 5 taken 2 times.
|
34373 | for(std::size_t dofIdxGlobal = 0; dofIdxGlobal < gridGeometry_->numDofs(); dofIdxGlobal++) |
| 368 | 57285 | sol_[dofIdxGlobal][saturationIdx] = associatedMass[dofIdxGlobal] / massCoeff[dofIdxGlobal]; | |
| 369 | } | ||
| 370 | |||
| 371 | // reset entries in adaptation map | ||
| 372 |
1/2✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
|
30 | adaptionMap_.resize( typename PersistentContainer::Value() ); |
| 373 | 15 | adaptionMap_.shrinkToFit(); | |
| 374 | 43 | adaptionMap_.fill( typename PersistentContainer::Value() ); | |
| 375 | |||
| 376 | //! TODO: fix adaptive simulations in parallel | ||
| 377 | //#if HAVE_MPI | ||
| 378 | // // communicate ghost data | ||
| 379 | // using SolutionTypes = typename GetProp<TypeTag, SolutionTypes>; | ||
| 380 | // using ElementMapper = typename SolutionTypes::ElementMapper; | ||
| 381 | // using DataHandle = VectorExchange<ElementMapper, std::vector<CellData> >; | ||
| 382 | // DataHandle dataHandle(problem.elementMapper(), this->cellDataGlobal()); | ||
| 383 | // problem.gridView().template communicate<DataHandle>(dataHandle, | ||
| 384 | // Dune::InteriorBorder_All_Interface, | ||
| 385 | // Dune::ForwardCommunication); | ||
| 386 | //#endif | ||
| 387 | 15 | } | |
| 388 | |||
| 389 | /*! | ||
| 390 | * \brief Stores sons entries into father element for averaging | ||
| 391 | * | ||
| 392 | * Sum up the adaptedValues (sons values) into father element. We store from leaf | ||
| 393 | * upwards, so sons are stored first, then cells on the next leaf (=fathers) | ||
| 394 | * can be averaged. | ||
| 395 | * | ||
| 396 | * \param adaptedValues Container for model-specific values to be adapted | ||
| 397 | * \param adaptedValuesFather Values to be adapted of father cell | ||
| 398 | */ | ||
| 399 | 76916 | static void storeAdaptionValues(AdaptedValues& adaptedValues, | |
| 400 | AdaptedValues& adaptedValuesFather) | ||
| 401 | { | ||
| 402 | // Add associated mass of the child to the one of the father | ||
| 403 | 89322 | adaptedValuesFather.associatedMass += adaptedValues.associatedMass; | |
| 404 | |||
| 405 | if(!isBox) | ||
| 406 | { | ||
| 407 | // add the child's primary variables to the ones of father | ||
| 408 | // we have to divide the child's ones in case it was composed | ||
| 409 | // of several children as well! | ||
| 410 | 76916 | auto values = adaptedValues.u[0]; | |
| 411 | 76916 | values /= adaptedValues.count; | |
| 412 | 76916 | adaptedValuesFather.u[0] += values; | |
| 413 | |||
| 414 | // keep track of the number of children that composed this father | ||
| 415 | 76916 | adaptedValuesFather.count += 1; | |
| 416 | |||
| 417 | // A father element is never leaf | ||
| 418 | 76916 | adaptedValuesFather.wasLeaf = false; | |
| 419 | } | ||
| 420 | else | ||
| 421 | { | ||
| 422 | // For the box scheme, scaling of primary variables by count is obsolete | ||
| 423 | // Thus, we always want count = 1 | ||
| 424 | 12406 | adaptedValuesFather.count = 1; | |
| 425 | |||
| 426 | // A father element is never leaf | ||
| 427 | 12406 | adaptedValuesFather.wasLeaf = false; | |
| 428 | } | ||
| 429 | 76916 | } | |
| 430 | |||
| 431 | std::shared_ptr<const Problem> problem_; | ||
| 432 | std::shared_ptr<GridGeometry> gridGeometry_; | ||
| 433 | std::shared_ptr<const GridVariables> gridVariables_; | ||
| 434 | SolutionVector& sol_; | ||
| 435 | PersistentContainer adaptionMap_; | ||
| 436 | }; | ||
| 437 | |||
| 438 | } // end namespace Dumux | ||
| 439 | |||
| 440 | #endif | ||
| 441 |