GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/constraintsolvers/immiscibleflash.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 100 126 79.4%
Functions: 5 12 41.7%
Branches: 61 102 59.8%
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 ConstraintSolvers
10 * \brief Determines the pressures and saturations of all fluid phases
11 * given the total mass of all components.
12 */
13 #ifndef DUMUX_IMMISCIBLE_FLASH_HH
14 #define DUMUX_IMMISCIBLE_FLASH_HH
15
16 #include <dune/common/fvector.hh>
17 #include <dune/common/fmatrix.hh>
18
19 #include <dumux/common/exceptions.hh>
20 #include <dumux/material/fluidmatrixinteractions/mp/mpadapter.hh>
21
22 namespace Dumux {
23
24 /*!
25 * \ingroup ConstraintSolvers
26 * \brief Determines the pressures and saturations of all fluid phases
27 * given the total mass of all components.
28 *
29 * In a N-phase, N-component context, we have the following
30 * unknowns if assuming immiscibility:
31 *
32 * - N pressures
33 * - N saturations
34 *
35 * This sums up to 2*N unknowns. On the equations side of things, we
36 * have:
37 *
38 * - N total component molarities
39 * - 1 The sum of all saturations is 1
40 * - N-1 Relations from capillary pressure
41 *
42 * this also sums up to 2*N. We include the capillary pressures and
43 * the sum of the saturations explicitly. This means that we only
44 * solve for the first pressure and N-1 saturations.
45 *
46 * If a fluid phase is incompressible, the pressure cannot determined
47 * by this, though. In this case the original pressure is kept, and
48 * the saturation of the phase is calculated by dividing the global
49 * molarity of the component by the phase density.
50 */
51 template <class Scalar, class FluidSystem>
52 class ImmiscibleFlash
53 {
54 static constexpr int numPhases = FluidSystem::numPhases;
55 static constexpr int numComponents = FluidSystem::numComponents;
56 static_assert(numPhases == numComponents,
57 "Immiscibility assumes that the number of phases"
58 " is equal to the number of components");
59
60 using ParameterCache = typename FluidSystem::ParameterCache;
61
62 static constexpr int numEq = numPhases;
63
64 using Matrix = Dune::FieldMatrix<Scalar, numEq, numEq>;
65 using Vector = Dune::FieldVector<Scalar, numEq>;
66
67 public:
68 using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
69
70 /*!
71 * \brief Construct a new flash
72 * \param wettingPhaseIdx the phase index of the wetting phase
73 */
74 4 explicit ImmiscibleFlash(int wettingPhaseIdx = 0)
75 4 : wettingPhaseIdx_(wettingPhaseIdx)
76 {}
77
78 /*!
79 * \brief Guess initial values for all quantities.
80 * \param fluidState Thermodynamic state of the fluids
81 * \param paramCache Container for cache parameters
82 * \param globalMolarities
83 */
84 template <class FluidState>
85 ✗ void guessInitial(FluidState &fluidState,
86 ParameterCache &paramCache,
87 const ComponentVector &globalMolarities)
88 {
89 // the sum of all molarities
90 ✗ Scalar sumMoles = 0;
91 ✗ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
92 sumMoles += globalMolarities[compIdx];
93
94
2/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 4 times.
 12 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
95 // pressure. use atmospheric pressure as initial guess
96 8 fluidState.setPressure(phaseIdx, 2e5);
97
98 // saturation. assume all fluids to be equally distributed
99 16 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
100 }
101 ✗ }
102
103 /*!
104 * \brief Calculates the chemical equilibrium from the component
105 * fugacities in a phase.
106 * \param fluidState Thermodynamic state of the fluids
107 * \param paramCache Container for cache parameters
108 * \param globalMolarities
109 * \param material The material law object
110 *
111 * The phase's fugacities must already be set.
112 */
113 template <class MaterialLaw, class FluidState>
114 4 void solve(FluidState& fluidState,
115 ParameterCache& paramCache,
116 const MaterialLaw& material,
117 const ComponentVector& globalMolarities)
118 {
119 /////////////////////////
120 // Check if all fluid phases are incompressible
121 /////////////////////////
122 4 bool allIncompressible = true;
123 4 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
124 4 if (FluidSystem::isCompressible(phaseIdx)) {
125 4 allIncompressible = false;
126 break;
127 }
128 }
129
130 if (allIncompressible) {
131 paramCache.updateAll(fluidState);
132 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
133 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
134 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
135 fluidState.setDensity(phaseIdx, rho);
136 fluidState.setMolarDensity(phaseIdx, rhoMolar);
137
138 Scalar saturation =
139 globalMolarities[/*compIdx=*/phaseIdx]
140 / fluidState.molarDensity(phaseIdx);
141 fluidState.setSaturation(phaseIdx, saturation);
142 }
143 }
144
145 // /////////////////////////
146 // Newton method
147 /////////////////////////
148
149 // Jacobian matrix
150 4 Matrix J;
151 // solution, i.e. phase composition
152 4 Vector deltaX;
153 // right hand side
154 4 Vector b;
155
156 4 completeFluidState_(fluidState, paramCache, material);
157
158 const int nMax = 50; // <- maximum number of newton iterations
159
1/2
✓ Branch 0 taken 42 times.
✗ Branch 1 not taken.
 42 for (int nIdx = 0; nIdx < nMax; ++nIdx) {
160 // calculate Jacobian matrix and right hand side
161 42 linearize_(J, b, fluidState, paramCache, material, globalMolarities);
162
163 // Solve J*x = b
164 42 deltaX = 0;
165
166 42 try { J.solve(deltaX, b); }
167 catch (Dune::FMatrixError& e)
168 {
169 /*
170 std::cout << "error: " << e << "\n";
171 std::cout << "b: " << b << "\n";
172 */
173
174 throw NumericalProblem(e.what());
175 }
176
177 /*
178 printFluidState_(fluidState);
179 std::cout << "J:\n";
180 for (int i = 0; i < numEq; ++i) {
181 for (int j = 0; j < numEq; ++j) {
182 std::ostringstream os;
183 os << J[i][j];
184
185 std::string s(os.str());
186 do {
187 s += " ";
188 } while (s.size() < 20);
189 std::cout << s;
190 }
191 std::cout << "\n";
192 }
193 std::cout << "residual: " << b << "\n";
194 std::cout << "deltaX: " << deltaX << "\n";
195 std::cout << "---------------\n";
196 */
197
198 // update the fluid quantities.
199 42 Scalar relError = update_(fluidState, paramCache, material, deltaX);
200
201
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 38 times.
 42 if (relError < 1e-9)
202 4 return;
203 }
204
205 /*
206 printFluidState_(fluidState);
207 std::cout << "globalMolarities: ";
208 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
209 std::cout << globalMolarities[compIdx] << " ";
210 std::cout << "\n";
211 */
212
213 ✗ DUNE_THROW(NumericalProblem,
214 "Flash calculation failed."
215 " {c_alpha^kappa} = {" << globalMolarities << "}, T = "
216 << fluidState.temperature(/*phaseIdx=*/0));
217 }
218
219
220 protected:
221 template <class FluidState>
222 void printFluidState_(const FluidState &fs)
223 {
224 std::cout << "saturations: ";
225 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
226 std::cout << fs.saturation(phaseIdx) << " ";
227 std::cout << "\n";
228
229 std::cout << "pressures: ";
230 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
231 std::cout << fs.pressure(phaseIdx) << " ";
232 std::cout << "\n";
233
234 std::cout << "densities: ";
235 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
236 std::cout << fs.density(phaseIdx) << " ";
237 std::cout << "\n";
238
239 std::cout << "global component molarities: ";
240 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
241 Scalar sum = 0;
242 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
243 sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
244 }
245 std::cout << sum << " ";
246 }
247 std::cout << "\n";
248 }
249
250 template <class MaterialLaw, class FluidState>
251 42 void linearize_(Matrix &J,
252 Vector &b,
253 FluidState &fluidState,
254 ParameterCache &paramCache,
255 const MaterialLaw& material,
256 const ComponentVector &globalMolarities)
257 {
258 42 FluidState origFluidState(fluidState);
259 ParameterCache origParamCache(paramCache);
260
261 42 Vector tmp;
262
263 // reset jacobian
264 42 J = 0;
265
266 calculateDefect_(b, fluidState, fluidState, globalMolarities);
267
268 // assemble jacobian matrix
269
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 42 times.
 126 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
270 ////////
271 // approximately calculate partial derivatives of the
272 // fugacity defect of all components in regard to the mole
273 // fraction of the i-th component. This is done via
274 // forward differences
275
276 // deviate the mole fraction of the i-th component
277 168 Scalar x_i = getQuantity_(fluidState, pvIdx);
278 168 const Scalar eps = 1e-10/quantityWeight_(fluidState, pvIdx);
279 84 setQuantity_<MaterialLaw>(fluidState, paramCache, material, pvIdx, x_i + eps);
280
281 // compute derivative of the defect
282 calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
283 tmp -= b;
284 tmp /= eps;
285 // store derivative in jacobian matrix
286
2/2
✓ Branch 0 taken 168 times.
✓ Branch 1 taken 84 times.
 252 for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
287 672 J[eqIdx][pvIdx] = tmp[eqIdx];
288
289 // fluid state and parameter cache to their original values
290 84 fluidState = origFluidState;
291 paramCache = origParamCache;
292
293 // end forward differences
294 ////////
295 }
296 42 }
297
298 template <class FluidState>
299 ✗ void calculateDefect_(Vector &b,
300 const FluidState &fluidStateEval,
301 const FluidState &fluidState,
302 const ComponentVector &globalMolarities)
303 {
304 // global molarities are given
305
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 84 times.
✓ Branch 4 taken 84 times.
✓ Branch 5 taken 168 times.
 378 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
306 504 b[phaseIdx] =
307 fluidState.saturation(phaseIdx)
308 756 * fluidState.molarity(phaseIdx, /*compIdx=*/phaseIdx);
309
310 756 b[phaseIdx] -= globalMolarities[/*compIdx=*/phaseIdx];
311 }
312 ✗ }
313
314 template <class MaterialLaw, class FluidState>
315 42 Scalar update_(FluidState &fluidState,
316 ParameterCache &paramCache,
317 const MaterialLaw& material,
318 const Vector &deltaX)
319 {
320 42 Scalar relError = 0;
321
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 42 times.
 126 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
322 168 Scalar tmp = getQuantity_(fluidState, pvIdx);
323
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 42 times.
 84 Scalar delta = deltaX[pvIdx];
324 using std::max;
325 using std::abs;
326 using std::min;
327
4/4
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 45 times.
✓ Branch 3 taken 39 times.
 126 relError = max(relError, abs(delta)*quantityWeight_(fluidState, pvIdx));
328
329
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 42 times.
 84 if (isSaturationIdx_(pvIdx)) {
330 // dampen to at most 20% change in saturation per
331 // iteration
332
4/4
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 40 times.
✓ Branch 3 taken 2 times.
 122 delta = min(0.2, max(-0.2, delta));
333 }
334 42 else if (isPressureIdx_(pvIdx)) {
335 // dampen to at most 30% change in pressure per
336 // iteration
337
5/6
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 24 times.
✓ Branch 4 taken 42 times.
✗ Branch 5 not taken.
 144 delta = min(0.30*fluidState.pressure(0), max(-0.30*fluidState.pressure(0), delta));
338 }
339
340 84 setQuantityRaw_(fluidState, pvIdx, tmp - delta);
341 }
342
343 42 completeFluidState_<MaterialLaw>(fluidState, paramCache, material);
344
345 42 return relError;
346 }
347
348 template <class MaterialLaw, class FluidState>
349 46 void completeFluidState_(FluidState &fluidState,
350 ParameterCache &paramCache,
351 const MaterialLaw& material)
352 {
353 // calculate the saturation of the last phase as a function of
354 // the other saturations
355 46 Scalar sumSat = 0.0;
356
2/2
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 46 times.
 92 for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
357 92 sumSat += fluidState.saturation(phaseIdx);
358
359
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 45 times.
 46 if (sumSat > 1.0) {
360 // make sure that the last saturation does not become
361 // negative
362
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
363 {
364 1 Scalar S = fluidState.saturation(phaseIdx);
365 2 fluidState.setSaturation(phaseIdx, S/sumSat);
366 }
367 sumSat = 1;
368 }
369
370 // set the last saturation
371 46 fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
372
373 // update the pressures using the material law (saturations
374 // and first pressure are already set because it is implicitly
375 // solved for.)
376 46 const auto pc = material.capillaryPressures(fluidState, wettingPhaseIdx_);
377
2/2
✓ Branch 0 taken 46 times.
✓ Branch 1 taken 46 times.
 92 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
378 46 fluidState.setPressure(phaseIdx,
379 fluidState.pressure(0)
380 184 + (pc[phaseIdx] - pc[0]));
381
382 // update the parameter cache
383 paramCache.updateAll(fluidState, /*except=*/ParameterCache::Temperature|ParameterCache::Composition);
384
385 // update all densities
386
2/2
✓ Branch 0 taken 92 times.
✓ Branch 1 taken 46 times.
 138 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
387 92 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
388 92 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
389 92 fluidState.setDensity(phaseIdx, rho);
390 184 fluidState.setMolarDensity(phaseIdx, rhoMolar);
391 }
392 46 }
393
394 ✗ bool isPressureIdx_(int pvIdx)
395 ✗ { return pvIdx == 0; }
396
397 ✗ bool isSaturationIdx_(int pvIdx)
398 ✗ { return 1 <= pvIdx && pvIdx < numPhases; }
399
400 // retrieves a quantity from the fluid state
401 template <class FluidState>
402 ✗ Scalar getQuantity_(const FluidState &fs, int pvIdx)
403 {
404 ✗ assert(pvIdx < numEq);
405
406 // first pressure
407
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 42 times.
✓ Branch 4 taken 42 times.
✓ Branch 5 taken 42 times.
 168 if (pvIdx < 1) {
408 84 int phaseIdx = 0;
409 168 return fs.pressure(phaseIdx);
410 }
411 // saturations
412 else {
413 assert(pvIdx < numPhases);
414 84 int phaseIdx = pvIdx - 1;
415 168 return fs.saturation(phaseIdx);
416 }
417 }
418
419 // set a quantity in the fluid state
420 template <class MaterialLaw, class FluidState>
421 84 void setQuantity_(FluidState &fs,
422 ParameterCache &paramCache,
423 const MaterialLaw& material,
424 int pvIdx,
425 Scalar value)
426 {
427
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 84 times.
 84 assert(0 <= pvIdx && pvIdx < numEq);
428
429
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 42 times.
 84 if (pvIdx < 1) {
430 // -> first pressure
431 42 Scalar delta = value - fs.pressure(0);
432
433 // set all pressures. here we assume that the capillary
434 // pressure does not depend on absolute pressure.
435
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 84 times.
 126 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
436 252 fs.setPressure(phaseIdx, fs.pressure(phaseIdx) + delta);
437 paramCache.updateAllPressures(fs);
438
439 // update all densities
440
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 42 times.
 126 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
441 84 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
442 84 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
443 84 fs.setDensity(phaseIdx, rho);
444 168 fs.setMolarDensity(phaseIdx, rhoMolar);
445 }
446 }
447 else {
448 assert(pvIdx < numPhases);
449
450 // -> first M-1 saturations
451 42 Scalar delta = value - fs.saturation(/*phaseIdx=*/pvIdx - 1);
452 42 fs.setSaturation(/*phaseIdx=*/pvIdx - 1, value);
453
454 // set last saturation (-> minus the change of the
455 // saturation of the other phases)
456 84 fs.setSaturation(/*phaseIdx=*/numPhases - 1,
457 84 fs.saturation(numPhases - 1) - delta);
458
459 // update all fluid pressures using the capillary pressure
460 // law
461 42 const auto pc = material.capillaryPressures(fs, wettingPhaseIdx_);
462
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 42 times.
 84 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
463 42 fs.setPressure(phaseIdx,
464 fs.pressure(0)
465 168 + (pc[phaseIdx] - pc[0]));
466 paramCache.updateAllPressures(fs);
467
468 // update all densities
469
2/2
✓ Branch 0 taken 84 times.
✓ Branch 1 taken 42 times.
 126 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
470 84 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
471 84 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
472 84 fs.setDensity(phaseIdx, rho);
473 168 fs.setMolarDensity(phaseIdx, rhoMolar);
474
475 }
476 }
477 84 }
478
479 // set a quantity in the fluid state
480 template <class FluidState>
481 ✗ void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
482 {
483 ✗ assert(pvIdx < numEq);
484
485 // first pressure
486 ✗ if (pvIdx < 1) {
487 ✗ int phaseIdx = 0;
488 ✗ fs.setPressure(phaseIdx, value);
489 }
490 // saturations
491 else {
492 assert(pvIdx < numPhases);
493 ✗ int phaseIdx = pvIdx - 1;
494
495 // make sure that the first M-1 saturations does not get
496 // negative
497 using std::max;
498 ✗ value = max(0.0, value);
499 ✗ fs.setSaturation(phaseIdx, value);
500 }
501 ✗ }
502
503 template <class FluidState>
504 ✗ Scalar quantityWeight_(const FluidState &fs, int pvIdx)
505 {
506 // first pressure
507
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 42 times.
✓ Branch 3 taken 42 times.
✓ Branch 4 taken 42 times.
✓ Branch 5 taken 42 times.
 168 if (pvIdx < 1)
508 ✗ return 1.0/1e5;
509 // first M - 1 saturations
510 else {
511 ✗ assert(pvIdx < numPhases);
512 ✗ return 1.0;
513 }
514 }
515
516 private:
517 int wettingPhaseIdx_ = 0; //!< the phase index of the wetting phase
518
519 };
520
521 } // end namespace Dumux
522
523 #endif
524