GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/constraintsolvers/ncpflash.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 186 210 88.6%
Functions: 16 26 61.5%
Branches: 140 190 73.7%
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 phase compositions, pressures and saturations
11 * given the total mass of all components.
12 */
13 #ifndef DUMUX_NCP_FLASH_HH
14 #define DUMUX_NCP_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 phase compositions, pressures and saturations
27 * given the total mass of all components.
28 *
29 * In a M-phase, N-component context, we have the following
30 * unknowns:
31 *
32 * - M pressures
33 * - M saturations
34 * - M*N mole fractions
35 *
36 * This sums up to M*(N + 2). On the equations side of things,
37 * we have:
38 *
39 * - (M - 1)*N equation stemming from the fact that the
40 * fugacity of any component is the same in all phases
41 * - 1 equation from the closure condition of all saturations
42 * (they sum up to 1)
43 * - M - 1 constraints from the capillary pressures
44 * \f$(-> p_\beta = p_\alpha + p_c\alpha,\beta)\f$
45 * - N constraints from the fact that the total mass of each
46 * component is given \f$(-> sum_\alpha rhoMolar_\alpha *
47 * x_\alpha^\kappa = const)\f$
48 * - M model constraints. Here we use the NCP constraints
49 * (-> 0 = min \f$ {S_\alpha, 1 - \sum_\kappa x_\alpha^\kappa}\f$)
50 *
51 * this also sums up to M*(N + 2).
52 *
53 * We use the following catches: Capillary pressures are taken
54 * into account explicitly, so that only the pressure of the first
55 * phase is solved implicitly, also the closure condition for the
56 * saturations is taken into account explicitly, which means, that
57 * we don't need to implicitly solve for the last
58 * saturation. These two measures reduce the number of unknowns to
59 * M*(N + 1), namely:
60 *
61 * - 1 pressure
62 * - M - 1 saturations
63 * - M*N mole fractions
64 */
65 template <class Scalar, class FluidSystem>
66 class NcpFlash
67 {
68 enum { numPhases = FluidSystem::numPhases };
69 enum { numComponents = FluidSystem::numComponents };
70
71 using ParameterCache = typename FluidSystem::ParameterCache;
72
73 static constexpr int numEq = numPhases*(numComponents + 1);
74
75 using Matrix = Dune::FieldMatrix<Scalar, numEq, numEq>;
76 using Vector = Dune::FieldVector<Scalar, numEq>;
77
78 public:
79 using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
80
81 /*!
82 * \brief Construct a new flash
83 * \param wettingPhaseIdx the phase index of the wetting phase
84 */
85 56 explicit NcpFlash(int wettingPhaseIdx = 0)
86 56 : wettingPhaseIdx_(wettingPhaseIdx)
87 {}
88
89 /*!
90 * \brief Guess initial values for all quantities.
91 * \param fluidState Thermodynamic state of the fluids
92 * \param paramCache Container for cache parameters
93 * \param globalMolarities
94 */
95 template <class FluidState>
96 5 void guessInitial(FluidState &fluidState,
97 ParameterCache &paramCache,
98 const ComponentVector &globalMolarities)
99 {
100 // the sum of all molarities
101 5 Scalar sumMoles = 0;
102
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 15 times.
 20 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
103 30 sumMoles += globalMolarities[compIdx];
104
105
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 5 times.
 16 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
106 // composition
107
2/2
✓ Branch 0 taken 37 times.
✓ Branch 1 taken 11 times.
 48 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
108 37 fluidState.setMoleFraction(phaseIdx,
109 compIdx,
110 74 globalMolarities[compIdx]/sumMoles);
111
112 // pressure. use atmospheric pressure as initial guess
113 11 fluidState.setPressure(phaseIdx, 1.0135e5);
114
115 // saturation. assume all fluids to be equally distributed
116 22 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
117 }
118
119 // set the fugacity coefficients of all components in all phases
120 paramCache.updateAll(fluidState);
121
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 5 times.
 16 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
122
2/2
✓ Branch 0 taken 37 times.
✓ Branch 1 taken 11 times.
 48 for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
123 37 Scalar phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
124 74 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
125 }
126 }
127 5 }
128
129 /*!
130 * \brief Calculates the chemical equilibrium from the component
131 * fugacities in a phase.
132 * \param fluidState Thermodynamic state of the fluids
133 * \param paramCache Container for cache parameters
134 * \param globalMolarities
135 * \param material The material law object
136 *
137 * The phase's fugacities must already be set.
138 */
139 template <class MaterialLaw, class FluidState>
140 664 void solve(FluidState &fluidState,
141 ParameterCache &paramCache,
142 const MaterialLaw& material,
143 const ComponentVector &globalMolarities)
144 {
145 /////////////////////////
146 // Newton method
147 /////////////////////////
148
149 // Jacobian matrix
150 664 Matrix J;
151 // solution, i.e. phase composition
152 664 Vector deltaX;
153 // right hand side
154 664 Vector b;
155
156 // make the fluid state consistent with the fluid system.
157 664 completeFluidState_(fluidState, paramCache, material);
158
159 /*
160 std::cout << "--------------------\n";
161 std::cout << "globalMolarities: ";
162 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
163 std::cout << globalMolarities[compIdx] << " ";
164 std::cout << "\n";
165 */
166
167 const int nMax = 50; // <- maximum number of newton iterations
168
1/2
✓ Branch 0 taken 1744 times.
✗ Branch 1 not taken.
 1744 for (int nIdx = 0; nIdx < nMax; ++nIdx) {
169 // calculate Jacobian matrix and right hand side
170 1744 linearize_(J, b, fluidState, paramCache, material, globalMolarities);
171
172 // Solve J*x = b
173 deltaX = 0;
174
175 try {
176
177 // simple preconditioning of the linear system to reduce condition number
178 int eqIdx = 0;
179
2/2
✓ Branch 0 taken 11798 times.
✓ Branch 1 taken 1744 times.
 13542 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
180 {
181
2/2
✓ Branch 0 taken 23432 times.
✓ Branch 1 taken 11798 times.
 35230 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
182 {
183 46864 J[eqIdx] *= 10e-7; // roughly the magnitude of the fugacities
184 23432 b[eqIdx] *= 10e-7;
185 23432 ++eqIdx;
186 }
187 }
188
189
1/2
✓ Branch 1 taken 1744 times.
✗ Branch 2 not taken.
 1744 J.solve(deltaX, b);
190
191 }
192 ✗ catch (Dune::FMatrixError &e)
193 {
194 /*
195 printFluidState_(fluidState);
196 std::cout << "error: " << e << "\n";
197 std::cout << "b: " << b << "\n";
198 std::cout << "J: " << J << "\n";
199 */
200
201 ✗ throw NumericalProblem(e.what());
202 }
203
204 /*
205 printFluidState_(fluidState);
206 std::cout << "J:\n";
207 for (int i = 0; i < numEq; ++i) {
208 for (int j = 0; j < numEq; ++j) {
209 std::ostringstream os;
210 os << J[i][j];
211
212 std::string s(os.str());
213 do {
214 s += " ";
215 } while (s.size() < 20);
216 std::cout << s;
217 }
218 std::cout << "\n";
219 }
220
221 std::cout << "deltaX: " << deltaX << "\n";
222 std::cout << "---------------\n";
223 */
224
225 // update the fluid quantities.
226 1744 Scalar relError = update_(fluidState, paramCache, material, deltaX);
227
228
2/2
✓ Branch 0 taken 664 times.
✓ Branch 1 taken 1080 times.
 1744 if (relError < 1e-9)
229 664 return;
230 }
231
232 /*
233 printFluidState_(fluidState);
234 std::cout << "globalMolarities: ";
235 for (int compIdx = 0; compIdx < numComponents; ++ compIdx)
236 std::cout << globalMolarities[compIdx] << " ";
237 std::cout << "\n";
238 */
239
240 ✗ DUNE_THROW(NumericalProblem,
241 "Flash calculation failed."
242 " {c_alpha^kappa} = {" << globalMolarities << "}, T = "
243 << fluidState.temperature(/*phaseIdx=*/0));
244 }
245
246
247 protected:
248 template <class FluidState>
249 void printFluidState_(const FluidState &fs)
250 {
251 std::cout << "saturations: ";
252 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
253 std::cout << fs.saturation(phaseIdx) << " ";
254 std::cout << "\n";
255
256 std::cout << "pressures: ";
257 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
258 std::cout << fs.pressure(phaseIdx) << " ";
259 std::cout << "\n";
260
261 std::cout << "densities: ";
262 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
263 std::cout << fs.density(phaseIdx) << " ";
264 std::cout << "\n";
265
266 std::cout << "molar densities: ";
267 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
268 std::cout << fs.molarDensity(phaseIdx) << " ";
269 std::cout << "\n";
270
271 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
272 std::cout << "composition " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
273 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
274 std::cout << fs.moleFraction(phaseIdx, compIdx) << " ";
275 }
276 std::cout << "\n";
277 }
278
279 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
280 std::cout << "fugacities " << FluidSystem::phaseName(phaseIdx) << "Phase: ";
281 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
282 std::cout << fs.fugacity(phaseIdx, compIdx) << " ";
283 }
284 std::cout << "\n";
285 }
286
287 std::cout << "global component molarities: ";
288 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
289 Scalar sum = 0;
290 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
291 sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
292 }
293 std::cout << sum << " ";
294 }
295 std::cout << "\n";
296 }
297
298 template <class MaterialLaw, class FluidState>
299 1744 void linearize_(Matrix &J,
300 Vector &b,
301 FluidState &fluidState,
302 ParameterCache &paramCache,
303 const MaterialLaw& material,
304 const ComponentVector &globalMolarities)
305 {
306 1744 FluidState origFluidState(fluidState);
307 1662 ParameterCache origParamCache(paramCache);
308
309 1744 Vector tmp;
310
311 // reset jacobian
312 1744 J = 0;
313
314 1744 calculateDefect_(b, fluidState, fluidState, globalMolarities);
315
316 // assemble jacobian matrix
317
2/2
✓ Branch 1 taken 40380 times.
✓ Branch 2 taken 1744 times.
 42124 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
318 ////////
319 // approximately calculate partial derivatives of the
320 // fugacity defect of all components in regard to the mole
321 // fraction of the i-th component. This is done via
322 // forward differences
323
324 // deviate the mole fraction of the i-th component
325 40380 Scalar x_i = getQuantity_(fluidState, pvIdx);
326 80760 const Scalar eps = 1e-8/quantityWeight_(fluidState, pvIdx);
327 40380 setQuantity_(fluidState, paramCache, material, pvIdx, x_i + eps);
328
329 // compute derivative of the defect
330 40380 calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
331 40380 tmp -= b;
332 tmp /= eps;
333
334 // store derivative in jacobian matrix
335
2/2
✓ Branch 0 taken 960264 times.
✓ Branch 1 taken 40380 times.
 1000644 for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
336 3841056 J[eqIdx][pvIdx] = tmp[eqIdx];
337
338 // fluid state and parameter cache to their original values
339 40380 fluidState = origFluidState;
340 39888 paramCache = origParamCache;
341
342 // end forward differences
343 ////////
344 }
345 1744 }
346
347 template <class FluidState>
348 42124 void calculateDefect_(Vector &b,
349 const FluidState &fluidStateEval,
350 const FluidState &fluidState,
351 const ComponentVector &globalMolarities)
352 {
353 42124 int eqIdx = 0;
354
355 // fugacity of any component must be equal in all phases
356
2/2
✓ Branch 0 taken 291998 times.
✓ Branch 1 taken 42124 times.
 334122 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
357
2/2
✓ Branch 0 taken 582848 times.
✓ Branch 1 taken 291998 times.
 874846 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
358 1165696 b[eqIdx] =
359 1748544 fluidState.fugacity(/*phaseIdx=*/0, compIdx) -
360 fluidState.fugacity(phaseIdx, compIdx);
361 582848 ++eqIdx;
362 }
363 }
364
365
1/2
✓ Branch 0 taken 42124 times.
✗ Branch 1 not taken.
 42124 assert(eqIdx == numComponents*(numPhases - 1));
