GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/2pncmin/isothermal/problem.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 98 102 96.1%
Functions: 14 18 77.8%
Branches: 101 160 63.1%
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 TwoPNCMinTests
10 * \brief Problem where water is injected in a for flushing precipitated salt clogging a gas reservoir.
11 */
12 #ifndef DUMUX_DISSOLUTION_PROBLEM_HH
13 #define DUMUX_DISSOLUTION_PROBLEM_HH
14
15 #include <dumux/common/properties.hh>
16 #include <dumux/common/parameters.hh>
17 #include <dumux/common/boundarytypes.hh>
18 #include <dumux/common/numeqvector.hh>
19
20 #include <dumux/porousmediumflow/problem.hh>
21
22 namespace Dumux {
23
24 /*!
25 * \ingroup TwoPNCMinTests
26 * \brief Problem where water is injected to flush precipitated salt in a gas reservoir clogged due to precipitated salt.
27 *
28 */
29 template <class TypeTag>
30 class DissolutionProblem : public PorousMediumFlowProblem<TypeTag>
31 {
32 using ParentType = PorousMediumFlowProblem<TypeTag>;
33 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
34 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
35 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
36 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
37 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
38 using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
39
40 enum
41 {
42 // primary variable indices
43 pressureIdx = Indices::pressureIdx,
44 switchIdx = Indices::switchIdx,
45
46 // component indices
47 // TODO: using xwNaClIdx as privaridx works here, but
48 // looks like magic. Can this be done differently??
49 xwNaClIdx = FluidSystem::NaClIdx,
50 precipNaClIdx = FluidSystem::numComponents,
51
52 // Indices of the components
53 H2OIdx = FluidSystem::H2OIdx,
54 NaClIdx = FluidSystem::NaClIdx,
55
56 // Indices of the phases
57 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
58 gasPhaseIdx = FluidSystem::gasPhaseIdx,
59
60 // index of the solid phase
61 sPhaseIdx = SolidSystem::comp0Idx,
62
63
64 // Index of the primary component of G and L phase
65 conti0EqIdx = Indices::conti0EqIdx,
66 precipNaClEqIdx = Indices::conti0EqIdx + FluidSystem::numComponents,
67
68 // Phase State
69 bothPhases = Indices::bothPhases,
70
71 // Grid and world dimension
72 dim = GridView::dimension,
73 dimWorld = GridView::dimensionworld,
74 };
75
76 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
77 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
78 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
79 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
80 using Element = typename GridView::template Codim<0>::Entity;
81 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
82 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
83 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
84 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
85 using GlobalPosition = typename SubControlVolume::GlobalPosition;
86
87 public:
88 2 DissolutionProblem(std::shared_ptr<const GridGeometry> gridGeometry)
89
8/24
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2 times.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
 6 : ParentType(gridGeometry)
90 {
91
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 outerSalinity_ = getParam<Scalar>("Problem.OuterSalinity");
92
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 reservoirPressure_ = getParam<Scalar>("Problem.ReservoirPressure");
93
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 initLiqSaturation_ = getParam<Scalar>("Problem.LiquidSaturation");
94
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 initPrecipitatedSaltBlock_ = getParam<Scalar>("Problem.InitPrecipitatedSaltBlock");
95
96
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 outerLiqSaturation_ = getParam<Scalar>("Problem.OuterLiqSaturation");
97
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 innerLiqSaturation_ = getParam<Scalar>("Problem.InnerLiqSaturation");
98
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 innerSalinity_ = getParam<Scalar>("Problem.InnerSalinity");
99
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 innerPressure_ = getParam<Scalar>("Problem.InnerPressure");
100
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 outerPressure_ = getParam<Scalar>("Problem.OuterPressure");
101
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 reservoirSaturation_ = getParam<Scalar>("Problem.reservoirSaturation");
102
103
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nTemperature_ = getParam<int>("FluidSystem.NTemperature");
104
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 nPressure_ = getParam<int>("FluidSystem.NPressure");
105
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 pressureLow_ = getParam<Scalar>("FluidSystem.PressureLow");
106
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 pressureHigh_ = getParam<Scalar>("FluidSystem.PressureHigh");
107
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 temperatureLow_ = getParam<Scalar>("FluidSystem.TemperatureLow");
108
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 temperatureHigh_ = getParam<Scalar>("FluidSystem.TemperatureHigh");
109
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 2 name_ = getParam<std::string>("Problem.Name");
110
111 2 unsigned int codim = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethods::box ? dim : 0;
112
4/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 10 taken 2 times.
✗ Branch 11 not taken.
 6 permeability_.resize(gridGeometry->gridView().size(codim));
113
114 4 FluidSystem::init(/*Tmin=*/temperatureLow_,
115 /*Tmax=*/temperatureHigh_,
116 2 /*nT=*/nTemperature_,
117 /*pmin=*/pressureLow_,
118 /*pmax=*/pressureHigh_,
119
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 /*np=*/nPressure_);
120 2 }
121
122 ✗ void setTime( Scalar time )
123 {
124
1/2
✓ Branch 1 taken 85 times.
✗ Branch 2 not taken.
 85 time_ = time;
125 ✗ }
126
127 ✗ void setTimeStepSize( Scalar timeStepSize )
128 {
129
1/2
✓ Branch 1 taken 85 times.
✗ Branch 2 not taken.
 85 timeStepSize_ = timeStepSize;
130 ✗ }
131
132
133 /*!
134 * \brief The problem name.
135 *
136 * This is used as a prefix for files generated by the simulation.
137 */
138 const std::string& name() const
139
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 { return name_; }
140
141 /*!
142 * \brief Specifies which kind of boundary condition should be
143 * used for which equation on a given boundary segment.
144 */
145 155642 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
146 {
147 155642 BoundaryTypes bcTypes;
148
149 466926 const Scalar rmax = this->gridGeometry().bBoxMax()[0];
150 466926 const Scalar rmin = this->gridGeometry().bBoxMin()[0];
151
152 // default to Neumann
153 155642 bcTypes.setAllNeumann();
154
155 // Constant pressure at reservoir boundary (Dirichlet condition)
156
4/4
✓ Branch 0 taken 27970 times.
✓ Branch 1 taken 127672 times.
✓ Branch 2 taken 27970 times.
✓ Branch 3 taken 127672 times.
 311284 if(globalPos[0] > rmax - eps_)
157 bcTypes.setAllDirichlet();
158
159 // Constant pressure at well (Dirichlet condition)
160
4/4
✓ Branch 0 taken 27970 times.
✓ Branch 1 taken 127672 times.
✓ Branch 2 taken 27970 times.
✓ Branch 3 taken 127672 times.
 311284 if(globalPos[0] < rmin + eps_)
161 bcTypes.setAllDirichlet();
162
163 155642 return bcTypes;
164 }
165
166 /*!
167 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
168 */
169 17860 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
170 {
171 17860 PrimaryVariables priVars(0.0);
172
2/2
✓ Branch 0 taken 8930 times.
✓ Branch 1 taken 8930 times.
 17860 priVars.setState(bothPhases);
173
174
6/6
✓ Branch 0 taken 8930 times.
✓ Branch 1 taken 8930 times.
✓ Branch 2 taken 8930 times.
✓ Branch 3 taken 8930 times.
✓ Branch 4 taken 8930 times.
✓ Branch 5 taken 8930 times.
 53580 const Scalar rmax = this->gridGeometry().bBoxMax()[0];
175
6/6
✓ Branch 0 taken 8930 times.
✓ Branch 1 taken 8930 times.
✓ Branch 2 taken 8930 times.
✓ Branch 3 taken 8930 times.
✓ Branch 4 taken 8930 times.
✓ Branch 5 taken 8930 times.
 53580 const Scalar rmin = this->gridGeometry().bBoxMin()[0];
176
177
4/4
✓ Branch 0 taken 8930 times.
✓ Branch 1 taken 8930 times.
✓ Branch 2 taken 8930 times.
✓ Branch 3 taken 8930 times.
 35720 if(globalPos[0] > rmax - eps_)
178 {
179 8930 priVars[pressureIdx] = outerPressure_ ; // Outer boundary pressure bar
180 8930 priVars[switchIdx] = outerLiqSaturation_; // Saturation outer boundary
181 8930 priVars[xwNaClIdx] = massToMoleFrac_(outerSalinity_);// mole fraction salt
182 17860 priVars[precipNaClIdx] = 0.0;// precipitated salt
183 }
184
185
4/4
✓ Branch 0 taken 8930 times.
✓ Branch 1 taken 8930 times.
✓ Branch 2 taken 8930 times.
✓ Branch 3 taken 8930 times.
 35720 if(globalPos[0] < rmin + eps_)
186 {
187 8930 priVars[pressureIdx] = innerPressure_ ; // Inner boundary pressure bar
188 8930 priVars[switchIdx] = innerLiqSaturation_; // Saturation inner boundary
189 8930 priVars[xwNaClIdx] = massToMoleFrac_(innerSalinity_);// mole fraction salt
190 17860 priVars[precipNaClIdx] = 0.0;// precipitated salt
191 }
192
193 17860 return priVars;
194 }
195
196 /*!
197 * \brief Evaluates the initial value for a control volume.
198 *
199 * \param globalPos The global position
200 *
201 */
202 431 PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
203 {
204 431 PrimaryVariables priVars(0.0);
205 431 priVars.setState(bothPhases);
206
207 431 priVars[pressureIdx] = reservoirPressure_;
208 431 priVars[switchIdx] = initLiqSaturation_; // Sw primary variable
209
2/2
✓ Branch 1 taken 378 times.
✓ Branch 2 taken 53 times.
 431 priVars[xwNaClIdx] = massToMoleFrac_(outerSalinity_); // mole fraction
210
8/8
✓ Branch 0 taken 378 times.
✓ Branch 1 taken 53 times.
✓ Branch 2 taken 378 times.
✓ Branch 3 taken 53 times.
✓ Branch 4 taken 157 times.
✓ Branch 5 taken 221 times.
✓ Branch 6 taken 157 times.
✓ Branch 7 taken 221 times.
 862 if(globalPos[0] > 5.0 - eps_ && globalPos[0] < 19.0 + eps_)
211 314 priVars[precipNaClIdx] = initPrecipitatedSaltBlock_; // [kg/m^3]
212 else
213 548 priVars[precipNaClIdx] = 0.0; // [kg/m^3]
214
215 431 return priVars;
216 }
217
218
219 /*!
220 * \brief Evaluates the source term for all phases within a given
221 * sub-controlvolume.
222 *
223 * This is the method for the case where the source term is
224 * potentially solution dependent and requires some quantities that
225 * are specific to the fully-implicit method.
226 *
227 * \param element The finite element
228 * \param fvGeometry The finite-volume geometry
229 * \param elemVolVars All volume variables for the element
230 * \param scv The subcontrolvolume
231 *
232 */
233 3791600 NumEqVector source(const Element &element,
234 const FVElementGeometry& fvGeometry,
235 const ElementVolumeVariables& elemVolVars,
236 const SubControlVolume &scv) const
237 {
238 3791600 NumEqVector source(0.0);
239
240 3791600 const auto& volVars = elemVolVars[scv];
241
242 3791600 Scalar moleFracNaCl_wPhase = volVars.moleFraction(liquidPhaseIdx, NaClIdx);
243 3791600 Scalar massFracNaCl_Max_wPhase = this->spatialParams().solubilityLimit();
244 3791600 Scalar moleFracNaCl_Max_wPhase = massToMoleFrac_(massFracNaCl_Max_wPhase);
245 3791600 Scalar saltPorosity = this->spatialParams().minimalPorosity(element, scv);
246
247 // liquid phase
248 using std::abs;
249 7583200 Scalar precipSalt = volVars.porosity() * volVars.molarDensity(liquidPhaseIdx)
250
2/2
✓ Branch 0 taken 1498298 times.
✓ Branch 1 taken 2293302 times.
 3791600 * volVars.saturation(liquidPhaseIdx)
251
2/2
✓ Branch 0 taken 1498298 times.
✓ Branch 1 taken 2293302 times.
 3791600 * abs(moleFracNaCl_wPhase - moleFracNaCl_Max_wPhase)
252 3791600 / timeStepSize_;
253
2/2
✓ Branch 0 taken 1498298 times.
✓ Branch 1 taken 2293302 times.
 3791600 if (moleFracNaCl_wPhase < moleFracNaCl_Max_wPhase)
254 1498298 precipSalt *= -1;
255
256 // make sure we don't dissolve more salt than previously precipitated
257
2/2
✓ Branch 2 taken 369065 times.
✓ Branch 3 taken 3422535 times.
 7583200 if (precipSalt*timeStepSize_ + volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)< 0)
258 738130 precipSalt = -volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)/timeStepSize_;
259
260
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 3791600 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3791600 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3791600 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 11374800 if (volVars.solidVolumeFraction(sPhaseIdx) >= this->spatialParams().referencePorosity(element, scv) - saltPorosity && precipSalt > 0)
261 precipSalt = 0;
262
263 7583200 source[conti0EqIdx + NaClIdx] += -precipSalt;
264 7583200 source[precipNaClEqIdx] += precipSalt;
265 3791600 return source;
266
267 }
268
269 /*!
270 * \brief Adds additional VTK output data to the VTKWriter.
271 *
272 * Function is called by the output module on every write.
273 */
274
275 const std::vector<Scalar>& getPermeability()
276 {
277
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return permeability_;
278 }
279
280 87 void updateVtkOutput(const SolutionVector& curSol)
281 {
282
2/4
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
 174 auto fvGeometry = localView(this->gridGeometry());
283
5/10
✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 87 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 87 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 17487 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 8800 times.
✗ Branch 14 not taken.
 35148 for (const auto& element : elements(this->gridGeometry().gridView()))
284 {
285
2/4
✓ Branch 1 taken 8800 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8800 times.
✗ Branch 5 not taken.
 34800 const auto elemSol = elementSolution(element, curSol, this->gridGeometry());
286
1/2
✓ Branch 1 taken 17400 times.
✗ Branch 2 not taken.
 17400 fvGeometry.bindElement(element);
287
288
4/4
✓ Branch 0 taken 43800 times.
✓ Branch 1 taken 17400 times.
✓ Branch 2 taken 35200 times.
✓ Branch 3 taken 8800 times.
 105200 for (auto&& scv : scvs(fvGeometry))
289 {
290 43800 VolumeVariables volVars;
291
1/2
✓ Branch 1 taken 43800 times.
✗ Branch 2 not taken.
 43800 volVars.update(elemSol, *this, element, scv);
292 43800 const auto dofIdxGlobal = scv.dofIndex();
293 87600 permeability_[dofIdxGlobal] = volVars.permeability();
294 }
295 }
296 87 }
297
298 private:
299
300 /*!
301 * \brief Returns the molality of NaCl (mol NaCl / kg water) for a given mole fraction.
302 *
303 * \param massFracNaClInWater The mass fraction of NaCl in liquid phase water (kg NaCl / kg solution)
304 */
305 3809891 static Scalar massToMoleFrac_(Scalar massFracNaClInWater)
306 {
307 3809891 const Scalar molarMassWater = FluidSystem::molarMass(H2OIdx); // in kg/mol
308 3809891 const Scalar molarMassNaCl = FluidSystem::molarMass(NaClIdx); // in kg/mol
309
310 3809891 const auto moleFracNaClInWater = -molarMassWater * massFracNaClInWater / ((molarMassNaCl - molarMassWater) * massFracNaClInWater - molarMassNaCl);
311 3809891 return moleFracNaClInWater;
312 }
313
314 int nTemperature_;
315 int nPressure_;
316 std::string name_;
317
318 Scalar pressureLow_, pressureHigh_;
319 Scalar temperatureLow_, temperatureHigh_;
320 Scalar outerSalinity_;
321 Scalar reservoirPressure_;
322 Scalar innerPressure_;
323 Scalar outerPressure_;
324 Scalar initLiqSaturation_;
325 Scalar outerLiqSaturation_;
326 Scalar innerLiqSaturation_;
327 Scalar initPrecipitatedSaltBlock_;
328 Scalar innerSalinity_;
329 Scalar time_ = 0.0;
330 Scalar timeStepSize_ = 0.0;
331 static constexpr Scalar eps_ = 1e-6;
332 Scalar reservoirSaturation_;
333 std::vector<double> permeability_;
334 };
335
336 } // end namespace Dumux
337
338 #endif
339