GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/material/constraintsolvers/misciblemultiphasecomposition.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 44 48 91.7%
Functions: 11 12 91.7%
Branches: 27 46 58.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 Computes the composition of all phases of a N-phase,
11 * N-component fluid system assuming that all N phases are
12 * present
13 */
14 #ifndef DUMUX_MISCIBLE_MULTIPHASE_COMPOSITION_HH
15 #define DUMUX_MISCIBLE_MULTIPHASE_COMPOSITION_HH
16
17 #include <dune/common/fvector.hh>
18 #include <dune/common/fmatrix.hh>
19
20 #include <dumux/common/exceptions.hh>
21
22 namespace Dumux {
23 /*!
24 * \ingroup ConstraintSolvers
25 * \brief Computes the composition of all phases of a N-phase,
26 * N-component fluid system assuming that all N phases are
27 * present
28 *
29 * The constraint solver assumes the following quantities to be set:
30 *
31 * - temperatures of *all* phases
32 * - saturations of *all* phases
33 * - pressures of *all* phases
34 *
35 * It also assumes that the mole/mass fractions of all phases sum up
36 * to 1. After calling the solve() method the following quantities
37 * are calculated in addition:
38 *
39 * - temperature of *all* phases
40 * - density, molar density, molar volume of *all* phases
41 * - composition in mole and mass fractions and molarities of *all* phases
42 * - mean molar masses of *all* phases
43 * - fugacity coefficients of *all* components in *all* phases
44 */
45 template <class Scalar, class FluidSystem>
46 class MiscibleMultiPhaseComposition
47 {
48 static constexpr int numPhases = FluidSystem::numPhases;
49 static constexpr int numComponents = FluidSystem::numComponents;
50 static const int numMajorComponents = FluidSystem::numPhases;
51
52 public:
53 /*!
54 * \brief @copybrief Dumux::MiscibleMultiPhaseComposition
55 *
56 * This function additionally considers a lowering of the saturation vapor pressure
57 * of the wetting phase by the Kelvin equation:
58 * \f[
59 * p^\textrm{w}_\textrm{sat,Kelvin}
60 * = p^\textrm{w}_\textrm{sat}
61 * \exp \left( -\frac{p_\textrm{c}}{\varrho_\textrm{w} R_\textrm{w} T} \right)
62 * \f]
63 *
64 * \param fluidState A container with the current (physical) state of the fluid
65 * \param paramCache A container for iterative calculation of fluid composition
66 * \param knownPhaseIdx The index of the phase with known properties
67 */
68 template <class FluidState, class ParameterCache>
69 13908138 static void solve(FluidState &fluidState,
70 ParameterCache &paramCache,
71 int knownPhaseIdx = 0)
72 {
73 #ifndef NDEBUG
74 // currently this solver can only handle fluid systems which
75 // assume ideal mixtures of all fluids. TODO: relax this
76 // (requires solving a non-linear system of equations, i.e. using
77 // Newton method.)
78 13908138 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
79 assert(FluidSystem::isIdealMixture(phaseIdx));
80
81 }
82 #endif
83
84 //get the known mole fractions from the fluidState
85 //in a 2pnc system the n>2 mole fractions are primary variables and are already set in the fluidstate
86 ✗ Dune::FieldVector<Scalar, numComponents-numMajorComponents> xKnown(0.0);
87
2/2
✓ Branch 0 taken 11222029 times.
✓ Branch 1 taken 11222029 times.
 22444058 for (int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
88 {
89 22444058 xKnown[knownCompIdx] = fluidState.moleFraction(knownPhaseIdx, knownCompIdx + numMajorComponents);
90 }
91
92 // compute all fugacity coefficients
93
2/2
✓ Branch 0 taken 27816276 times.
✓ Branch 1 taken 13908138 times.
 41724414 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
94 paramCache.updatePhase(fluidState, phaseIdx);
95
96 // since we assume ideal mixtures, the fugacity
97 // coefficients of the components cannot depend on
98 // composition, i.e. the parameters in the cache are valid
99
2/2
✓ Branch 0 taken 78076610 times.
✓ Branch 1 taken 27816276 times.
 105892886 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
100
2/2
✓ Branch 0 taken 27207726 times.
✓ Branch 1 taken 27207726 times.
 78076610 Scalar fugCoeff = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
101 156153220 fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
102 }
103 }
104
105
106 // create the linear system of equations which defines the
107 // mole fractions
108 13908138 Dune::FieldMatrix<Scalar, numComponents*numPhases, numComponents*numPhases> M(0.0);
109 13908138 Dune::FieldVector<Scalar, numComponents*numPhases> x(0.0);
110 13908138 Dune::FieldVector<Scalar, numComponents*numPhases> b(0.0);
111
112 // assemble the equations expressing the assumption that the
113 // sum of all mole fractions in each phase must be 1
114
2/2
✓ Branch 0 taken 27816276 times.
✓ Branch 1 taken 13908138 times.
 41724414 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
115 27816276 int rowIdx = numComponents*(numPhases - 1) + phaseIdx;
116 27816276 b[rowIdx] = 1.0;
117
118
2/2
✓ Branch 0 taken 78076610 times.
✓ Branch 1 taken 27816276 times.
 105892886 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
119 78076610 int colIdx = phaseIdx*numComponents + compIdx;
120
121 234229830 M[rowIdx][colIdx] = 1.0;
122 }
123 }
124
125 // set the additional equations for the numComponents-numMajorComponents
126 // this is only relevant if numphases != numcomponents e.g. in a 2pnc system
127 // Components, of which the mole fractions are known, set to molefraction(knownCompIdx)=xKnown
128
2/2
✓ Branch 0 taken 11222029 times.
✓ Branch 1 taken 11222029 times.
 22444058 for(int knownCompIdx = 0; knownCompIdx < numComponents-numMajorComponents; ++knownCompIdx)