366
367 // the fact saturations must sum up to 1 is included explicitly!
368
369 // capillary pressures are explicitly included!
370
371 // global molarities are given
372
2/2
✓ Branch 0 taken 42124 times.
✓ Branch 1 taken 291998 times.
 334122 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
373 291998 b[eqIdx] = 0.0;
374
2/2
✓ Branch 0 taken 874846 times.
✓ Branch 1 taken 291998 times.
 1166844 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
375 874846 b[eqIdx] +=
376 fluidState.saturation(phaseIdx)
377 2624538 * fluidState.molarity(phaseIdx, compIdx);
378 }
379
380 583996 b[eqIdx] -= globalMolarities[compIdx];
381 291998 ++eqIdx;
382 }
383
384 // model assumptions (-> non-linear complementarity functions)
385 // must be adhered
386
2/2
✓ Branch 0 taken 125798 times.
✓ Branch 1 taken 42124 times.
 167922 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
387 Scalar sumMoleFracEval = 0.0;
388
2/2
✓ Branch 0 taken 874846 times.
✓ Branch 1 taken 125798 times.
 1000644 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
389 1749692 sumMoleFracEval += fluidStateEval.moleFraction(phaseIdx, compIdx);
390
391
4/4
✓ Branch 0 taken 82760 times.
✓ Branch 1 taken 43038 times.
✓ Branch 2 taken 82760 times.
✓ Branch 3 taken 43038 times.
 251596 if (1.0 - sumMoleFracEval > fluidStateEval.saturation(phaseIdx)) {
392 129114 b[eqIdx] = fluidState.saturation(phaseIdx);
393 }
394 else {
395 Scalar sumMoleFrac = 0.0;
396
2/2
✓ Branch 0 taken 574770 times.
✓ Branch 1 taken 82760 times.
 657530 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
397 1149540 sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
398 165520 b[eqIdx] = 1.0 - sumMoleFrac;
399 }
400
401 125798 ++eqIdx;
402 }
403 42124 }
404
405 template <class MaterialLaw, class FluidState>
406 1744 Scalar update_(FluidState &fluidState,
407 ParameterCache &paramCache,
408 const MaterialLaw& material,
409 const Vector &deltaX)
410 {
411 1744 Scalar relError = 0;
412
2/2
✓ Branch 0 taken 40380 times.
✓ Branch 1 taken 1744 times.
 42124 for (int pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
413 40380 Scalar tmp = getQuantity_(fluidState, pvIdx);
414
2/2
✓ Branch 0 taken 38636 times.
✓ Branch 1 taken 1744 times.
 40380 Scalar delta = deltaX[pvIdx];
415 using std::max;
416 using std::abs;
417 using std::min;
418
4/4
✓ Branch 0 taken 38636 times.
✓ Branch 1 taken 1744 times.
✓ Branch 2 taken 2418 times.
✓ Branch 3 taken 37962 times.
 79016 relError = max(relError, abs(delta)*quantityWeight_(fluidState, pvIdx));
419
420
2/2
✓ Branch 0 taken 82 times.
✓ Branch 1 taken 410 times.
 80268 if (isSaturationIdx_(pvIdx)) {
421 // dampen to at most 20% change in saturation per
422 // iteration
423
4/4
✓ Branch 0 taken 3394 times.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 3398 times.
✓ Branch 3 taken 8 times.
 10198 delta = min(0.2, max(-0.2, delta));
424 }
425 73948 else if (isMoleFracIdx_(pvIdx)) {
426 // dampen to at most 15% change in mole fraction per
427 // iteration
428
4/4
✓ Branch 0 taken 35146 times.
✓ Branch 1 taken 84 times.
✓ Branch 2 taken 35072 times.
✓ Branch 3 taken 158 times.
 105448 delta = min(0.15, max(-0.15, delta));
429 }
430
1/2
✓ Branch 0 taken 1744 times.
✗ Branch 1 not taken.
 1744 else if (isPressureIdx_(pvIdx)) {
431 // dampen to at most 15% change in pressure per
432 // iteration
433
8/8
✓ Branch 0 taken 1534 times.
✓ Branch 1 taken 210 times.
✓ Branch 2 taken 1534 times.
✓ Branch 3 taken 210 times.
✓ Branch 4 taken 1641 times.
✓ Branch 5 taken 103 times.
✓ Branch 6 taken 1641 times.
✓ Branch 7 taken 103 times.
 8197 delta = min(0.15*fluidState.pressure(0), max(-0.15*fluidState.pressure(0), delta));
434 }
435
436 40380 setQuantityRaw_(fluidState, pvIdx, tmp - delta);
437 }
438
439 // make sure all saturations, pressures and mole fractions are non-negative
440 Scalar sumSat = 0;
441
2/2
✓ Branch 0 taken 5150 times.
✓ Branch 1 taken 1744 times.
 6894 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
442
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 5148 times.
 5150 Scalar value = fluidState.saturation(phaseIdx);
443
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 5148 times.
 5150 if (value < -0.05) {
444 2 value = -0.05;
445 2 fluidState.setSaturation(phaseIdx, value);
446 }
447 5150 sumSat += value;
448
449
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5150 times.
 5150 value = fluidState.pressure(phaseIdx);
450
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5150 times.
 5150 if (value < 0)
451 ✗ fluidState.setPressure(phaseIdx, 0.0);
452
453 bool allMoleFractionsAreZero = true;
454
2/2
✓ Branch 0 taken 35230 times.
✓ Branch 1 taken 5150 times.
 40380 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
455
2/2
✓ Branch 0 taken 3334 times.
✓ Branch 1 taken 31896 times.
 35230 value = fluidState.moleFraction(phaseIdx, compIdx);
456
2/2
✓ Branch 0 taken 3334 times.
✓ Branch 1 taken 31896 times.
 35230 if (value > 0)
457 allMoleFractionsAreZero = false;
458
2/2
✓ Branch 0 taken 2433 times.
✓ Branch 1 taken 901 times.
 3334 else if (value < 0)
459 2433 fluidState.setMoleFraction(phaseIdx, compIdx, 0.0);
460 }
461
462 // set at least one mole fraction to a positive value
463 // to prevent breakdown of the linear solver in completeFluidState_
464
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 5147 times.
 5150 if (allMoleFractionsAreZero)
