GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/porousmediumflow/compositional/localresidual.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 52 52 100.0%
Functions: 188 190 98.9%
Branches: 72 103 69.9%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-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 Element-wise calculation of the local residual for problems
11 * using compositional fully implicit model.
12 */
13
14 #ifndef DUMUX_COMPOSITIONAL_LOCAL_RESIDUAL_HH
15 #define DUMUX_COMPOSITIONAL_LOCAL_RESIDUAL_HH
16
17 #include <vector>
18 #include <dune/common/exceptions.hh>
19 #include <dumux/common/properties.hh>
20 #include <dumux/common/numeqvector.hh>
21 #include <dumux/discretization/method.hh>
22 #include <dumux/discretization/defaultlocaloperator.hh>
23 #include <dumux/flux/referencesystemformulation.hh>
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup PorousmediumflowModels
29 * \brief Element-wise calculation of the local residual for problems
30 * using compositional fully implicit model.
31 */
32 template<class TypeTag>
33 class CompositionalLocalResidual
34 : public DiscretizationDefaultLocalOperator<TypeTag>
35 {
36 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
37 using ParentType = DiscretizationDefaultLocalOperator<TypeTag>;
38 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
39 using Problem = GetPropType<TypeTag, Properties::Problem>;
40 using FVElementGeometry = typename GridGeometry::LocalView;
41 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
42 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
43 using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
44 using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
45 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
46 using GridView = typename GridGeometry::GridView;
47 using Element = typename GridView::template Codim<0>::Entity;
48 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
49 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
50 using EnergyLocalResidual = GetPropType<TypeTag, Properties::EnergyLocalResidual>;
51 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
52 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
53 using Indices = typename ModelTraits::Indices;
54
55 static constexpr int numPhases = ModelTraits::numFluidPhases();
56 static constexpr int numComponents = ModelTraits::numFluidComponents();
57 static constexpr bool useMoles = ModelTraits::useMoles();
58
59 enum { conti0EqIdx = Indices::conti0EqIdx };
60
61 //! The index of the component balance equation that gets replaced with the total mass balance
62 static constexpr int replaceCompEqIdx = ModelTraits::replaceCompEqIdx();
63 static constexpr bool useTotalMoleOrMassBalance = replaceCompEqIdx < numComponents;
64
65 public:
66
3/6
✓ Branch 1 taken 8670888 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 118658 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 516152 times.
✗ Branch 8 not taken.
 9998504 using ParentType::ParentType;
67
68 /*!
69 * \brief Evaluates the amount of all conservation quantities
70 * (e.g. phase mass) within a sub-control volume.
71 *
72 * The result should be averaged over the volume (e.g. phase mass
73 * inside a sub control volume divided by the volume)
74 *
75 * \param problem The problem
76 * \param scv The sub control volume
77 * \param volVars The current or previous volVars
78 */
79 424410848 NumEqVector computeStorage(const Problem& problem,
80 const SubControlVolume& scv,
81 const VolumeVariables& volVars) const
82 {
83 424410848 NumEqVector storage(0.0);
84
85 1691216920 const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
86 1707710416 { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
87
88 1674281224 const auto massOrMoleFraction= [](const auto& volVars, const int phaseIdx, const int compIdx)
89 1674281224 { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
90
91 // compute storage term of all components within all phases
92
8/8
✓ Branch 0 taken 745724842 times.
✓ Branch 1 taken 400768986 times.
✓ Branch 2 taken 6380106 times.
✓ Branch 3 taken 6380106 times.
✓ Branch 4 taken 3466380 times.
✓ Branch 5 taken 3466380 times.
✓ Branch 6 taken 3289104 times.
✓ Branch 7 taken 3289104 times.
 1184108342 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
93 {
94
8/8
✓ Branch 0 taken 1659518012 times.
✓ Branch 1 taken 745724842 times.
✓ Branch 2 taken 13668892 times.
✓ Branch 3 taken 6380106 times.
✓ Branch 4 taken 6932760 times.
✓ Branch 5 taken 3466380 times.
✓ Branch 6 taken 6578208 times.
✓ Branch 7 taken 3289104 times.
 2477077120 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
95 {
96 1559984896 auto eqIdx = conti0EqIdx + compIdx;
97
2/2
✓ Branch 0 taken 51073992 times.
✓ Branch 1 taken 33429192 times.
 84503184 if (eqIdx != replaceCompEqIdx)
98 1856410456 storage[eqIdx] += volVars.porosity()
99 1674281224 * volVars.saturation(phaseIdx)
100 1674281224 * massOrMoleDensity(volVars, phaseIdx)
101 1674281224 * massOrMoleFraction(volVars, phaseIdx, compIdx);
102 }
103
104 // in case one balance is substituted by the total mole balance
105 if (useTotalMoleOrMassBalance)
106
2/2
✓ Branch 0 taken 35289600 times.
✓ Branch 1 taken 17644800 times.
 68718792 storage[replaceCompEqIdx] += massOrMoleDensity(volVars, phaseIdx)
107 33429192 * volVars.porosity()
108 33429192 * volVars.saturation(phaseIdx);
109
110 //! The energy storage in the fluid phase with index phaseIdx
111 456184884 EnergyLocalResidual::fluidPhaseStorage(storage, problem, scv, volVars, phaseIdx);
112 }
113
114 //! The energy storage in the solid matrix
115 249909332 EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
116
117 421941096 return storage;
118 }
119
120 /*!
121 * \brief Evaluates the total flux of all conservation quantities
122 * over a face of a sub-control volume.
123 *
124 * \param problem The problem
125 * \param element The current element.
126 * \param fvGeometry The finite-volume geometry
127 * \param elemVolVars The volume variables of the current element
128 * \param scvf The sub control volume face to compute the flux on
129 * \param elemFluxVarsCache The cache related to flux computation
130 */
131 416929764 NumEqVector computeFlux(const Problem& problem,
132 const Element& element,
133 const FVElementGeometry& fvGeometry,
134 const ElementVolumeVariables& elemVolVars,
135 const SubControlVolumeFace& scvf,
136 const ElementFluxVariablesCache& elemFluxVarsCache) const
137 {
138 416929764 FluxVariables fluxVars;
139 416929764 fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
140 static constexpr auto referenceSystemFormulationDiffusion = FluxVariables::MolecularDiffusionType::referenceSystemFormulation();
141 // get upwind weights into local scope
142 416929764 NumEqVector flux(0.0);
143
144 2952783192 const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
145 2982830040 { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
146
147
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 81734568 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 130931032 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3052536 times.
 2890608536 const auto massOrMoleFraction = [](const auto& volVars, const int phaseIdx, const int compIdx)
148
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13443336 times.
 2897929784 { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
149
150 // Check for the reference system and adapt units of the flux accordingly.
151 1575748460 const auto adaptFluxUnits = [](const Scalar referenceFlux, const Scalar molarMass,
152 const ReferenceSystemFormulation referenceSystemFormulation)
153 {
154 if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged)
155 return useMoles ? referenceFlux/molarMass
156 448550828 : referenceFlux;
157 else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged)
158 return useMoles ? referenceFlux
159 : referenceFlux*molarMass;
160 else
161 DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented");
162 };
163
164
2/2
✓ Branch 0 taken 723658610 times.
✓ Branch 1 taken 378593292 times.
 1178924846 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
165 {
166
1/2
✓ Branch 1 taken 366780850 times.
✗ Branch 2 not taken.
 761995082 const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
167 3044800 auto dispersiveFluxes = decltype(diffusiveFluxes)(0.0);
168 if constexpr (ModelTraits::enableCompositionalDispersion())
169
1/2
✓ Branch 1 taken 961600 times.
✗ Branch 2 not taken.
 10276800 dispersiveFluxes = fluxVars.compositionalDispersionFlux(phaseIdx);
170
171
3/3
✓ Branch 0 taken 769361712 times.
✓ Branch 1 taken 1099697220 times.
✓ Branch 2 taken 353675194 times.
 2337743542 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
172 {
173 // get equation index
174 1575748460 const auto eqIdx = conti0EqIdx + compIdx;
175
176 // the physical quantities for which we perform upwinding
177
2/4
✓ Branch 0 taken 8703666 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 6319758 times.
✗ Branch 3 not taken.
 4433112772 const auto upwindTerm = [&massOrMoleDensity, &massOrMoleFraction, phaseIdx, compIdx] (const auto& volVars)
178
16/26
✓ Branch 0 taken 8703666 times.
✓ Branch 1 taken 92361456 times.
✓ Branch 2 taken 8703666 times.
✓ Branch 3 taken 130931032 times.
✓ Branch 4 taken 8703666 times.
✓ Branch 5 taken 3052536 times.
✓ Branch 6 taken 8703666 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 8703666 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 8703666 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 6319758 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 6319758 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6319758 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 6319758 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 6319758 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 6319758 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 2816448 times.
 2973046904 { return massOrMoleDensity(volVars, phaseIdx)*massOrMoleFraction(volVars, phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
179
180 // advective fluxes (only for the component balances)
181
2/2
✓ Branch 0 taken 35121760 times.
✓ Branch 1 taken 23067680 times.
 66209088 if (eqIdx != replaceCompEqIdx)
182
4/5
✓ Branch 1 taken 1408004368 times.
✓ Branch 2 taken 16917600 times.
✓ Branch 3 taken 2895 times.
✓ Branch 4 taken 1189066215 times.
✗ Branch 5 not taken.
 1548670956 flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
183
184 // diffusive and dispersive fluxes
185
3/3
✓ Branch 1 taken 1209170600 times.
✓ Branch 2 taken 12054080 times.
✓ Branch 0 taken 11013600 times.
 1575748460 auto diffusiveAndDispersiveFluxes = adaptFluxUnits(diffusiveFluxes[compIdx],
186 FluidSystem::molarMass(compIdx),
187 referenceSystemFormulationDiffusion);
188 if constexpr (ModelTraits::enableCompositionalDispersion())
189 {
190 static constexpr auto referenceSystemFormulationDispersion =
191 FluxVariables::DispersionFluxType::referenceSystemFormulation();
192 20553600 diffusiveAndDispersiveFluxes += adaptFluxUnits(dispersiveFluxes[compIdx],
193 FluidSystem::molarMass(compIdx),
194 referenceSystemFormulationDispersion);
195 }
196
2/2
✓ Branch 0 taken 35121760 times.
✓ Branch 1 taken 23067680 times.
 66209088 if(eqIdx != replaceCompEqIdx)
197 1548670956 flux[eqIdx] += diffusiveAndDispersiveFluxes;
198 if (useTotalMoleOrMassBalance)
199 flux[replaceCompEqIdx] += diffusiveAndDispersiveFluxes;
200 }
201
202 // in case one balance is substituted by the total mole balance
203 if (useTotalMoleOrMassBalance)
204 {
205 // the physical quantities for which we perform upwinding
206
2/4
✓ Branch 0 taken 8703666 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 6319758 times.
✗ Branch 3 not taken.
 31087328 const auto upwindTerm = [&massOrMoleDensity, phaseIdx] (const auto& volVars)
207
4/8
✓ Branch 0 taken 8703666 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8703666 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6319758 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6319758 times.
✗ Branch 7 not taken.
 27077504 { return massOrMoleDensity(volVars, phaseIdx)*volVars.mobility(phaseIdx); };
208
209
2/4
✓ Branch 1 taken 23067680 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 400080 times.
✗ Branch 5 not taken.
 27077504 flux[replaceCompEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
210 }
211
212 //! Add advective phase energy fluxes. For isothermal model the contribution is zero.
213
1/2
✓ Branch 1 taken 301625597 times.
✗ Branch 2 not taken.
 327744579 EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
214 }
215
216 //! Add diffusive and dispersive energy fluxes. For isothermal model the contribution is zero.
217 416929764 EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
218 5622400 EnergyLocalResidual::heatDispersionFlux(flux, fluxVars);
219
220 416929764 return flux;
221 }
222 };
223
224 } // end namespace Dumux
225
226 #endif
227