GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/2pncmin/isothermal/problem.hh
Date: 2024-05-04 19:09:25
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 3 DissolutionProblem(std::shared_ptr<const GridGeometry> gridGeometry)
89
8/24
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 3 times.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 3 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.
 9 : ParentType(gridGeometry)
90 {
91
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 outerSalinity_ = getParam<Scalar>("Problem.OuterSalinity");
92
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 reservoirPressure_ = getParam<Scalar>("Problem.ReservoirPressure");
93
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 initLiqSaturation_ = getParam<Scalar>("Problem.LiquidSaturation");
94
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 initPrecipitatedSaltBlock_ = getParam<Scalar>("Problem.InitPrecipitatedSaltBlock");
95
96
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 outerLiqSaturation_ = getParam<Scalar>("Problem.OuterLiqSaturation");
97
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 innerLiqSaturation_ = getParam<Scalar>("Problem.InnerLiqSaturation");
98
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 innerSalinity_ = getParam<Scalar>("Problem.InnerSalinity");
99
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 innerPressure_ = getParam<Scalar>("Problem.InnerPressure");
100
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 outerPressure_ = getParam<Scalar>("Problem.OuterPressure");
101
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 reservoirSaturation_ = getParam<Scalar>("Problem.reservoirSaturation");
102
103
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 nTemperature_ = getParam<int>("FluidSystem.NTemperature");
104
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 nPressure_ = getParam<int>("FluidSystem.NPressure");
105
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 pressureLow_ = getParam<Scalar>("FluidSystem.PressureLow");
106
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 pressureHigh_ = getParam<Scalar>("FluidSystem.PressureHigh");
107
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 temperatureLow_ = getParam<Scalar>("FluidSystem.TemperatureLow");
108
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 temperatureHigh_ = getParam<Scalar>("FluidSystem.TemperatureHigh");
109
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
 3 name_ = getParam<std::string>("Problem.Name");
110
111 3 unsigned int codim = GetPropType<TypeTag, Properties::GridGeometry>::discMethod == DiscretizationMethods::box ? dim : 0;
112
4/8
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
 9 permeability_.resize(gridGeometry->gridView().size(codim));
113
114 6 FluidSystem::init(/*Tmin=*/temperatureLow_,
115 /*Tmax=*/temperatureHigh_,
116 3 /*nT=*/nTemperature_,
117 /*pmin=*/pressureLow_,
118 /*pmax=*/pressureHigh_,
119
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 /*np=*/nPressure_);
120 3 }
121
122 ✗ void setTime( Scalar time )
123 {
124
1/2
✓ Branch 1 taken 90 times.
✗ Branch 2 not taken.
 90 time_ = time;
125 ✗ }
126
127 ✗ void setTimeStepSize( Scalar timeStepSize )
128 {
129
1/2
✓ Branch 1 taken 90 times.
✗ Branch 2 not taken.
 90 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 3 times.
✗ Branch 2 not taken.
 3 { 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 162818 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
146 {
147 162818 BoundaryTypes bcTypes;
148
149 488454 const Scalar rmax = this->gridGeometry().bBoxMax()[0];
150 488454 const Scalar rmin = this->gridGeometry().bBoxMin()[0];
151
152 // default to Neumann
153 162818 bcTypes.setAllNeumann();
154
155 // Constant pressure at reservoir boundary (Dirichlet condition)
156
4/4
✓ Branch 0 taken 29430 times.
✓ Branch 1 taken 133388 times.
✓ Branch 2 taken 29430 times.
✓ Branch 3 taken 133388 times.
 325636 if(globalPos[0] > rmax - eps_)
157 bcTypes.setAllDirichlet();
158
159 // Constant pressure at well (Dirichlet condition)
160
4/4
✓ Branch 0 taken 29430 times.
✓ Branch 1 taken 133388 times.
✓ Branch 2 taken 29430 times.
✓ Branch 3 taken 133388 times.
 325636 if(globalPos[0] < rmin + eps_)
161 bcTypes.setAllDirichlet();
162
163 162818 return bcTypes;
164 }
165
166 /*!
167 * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
168 */
169 19420 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
170 {
171 19420 PrimaryVariables priVars(0.0);
172
2/2
✓ Branch 0 taken 9710 times.
✓ Branch 1 taken 9710 times.
 19420 priVars.setState(bothPhases);
173
174
6/6
✓ Branch 0 taken 9710 times.
✓ Branch 1 taken 9710 times.
✓ Branch 2 taken 9710 times.
✓ Branch 3 taken 9710 times.
✓ Branch 4 taken 9710 times.
✓ Branch 5 taken 9710 times.
 58260 const Scalar rmax = this->gridGeometry().bBoxMax()[0];
175
6/6
✓ Branch 0 taken 9710 times.
✓ Branch 1 taken 9710 times.
✓ Branch 2 taken 9710 times.
✓ Branch 3 taken 9710 times.
✓ Branch 4 taken 9710 times.
✓ Branch 5 taken 9710 times.
 58260 const Scalar rmin = this->gridGeometry().bBoxMin()[0];
176
177
4/4
✓ Branch 0 taken 9710 times.
✓ Branch 1 taken 9710 times.
✓ Branch 2 taken 9710 times.
✓ Branch 3 taken 9710 times.
 38840 if(globalPos[0] > rmax - eps_)
178 {
179 9710 priVars[pressureIdx] = outerPressure_ ; // Outer boundary pressure bar
180 9710 priVars[switchIdx] = outerLiqSaturation_; // Saturation outer boundary
181 9710 priVars[xwNaClIdx] = massToMoleFrac_(outerSalinity_);// mole fraction salt
182 19420 priVars[precipNaClIdx] = 0.0;// precipitated salt
183 }
184
185
4/4
✓ Branch 0 taken 9710 times.
✓ Branch 1 taken 9710 times.
✓ Branch 2 taken 9710 times.
✓ Branch 3 taken 9710 times.
 38840 if(globalPos[0] < rmin + eps_)
186 {
187 9710 priVars[pressureIdx] = innerPressure_ ; // Inner boundary pressure bar
188 9710 priVars[switchIdx] = innerLiqSaturation_; // Saturation inner boundary
189 9710 priVars[xwNaClIdx] = massToMoleFrac_(innerSalinity_);// mole fraction salt
190 19420 priVars[precipNaClIdx] = 0.0;// precipitated salt
191 }
192
193 19420 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 4254000 NumEqVector source(const Element &element,
234 const FVElementGeometry& fvGeometry,
235 const ElementVolumeVariables& elemVolVars,
236 const SubControlVolume &scv) const
237 {
238 4254000 NumEqVector source(0.0);
239
240 4254000 const auto& volVars = elemVolVars[scv];
241
242 4254000 Scalar moleFracNaCl_wPhase = volVars.moleFraction(liquidPhaseIdx, NaClIdx);
243 4254000 Scalar massFracNaCl_Max_wPhase = this->spatialParams().solubilityLimit();
244 4254000 Scalar moleFracNaCl_Max_wPhase = massToMoleFrac_(massFracNaCl_Max_wPhase);
245 4254000 Scalar saltPorosity = this->spatialParams().minimalPorosity(element, scv);
246
247 // liquid phase
248 using std::abs;
249 8508000 Scalar precipSalt = volVars.porosity() * volVars.molarDensity(liquidPhaseIdx)
250
2/2
✓ Branch 0 taken 1879302 times.
✓ Branch 1 taken 2374698 times.
 4254000 * volVars.saturation(liquidPhaseIdx)
251
2/2
✓ Branch 0 taken 1879302 times.
✓ Branch 1 taken 2374698 times.
 4254000 * abs(moleFracNaCl_wPhase - moleFracNaCl_Max_wPhase)
252 4254000 / timeStepSize_;
253
2/2
✓ Branch 0 taken 1879302 times.
✓ Branch 1 taken 2374698 times.
 4254000 if (moleFracNaCl_wPhase < moleFracNaCl_Max_wPhase)
254 1879302 precipSalt *= -1;
255
256 // make sure we don't dissolve more salt than previously precipitated
257
2/2
✓ Branch 2 taken 441621 times.
✓ Branch 3 taken 3812379 times.
 8508000 if (precipSalt*timeStepSize_ + volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)< 0)
258 883242 precipSalt = -volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)/timeStepSize_;
259
260
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 4254000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4254000 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4254000 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 12762000 if (volVars.solidVolumeFraction(sPhaseIdx) >= this->spatialParams().referencePorosity(element, scv) - saltPorosity && precipSalt > 0)
261 precipSalt = 0;
262
263 8508000 source[conti0EqIdx + NaClIdx] += -precipSalt;
264 8508000 source[precipNaClEqIdx] += precipSalt;
265 4254000 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 3 times.
✗ Branch 2 not taken.
 3 return permeability_;
278 }
279
280 93 void updateVtkOutput(const SolutionVector& curSol)
281 {
282
2/4
✓ Branch 1 taken 93 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 93 times.
✗ Branch 5 not taken.
 186 auto fvGeometry = localView(this->gridGeometry());
283
5/10
✓ Branch 1 taken 93 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 93 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 93 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 18693 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 10000 times.
✗ Branch 14 not taken.
 37572 for (const auto& element : elements(this->gridGeometry().gridView()))
284 {
285
2/4
✓ Branch 1 taken 10000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10000 times.
✗ Branch 5 not taken.
 37200 const auto elemSol = elementSolution(element, curSol, this->gridGeometry());
286
1/2
✓ Branch 1 taken 18600 times.
✗ Branch 2 not taken.
 18600 fvGeometry.bindElement(element);
287
288
4/4
✓ Branch 0 taken 48600 times.
✓ Branch 1 taken 18600 times.
✓ Branch 2 taken 40000 times.
✓ Branch 3 taken 10000 times.
 117200 for (auto&& scv : scvs(fvGeometry))
289 {
290 48600 VolumeVariables volVars;
291
1/2
✓ Branch 1 taken 48600 times.
✗ Branch 2 not taken.
 48600 volVars.update(elemSol, *this, element, scv);
292 48600 const auto dofIdxGlobal = scv.dofIndex();
293 97200 permeability_[dofIdxGlobal] = volVars.permeability();
294 }
295 }
296 93 }
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 4273851 static Scalar massToMoleFrac_(Scalar massFracNaClInWater)
306 {
307 4273851 const Scalar molarMassWater = FluidSystem::molarMass(H2OIdx); // in kg/mol
308 4273851 const Scalar molarMassNaCl = FluidSystem::molarMass(NaClIdx); // in kg/mol
309
310 4273851 const auto moleFracNaClInWater = -molarMassWater * massFracNaClInWater / ((molarMassNaCl - molarMassWater) * massFracNaClInWater - molarMassNaCl);
311 4273851 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