129 {
130 11222029 int rowIdx = numComponents + numPhases + knownCompIdx;
131 11222029 int colIdx = knownPhaseIdx*numComponents + knownCompIdx + numMajorComponents;
132 22444058 M[rowIdx][colIdx] = 1.0;
133 22444058 b[rowIdx] = xKnown[knownCompIdx];
134 }
135
136 // assemble the equations expressing the fact that the
137 // fugacities of each component are equal in all phases
138
2/2
✓ Branch 0 taken 39038305 times.
✓ Branch 1 taken 13908138 times.
 52946443 for (int compIdx = 0; compIdx < numComponents; ++compIdx)
139 {
140 39038305 int col1Idx = compIdx;
141 78076610 const auto entryPhase0 = fluidState.fugacityCoefficient(0, compIdx)*fluidState.pressure(0);
142
143
2/2
✓ Branch 0 taken 39038305 times.
✓ Branch 1 taken 39038305 times.
 78076610 for (unsigned int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
144 {
145 39038305 int rowIdx = (phaseIdx - 1)*numComponents + compIdx;
146 39038305 int col2Idx = phaseIdx*numComponents + compIdx;
147 78076610 M[rowIdx][col1Idx] = entryPhase0;
148 195191525 M[rowIdx][col2Idx] = -fluidState.fugacityCoefficient(phaseIdx, compIdx)*fluidState.pressure(phaseIdx);
149 }
150 }
151
152 // preconditioning of M to reduce condition number
153
2/2
✓ Branch 0 taken 39038305 times.
✓ Branch 1 taken 13908138 times.
 52946443 for (int compIdx = 0; compIdx < numComponents; compIdx++)
154 {
155 // Multiply row of main component (Raoult's Law) with 10e-5 (order of magn. of pressure)
156
2/2
✓ Branch 0 taken 22444058 times.
✓ Branch 1 taken 11222029 times.
 33666087 if (compIdx < numMajorComponents)
157 55632552 M[compIdx] *= 10e-5;
158
159 // Multiply row of sec. components (Henry's Law) with 10e-9 (order of magn. of Henry constant)
160 else
161 ✗ M[compIdx] *= 10e-9;
162 }
163
164 // solve for all mole fractions
165
1/2
✓ Branch 1 taken 13908138 times.
✗ Branch 2 not taken.
 13908138 try { M.solve(x, b); }
166 ✗ catch (const Dune::FMatrixError& e) {
167 ✗ DUNE_THROW(NumericalProblem,
168 "MiscibleMultiPhaseComposition: Failed to solve the linear equation system for the phase composition.");
169 }
170
171 // set all mole fractions and the the additional quantities in
172 // the fluid state
173
2/2
✓ Branch 0 taken 27816276 times.
✓ Branch 1 taken 13908138 times.
 41724414 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
174
2/2
✓ Branch 0 taken 78076610 times.
✓ Branch 1 taken 27816276 times.
 105892886 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
175 78076610 int rowIdx = phaseIdx*numComponents + compIdx;
176 215137642 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
177 }
178 27816276 paramCache.updateComposition(fluidState, phaseIdx);
179
180 27816276 Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
181 27816276 fluidState.setDensity(phaseIdx, value);
182
183 27816276 value = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
184 55632552 fluidState.setMolarDensity(phaseIdx, value);
185 }
186 13908138 }
187 };
188
189 } // end namespace Dumux
190
191 #endif
192