| 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 CCMpfaDiscretization | ||
| 10 | * \brief Class providing functionality for the reconstruction of the | ||
| 11 | * gradients in the sub-control volumes involved in mpfa schemes. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_CC_MPFA_SCV_GRADIENTS_HH | ||
| 14 | #define DUMUX_CC_MPFA_SCV_GRADIENTS_HH | ||
| 15 | |||
| 16 | #include <vector> | ||
| 17 | #include <utility> | ||
| 18 | #include <type_traits> | ||
| 19 | |||
| 20 | #include <dune/common/fvector.hh> | ||
| 21 | |||
| 22 | #include <dumux/common/math.hh> | ||
| 23 | #include <dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh> | ||
| 24 | |||
| 25 | namespace Dumux { | ||
| 26 | |||
| 27 | /*! | ||
| 28 | * \ingroup CCMpfaDiscretization | ||
| 29 | * \brief Class providing functionality for the reconstruction of the | ||
| 30 | * gradients in the sub-control volumes involved in mpfa schemes. | ||
| 31 | */ | ||
| 32 | class CCMpfaScvGradients | ||
| 33 | { | ||
| 34 | //! Return type of the gradient computation function (pair of scv centers and corresponding gradients) | ||
| 35 | template<class GridGeometry, class Scalar> | ||
| 36 | using ResultPair = std::pair< std::vector<typename GridGeometry::SubControlVolume::GlobalPosition>, | ||
| 37 | std::vector<Dune::FieldVector<Scalar, GridGeometry::GridView::dimension>> >; | ||
| 38 | |||
| 39 | public: | ||
| 40 | |||
| 41 | /*! | ||
| 42 | * \brief Computes the phase velocities in the scvs of the grid. | ||
| 43 | * | ||
| 44 | * \param gridGeometry The finite volume grid geometry | ||
| 45 | * \param gridVariables The variables living on the grid | ||
| 46 | * \param x The vector containing the solution | ||
| 47 | * \param phaseIdx The index of the fluid phase to be considered | ||
| 48 | */ | ||
| 49 | template<class GridGeometry, class GridVariables, class SolutionVector> | ||
| 50 | static ResultPair<GridGeometry, typename GridVariables::Scalar> | ||
| 51 | computeVelocities(const GridGeometry& gridGeometry, | ||
| 52 | const GridVariables& gridVariables, | ||
| 53 | const SolutionVector& x, | ||
| 54 | unsigned int phaseIdx) | ||
| 55 | { | ||
| 56 | 801 | auto gradToVelocity = [phaseIdx] (const auto& grad, const auto& volVars) | |
| 57 | { | ||
| 58 | 400 | auto vel = mv(volVars.permeability(), grad); | |
| 59 | 800 | vel *= volVars.mobility(phaseIdx); | |
| 60 | 400 | return vel; | |
| 61 | }; | ||
| 62 | |||
| 63 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | return computePressureGradients(gridGeometry, gridVariables, x, phaseIdx, gradToVelocity); |
| 64 | } | ||
| 65 | |||
| 66 | /*! | ||
| 67 | * \brief Computes the pressure gradients in the scvs of the grid. | ||
| 68 | * | ||
| 69 | * \param gridGeometry The finite volume grid geometry | ||
| 70 | * \param gridVariables The variables living on the grid | ||
| 71 | * \param x The vector containing the solution | ||
| 72 | * \param phaseIdx The index of the fluid phase to be considered | ||
| 73 | */ | ||
| 74 | template<class GridGeometry, class GridVariables, class SolutionVector> | ||
| 75 | static ResultPair<GridGeometry, typename GridVariables::Scalar> | ||
| 76 | computePressureGradients(const GridGeometry& gridGeometry, | ||
| 77 | const GridVariables& gridVariables, | ||
| 78 | const SolutionVector& x, | ||
| 79 | unsigned int phaseIdx) | ||
| 80 | { | ||
| 81 | auto f = [] (const auto& grad, const auto& volVars) { return grad; }; | ||
| 82 | return computePressureGradients(gridGeometry, gridVariables, x, phaseIdx, f); | ||
| 83 | } | ||
| 84 | |||
| 85 | /*! | ||
| 86 | * \brief Computes the pressure gradients in the scvs of the grid. | ||
| 87 | * | ||
| 88 | * \param gridGeometry The finite volume grid geometry | ||
| 89 | * \param gridVariables The variables living on the grid | ||
| 90 | * \param x The vector containing the solution | ||
| 91 | * \param phaseIdx The index of the fluid phase to be considered | ||
| 92 | * \param f a function which is applied to the gradients and which | ||
| 93 | * has the following signature: | ||
| 94 | * | ||
| 95 | * Gradient f(const Gradient& gradient, const VolumeVariables& volVars) | ||
| 96 | * | ||
| 97 | * It receives the scv-gradient and the corresponding volume | ||
| 98 | * variables and returns the modified gradient. This can be | ||
| 99 | * used e.g. to turn the pressure gradients into velocities. | ||
| 100 | */ | ||
| 101 | template<class GridGeometry, class GridVariables, class SolutionVector, class F> | ||
| 102 | static ResultPair<GridGeometry, typename GridVariables::Scalar> | ||
| 103 | computePressureGradients(const GridGeometry& gridGeometry, | ||
| 104 | const GridVariables& gridVariables, | ||
| 105 | const SolutionVector& x, | ||
| 106 | unsigned int phaseIdx, | ||
| 107 | F& f) | ||
| 108 | { | ||
| 109 | using ElemVolVars = typename GridVariables::GridVolumeVariables::LocalView; | ||
| 110 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 111 | |||
| 112 | 243 | auto handleFunction = [&] (auto& result, | |
| 113 | const auto& handle, | ||
| 114 | const auto& iv, | ||
| 115 | const FVElementGeometry& fvGeometry, | ||
| 116 | 121 | const ElemVolVars& elemVolVars) | |
| 117 | { | ||
| 118 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 121 times.
|
242 | handle.advectionHandle().setPhaseIndex(phaseIdx); |
| 119 |
1/2✓ Branch 2 taken 121 times.
✗ Branch 3 not taken.
|
484 | const auto grads = InteractionVolumeAssemblerHelper::assembleScvGradients(handle.advectionHandle(), iv); |
| 120 |
4/4✓ Branch 0 taken 400 times.
✓ Branch 1 taken 121 times.
✓ Branch 2 taken 400 times.
✓ Branch 3 taken 121 times.
|
1042 | for (unsigned int i = 0; i < iv.numScvs(); ++i) |
| 121 | { | ||
| 122 | 800 | const auto& volVars = elemVolVars[iv.localScv(i).gridScvIndex()]; | |
| 123 |
1/2✓ Branch 1 taken 400 times.
✗ Branch 2 not taken.
|
400 | const auto scvGeometry = iv.getScvGeometry(i, fvGeometry); |
| 124 |
1/2✓ Branch 2 taken 400 times.
✗ Branch 3 not taken.
|
400 | result.first.push_back(scvGeometry.center()); |
| 125 |
1/4✓ Branch 3 taken 400 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
800 | result.second.push_back( f(grads[i], volVars) ); |
| 126 | } | ||
| 127 | }; | ||
| 128 | |||
| 129 | 1 | return computeGradients_(gridGeometry, gridVariables, x, handleFunction); | |
| 130 | } | ||
| 131 | |||
| 132 | private: | ||
| 133 | /*! | ||
| 134 | * \brief Computes the gradients executing the provided function | ||
| 135 | * on an interaction volume / data handle pair. | ||
| 136 | */ | ||
| 137 | template<class GridGeometry, class GridVariables, class SolutionVector, class HandleFunction> | ||
| 138 | static ResultPair<GridGeometry, typename GridVariables::Scalar> | ||
| 139 | 1 | computeGradients_(const GridGeometry& gridGeometry, | |
| 140 | const GridVariables& gridVariables, | ||
| 141 | const SolutionVector& x, | ||
| 142 | const HandleFunction& handleFunction) | ||
| 143 | { | ||
| 144 | using GridView = typename GridGeometry::GridView; | ||
| 145 | static constexpr int dim = GridView::dimension; | ||
| 146 | |||
| 147 | // first, find out how many scvs live on this grid | ||
| 148 | 1 | std::size_t numScvs = 0; | |
| 149 | 1 | const auto& gridView = gridGeometry.gridView(); | |
| 150 |
1/4✓ Branch 1 taken 101 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
101 | for (const auto& element : elements(gridView)) |
| 151 | 200 | numScvs += element.subEntities(dim); | |
| 152 | |||
| 153 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | ResultPair<GridGeometry, typename GridVariables::Scalar> result; |
| 154 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | result.first.reserve(numScvs); |
| 155 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | result.second.reserve(numScvs); |
| 156 |
4/10✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
|
2 | std::vector<bool> vertexHandled(gridView.size(dim), false); |
| 157 | |||
| 158 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | auto fvGeometry = localView(gridGeometry); |
| 159 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
3 | auto elemVolVars = localView(gridVariables.curGridVolVars()); |
| 160 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache()); |
| 161 | |||
| 162 |
2/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 101 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
102 | for (const auto& element : elements(gridView)) |
| 163 | { | ||
| 164 | bool allFinished = true; | ||
| 165 |
4/5✓ Branch 0 taken 400 times.
✓ Branch 1 taken 100 times.
✓ Branch 2 taken 400 times.
✓ Branch 3 taken 100 times.
✗ Branch 4 not taken.
|
600 | for (int i = 0; i < element.subEntities(dim); ++i) |
| 166 |
8/10✓ Branch 1 taken 400 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 400 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 121 times.
✓ Branch 7 taken 279 times.
✓ Branch 8 taken 121 times.
✓ Branch 9 taken 279 times.
✓ Branch 10 taken 121 times.
✓ Branch 11 taken 279 times.
|
800 | if (!vertexHandled[gridGeometry.vertexMapper().subIndex(element, i, dim)]) |
| 167 | 121 | allFinished = false; | |
| 168 | |||
| 169 | // bind views only if there is unfinished business | ||
| 170 |
1/2✓ Branch 0 taken 100 times.
✗ Branch 1 not taken.
|
100 | if (allFinished) |
| 171 | continue; | ||
| 172 | |||
| 173 | // compute gradients in all scvs of all interaction volumes in this element | ||
| 174 |
1/2✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
|
100 | fvGeometry.bind(element); |
| 175 |
1/2✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
|
100 | elemVolVars.bind(element, fvGeometry, x); |
| 176 |
1/2✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
|
100 | elemFluxVarsCache.bind(element, fvGeometry, elemVolVars); |
| 177 | |||
| 178 |
4/4✓ Branch 0 taken 800 times.
✓ Branch 1 taken 100 times.
✓ Branch 2 taken 800 times.
✓ Branch 3 taken 100 times.
|
1000 | for (const auto& scvf : scvfs(fvGeometry)) |
| 179 | { | ||
| 180 |
6/6✓ Branch 0 taken 121 times.
✓ Branch 1 taken 679 times.
✓ Branch 2 taken 121 times.
✓ Branch 3 taken 679 times.
✓ Branch 4 taken 121 times.
✓ Branch 5 taken 679 times.
|
2400 | if (vertexHandled[scvf.vertexIndex()]) |
| 181 | continue; | ||
| 182 | |||
| 183 | 121 | if (gridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) | |
| 184 | { | ||
| 185 | const auto& iv = elemFluxVarsCache.secondaryInteractionVolume(scvf); | ||
| 186 | const auto& handle = elemFluxVarsCache.secondaryDataHandle(scvf); | ||
| 187 | handleFunction(result, handle, iv, fvGeometry, elemVolVars); | ||
| 188 | } | ||
| 189 | else | ||
| 190 | { | ||
| 191 | 121 | const auto& iv = elemFluxVarsCache.primaryInteractionVolume(scvf); | |
| 192 | 121 | const auto& handle = elemFluxVarsCache.primaryDataHandle(scvf); | |
| 193 |
1/2✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
|
121 | handleFunction(result, handle, iv, fvGeometry, elemVolVars); |
| 194 | } | ||
| 195 | |||
| 196 | 363 | vertexHandled[scvf.vertexIndex()] = true; | |
| 197 | } | ||
| 198 | } | ||
| 199 | |||
| 200 | 1 | return result; | |
| 201 | } | ||
| 202 | }; | ||
| 203 | |||
| 204 | } // end namespace Dumux | ||
| 205 | |||
| 206 | #endif | ||
| 207 |