465 3 fluidState.setMoleFraction(phaseIdx, /*compIdx=*/0, 1e-10);
466 }
467
468 // last saturation
469
2/2
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 1727 times.
 1744 if (sumSat > 1.05) {
470
2/2
✓ Branch 0 taken 45 times.
✓ Branch 1 taken 17 times.
 62 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
471 45 Scalar value = fluidState.saturation(phaseIdx)/(0.95*sumSat);
472 90 fluidState.setSaturation(phaseIdx, value);
473 }
474 }
475
476 1744 completeFluidState_(fluidState, paramCache, material);
477
478 1744 return relError;
479 }
480
481 template <class MaterialLaw, class FluidState>
482 2408 void completeFluidState_(FluidState &fluidState,
483 ParameterCache &paramCache,
484 const MaterialLaw& material)
485 {
486 // calculate the saturation of the last phase as a function of
487 // the other saturations
488 2408 Scalar sumSat = 0.0;
489
2/2
✓ Branch 0 taken 4730 times.
✓ Branch 1 taken 2408 times.
 7138 for (int phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
490 9460 sumSat += fluidState.saturation(phaseIdx);
491 2408 fluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
492
493 // update the pressures using the material law (saturations
494 // and first pressure are already set because it is implicitly
495 // solved for.)
496 2408 const auto pc = material.capillaryPressures(fluidState, wettingPhaseIdx_);
497
2/2
✓ Branch 0 taken 4730 times.
✓ Branch 1 taken 2408 times.
 7138 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
498 4730 fluidState.setPressure(phaseIdx,
499 fluidState.pressure(0)
500 18920 + (pc[phaseIdx] - pc[0]));
501
502 // update the parameter cache
503 paramCache.updateAll(fluidState, /*except=*/ParameterCache::Temperature);
504
505 // update all densities and fugacity coefficients
506
2/2
✓ Branch 0 taken 7138 times.
✓ Branch 1 taken 2408 times.
 9546 for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
507
508 7138 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
509
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6966 times.
 7138 Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
510 7138 fluidState.setDensity(phaseIdx, rho);
511 7138 fluidState.setMolarDensity(phaseIdx, rhoMolar);
512
513
2/2
✓ Branch 0 taken 49106 times.
✓ Branch 1 taken 7138 times.
 56244 for (int compIdx = 0; compIdx < numComponents; ++ compIdx) {
514 49106 Scalar phi = FluidSystem::fugacityCoefficient( fluidState, paramCache, phaseIdx, compIdx);
515 98212 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
516 }
517 }
518 2408 }
519
520 ✗ bool isPressureIdx_(int pvIdx)
521 ✗ { return pvIdx == 0; }
522
523 ✗ bool isSaturationIdx_(int pvIdx)
524
2/2
✓ Branch 0 taken 3324 times.
✓ Branch 1 taken 36564 times.
 39888 { return 1 <= pvIdx && pvIdx < numPhases; }
525
526 ✗ bool isMoleFracIdx_(int pvIdx)
527
2/2
✓ Branch 0 taken 35230 times.
✓ Branch 1 taken 1744 times.
 36974 { return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
528
529 // retrieves a quantity from the fluid state
530 template <class FluidState>
531 ✗ Scalar getQuantity_(const FluidState &fs, int pvIdx)
532 {
533 ✗ assert(pvIdx < numEq);
534
535 // first pressure
536 ✗ if (pvIdx < 1) {
537 ✗ int phaseIdx = 0;
538 ✗ return fs.pressure(phaseIdx);
539 }
540 // first M - 1 saturations
541 ✗ else if (pvIdx < numPhases) {
542 ✗ int phaseIdx = pvIdx - 1;
543 ✗ return fs.saturation(phaseIdx);
544 }
545 // mole fractions
546 else // if (pvIdx < numPhases + numPhases*numComponents)
547 {
548 ✗ int phaseIdx = (pvIdx - numPhases)/numComponents;
549 ✗ int compIdx = (pvIdx - numPhases)%numComponents;
550 ✗ return fs.moleFraction(phaseIdx, compIdx);
551 }
552 }
553
554 // set a quantity in the fluid state
555 template <class MaterialLaw, class FluidState>
556 40380 void setQuantity_(FluidState &fs,
557 ParameterCache &paramCache,
558 const MaterialLaw& material,
559 int pvIdx,
560 Scalar value)
561 {
562
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40380 times.
 40380 assert(0 <= pvIdx && pvIdx < numEq);
563
564
2/2
✓ Branch 0 taken 1744 times.
✓ Branch 1 taken 38636 times.
 40380 if (pvIdx < 1) { // <- first pressure
565 1744 Scalar delta = value - fs.pressure(0);
566
567 // set all pressures. here we assume that the capillary
568 // pressure does not depend on absolute pressure.
569
2/2
✓ Branch 0 taken 5068 times.
✓ Branch 1 taken 1826 times.
 6894 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
570 15450 fs.setPressure(phaseIdx, fs.pressure(phaseIdx) + delta);
571 paramCache.updateAllPressures(fs);
572
573 // update all densities and fugacity coefficients
574
2/2
✓ Branch 0 taken 5150 times.
✓ Branch 1 taken 1744 times.
 6894 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
575 5150 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
576
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4986 times.
 5150 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
577 5150 fs.setDensity(phaseIdx, rho);
578 5150 fs.setMolarDensity(phaseIdx, rhoMolar);
579
580
2/2
✓ Branch 0 taken 35230 times.
✓ Branch 1 taken 5150 times.
 40380 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
581 35230 Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
582 70460 fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
583 }
584 }
585 }
586
2/2
✓ Branch 0 taken 3406 times.
✓ Branch 1 taken 35230 times.
 38636 else if (pvIdx < numPhases) { // <- first M - 1 saturations
587 3406 Scalar delta = value - fs.saturation(/*phaseIdx=*/pvIdx - 1);
588 3406 fs.setSaturation(/*phaseIdx=*/pvIdx - 1, value);
589
590 // set last saturation (-> minus the change of the
591 // satuation of the other phase)
592 6812 fs.setSaturation(/*phaseIdx=*/numPhases - 1,
593 6812 fs.saturation(numPhases - 1) - delta);
594
595 // update all fluid pressures using the capillary pressure
596 // law
597 3406 const auto pc = material.capillaryPressures(fs, wettingPhaseIdx_);
598
2/2
✓ Branch 0 taken 6730 times.
✓ Branch 1 taken 3406 times.
 10136 for (int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
599 6730 fs.setPressure(phaseIdx,
600 fs.pressure(0)
601 26920 + (pc[phaseIdx] - pc[0]));
602 paramCache.updateAllPressures(fs);
603
604 // update all densities and fugacity coefficients
605
2/2
✓ Branch 0 taken 10136 times.
✓ Branch 1 taken 3406 times.
 13542 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
606 10136 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
607
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9972 times.
 10136 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
608 10136 fs.setDensity(phaseIdx, rho);
609 10136 fs.setMolarDensity(phaseIdx, rhoMolar);
610
611
2/2
✓ Branch 0 taken 70132 times.
✓ Branch 1 taken 10136 times.
 80268 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
612 70132 Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx);
613 140264 fs.setFugacityCoefficient(phaseIdx, compIdx, phi);
614 }
615 }
616 }
617 else if (pvIdx < numPhases + numPhases*numComponents) // <- mole fractions
618 {
619 35230 int phaseIdx = (pvIdx - numPhases)/numComponents;
620 35230 int compIdx = (pvIdx - numPhases)%numComponents;
621
622 35230 fs.setMoleFraction(phaseIdx, compIdx, value);
623 35230 paramCache.updateSingleMoleFraction(fs, phaseIdx, compIdx);
624
625 // update the density of the phase
626 35230 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
627
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34902 times.
 35230 Scalar rhoMolar = FluidSystem::molarDensity(fs, paramCache, phaseIdx);
628
2/2
✓ Branch 0 taken 23268 times.
✓ Branch 1 taken 11634 times.
 35230 fs.setDensity(phaseIdx, rho);
629
2/2
✓ Branch 0 taken 23268 times.
✓ Branch 1 taken 11634 times.
 35230 fs.setMolarDensity(phaseIdx, rhoMolar);
630
631 // if the phase's fugacity coefficients are composition
632 // dependent, update them as well.
633
2/2
✓ Branch 0 taken 23268 times.
✓ Branch 1 taken 11634 times.
 35230 if (!FluidSystem::isIdealMixture(phaseIdx)) {
634
2/2
✓ Branch 0 taken 162876 times.
✓ Branch 1 taken 23268 times.
 186144 for (int compIdx2 = 0; compIdx2 < numComponents; ++compIdx2) {
635 162876 Scalar phi = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx2);
636 325752 fs.setFugacityCoefficient(phaseIdx, compIdx2, phi);
637 }
638 }
639 }
640 else {
641 assert(false);
642 }
643 40380 }
644
645 // set a quantity in the fluid state
646 template <class FluidState>
647 40380 void setQuantityRaw_(FluidState &fs, int pvIdx, Scalar value)
648 {
649
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40380 times.
 40380 assert(pvIdx < numEq);
650
651 // first pressure
652
2/2
✓ Branch 0 taken 1744 times.
✓ Branch 1 taken 38636 times.
 40380 if (pvIdx < 1) {
653 1744 int phaseIdx = 0;
654 1744 fs.setPressure(phaseIdx, value);
655 }
656 // first M - 1 saturations
657
2/2
✓ Branch 0 taken 3406 times.
✓ Branch 1 taken 35230 times.
 38636 else if (pvIdx < numPhases) {
658 3406 int phaseIdx = pvIdx - 1;
659 3406 fs.setSaturation(phaseIdx, value);
660 }
661 // mole fractions
662 else // if (pvIdx < numPhases + numPhases*numComponents)
663 {
664 35230 int phaseIdx = (pvIdx - numPhases)/numComponents;
665 35230 int compIdx = (pvIdx - numPhases)%numComponents;
666 35230 fs.setMoleFraction(phaseIdx, compIdx, value);
667 }
668 40380 }
669
670 template <class FluidState>
671 ✗ Scalar quantityWeight_(const FluidState &fs, int pvIdx)
672 {
673 // first pressure
674
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 38636 times.
✓ Branch 3 taken 1744 times.
✓ Branch 4 taken 38636 times.
✓ Branch 5 taken 1744 times.
 80760 if (pvIdx < 1)
675 ✗ return 1.0/1e5;
676 // first M - 1 saturations
677 ✗ else if (pvIdx < numPhases)
678 ✗ return 1.0;
679 // mole fractions
680 else // if (pvIdx < numPhases + numPhases*numComponents)
681 ✗ return 1.0;
682 }
683
684 private:
685 int wettingPhaseIdx_ = 0; //!< the phase index of the wetting phase
686 };
687
688 } // end namespace Dumux
689
690 #endif
691