| 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 fluid composition given the component | ||
| 11 | * fugacities and an arbitrary equation of state. | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_COMPOSITION_FROM_FUGACITIES_HH | ||
| 14 | #define DUMUX_COMPOSITION_FROM_FUGACITIES_HH | ||
| 15 | |||
| 16 | #include <dune/common/fvector.hh> | ||
| 17 | #include <dune/common/fmatrix.hh> | ||
| 18 | |||
| 19 | #include <dumux/common/exceptions.hh> | ||
| 20 | |||
| 21 | namespace Dumux { | ||
| 22 | |||
| 23 | /*! | ||
| 24 | * \ingroup ConstraintSolvers | ||
| 25 | * \brief Calculates the chemical equilibrium from the component | ||
| 26 | * fugacities \f$ f^\kappa \f$ in the phase \f$ \alpha \f$. | ||
| 27 | * | ||
| 28 | * This constraint solver takes the component fugacity \f$f^\kappa\f$ of | ||
| 29 | * of component \f$\kappa\f$, the temperature \f$ T_\alpha \f$, the pressure | ||
| 30 | * \f$p_\alpha\f$ and the composition \f$x^\lambda_\alpha\f$ of a phase | ||
| 31 | * \f$\alpha\f$ as input and calculates the mole fraction of component | ||
| 32 | * \f$\kappa\f$ in that fluid phase \f$x^\kappa_\alpha\f$. This means | ||
| 33 | * that the thermodynamic constraints used by this solver are | ||
| 34 | * | ||
| 35 | * \f$ f^\kappa = \Phi^\kappa_\alpha(\{x^\lambda_\alpha \}, T_\alpha, p_\alpha) p_\alpha x^\kappa_\alpha\; \f$, | ||
| 36 | * | ||
| 37 | * where \f${f^\kappa}\f$, \f$ T_\alpha \f$ and \f$ p_\alpha \f$ are fixed values. | ||
| 38 | */ | ||
| 39 | template <class Scalar, class FluidSystem> | ||
| 40 | class CompositionFromFugacities | ||
| 41 | { | ||
| 42 | static constexpr int numComponents = FluidSystem::numComponents; | ||
| 43 | using ParameterCache = typename FluidSystem::ParameterCache; | ||
| 44 | |||
| 45 | public: | ||
| 46 | using ComponentVector = Dune::FieldVector<Scalar, numComponents>; | ||
| 47 | |||
| 48 | /*! | ||
| 49 | * \brief Guess an initial value for the composition of the phase. | ||
| 50 | * \param fluidState Thermodynamic state of the fluids | ||
| 51 | * \param paramCache Container for cache parameters | ||
| 52 | * \param phaseIdx The phase index | ||
| 53 | * \param fugVec fugacity vector of the component | ||
| 54 | */ | ||
| 55 | template <class FluidState> | ||
| 56 | ✗ | static void guessInitial(FluidState &fluidState, | |
| 57 | ParameterCache ¶mCache, | ||
| 58 | int phaseIdx, | ||
| 59 | const ComponentVector &fugVec) | ||
| 60 | { | ||
| 61 |
2/5✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 15648930 times.
✓ Branch 4 taken 1943 times.
|
18038473 | if (!FluidSystem::isIdealMixture(phaseIdx)) |
| 62 | { | ||
| 63 | // Pure component fugacities | ||
| 64 | for (unsigned int i = 0; i < numComponents; ++ i) | ||
| 65 | { | ||
| 66 | fluidState.setMoleFraction(phaseIdx, | ||
| 67 | i, | ||
| 68 | 1.0/numComponents); | ||
| 69 | } | ||
| 70 | } | ||
| 71 | ✗ | } | |
| 72 | |||
| 73 | /*! | ||
| 74 | * \brief Calculates the chemical equilibrium from the component | ||
| 75 | * fugacities in a phase. | ||
| 76 | * \param fluidState Thermodynamic state of the fluids | ||
| 77 | * \param paramCache Container for cache parameters | ||
| 78 | * \param phaseIdx The phase index | ||
| 79 | * \param targetFug target fugacity | ||
| 80 | * | ||
| 81 | * The phase's fugacities must already be set. | ||
| 82 | */ | ||
| 83 | template <class FluidState> | ||
| 84 | 18038473 | static void solve(FluidState &fluidState, | |
| 85 | ParameterCache ¶mCache, | ||
| 86 | int phaseIdx, | ||
| 87 | const ComponentVector &targetFug) | ||
| 88 | { | ||
| 89 | // use a much more efficient method in case the phase is an | ||
| 90 | // ideal mixture | ||
| 91 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18038473 times.
|
18038473 | if (FluidSystem::isIdealMixture(phaseIdx)) { |
| 92 | 18038473 | solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug); | |
| 93 | return; | ||
| 94 | } | ||
| 95 | |||
| 96 | // save initial composition in case something goes wrong | ||
| 97 | Dune::FieldVector<Scalar, numComponents> xInit; | ||
| 98 | for (int i = 0; i < numComponents; ++i) { | ||
| 99 | xInit[i] = fluidState.moleFraction(phaseIdx, i); | ||
| 100 | } | ||
| 101 | |||
| 102 | ///////////////////////// | ||
| 103 | // Newton method | ||
| 104 | ///////////////////////// | ||
| 105 | |||
| 106 | // Jacobian matrix | ||
| 107 | Dune::FieldMatrix<Scalar, numComponents, numComponents> J; | ||
| 108 | // solution, i.e. phase composition | ||
| 109 | Dune::FieldVector<Scalar, numComponents> x; | ||
| 110 | // right hand side | ||
| 111 | Dune::FieldVector<Scalar, numComponents> b; | ||
| 112 | |||
| 113 | paramCache.updatePhase(fluidState, phaseIdx); | ||
| 114 | |||
| 115 | // maximum number of iterations | ||
| 116 | const int nMax = 25; | ||
| 117 | for (int nIdx = 0; nIdx < nMax; ++nIdx) { | ||
| 118 | // calculate Jacobian matrix and right hand side | ||
| 119 | linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug); | ||
| 120 | |||
| 121 | // Solve J*x = b | ||
| 122 | x = 0; | ||
| 123 | try { J.solve(x, b); } | ||
| 124 | catch (Dune::FMatrixError &e) | ||
| 125 | { throw NumericalProblem(e.what()); } | ||
| 126 | |||
| 127 | // update the fluid composition. b is also used to store | ||
| 128 | // the defect for the next iteration. | ||
| 129 | Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug); | ||
| 130 | |||
| 131 | if (relError < 1e-9) { | ||
| 132 | Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx); | ||
| 133 | Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx); | ||
| 134 | fluidState.setDensity(phaseIdx, rho); | ||
| 135 | fluidState.setMolarDensity(phaseIdx, rhoMolar); | ||
| 136 | |||
| 137 | //std::cout << "num iterations: " << nIdx << "\n"; | ||
| 138 | return; | ||
| 139 | } | ||
| 140 | } | ||
| 141 | |||
| 142 | DUNE_THROW(NumericalProblem, | ||
| 143 | "Calculating the " << FluidSystem::phaseName(phaseIdx) | ||
| 144 | << "Phase composition failed. Initial {x} = {" | ||
| 145 | << xInit | ||
| 146 | << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx) | ||
| 147 | << ", T = " << fluidState.temperature(phaseIdx)); | ||
| 148 | } | ||
| 149 | |||
| 150 | |||
| 151 | protected: | ||
| 152 | // update the phase composition in case the phase is an ideal | ||
| 153 | // mixture, i.e. the component's fugacity coefficients are | ||
| 154 | // independent of the phase's composition. | ||
| 155 | template <class FluidState> | ||
| 156 | 18038473 | static void solveIdealMix_(FluidState &fluidState, | |
| 157 | ParameterCache ¶mCache, | ||
| 158 | int phaseIdx, | ||
| 159 | const ComponentVector &fugacities) | ||
| 160 | { | ||
| 161 |
2/2✓ Branch 0 taken 54699773 times.
✓ Branch 1 taken 18038473 times.
|
72738246 | for (int i = 0; i < numComponents; ++ i) { |
| 162 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
54699773 | Scalar phi = FluidSystem::fugacityCoefficient(fluidState, |
| 163 | paramCache, | ||
| 164 | phaseIdx, | ||
| 165 | i); | ||
| 166 | 54699773 | Scalar gamma = phi * fluidState.pressure(phaseIdx); | |
| 167 | 54699773 | fluidState.setFugacityCoefficient(phaseIdx, i, phi); | |
| 168 | 139329720 | fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma); | |
| 169 | } | ||
| 170 | |||
| 171 | 18038473 | paramCache.updatePhase(fluidState, phaseIdx); | |
| 172 | |||
| 173 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
18038473 | Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx); |
| 174 | 18038473 | fluidState.setDensity(phaseIdx, rho); | |
| 175 | 18038473 | Scalar rhoMolar = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx); | |
| 176 | 36076946 | fluidState.setMolarDensity(phaseIdx, rhoMolar); | |
| 177 | 18038473 | } | |
| 178 | |||
| 179 | template <class FluidState> | ||
| 180 | ✗ | static void linearize_(Dune::FieldMatrix<Scalar, numComponents, numComponents> &J, | |
| 181 | Dune::FieldVector<Scalar, numComponents> &defect, | ||
| 182 | FluidState &fluidState, | ||
| 183 | ParameterCache ¶mCache, | ||
| 184 | int phaseIdx, | ||
| 185 | const ComponentVector &targetFug) | ||
| 186 | { | ||
| 187 | // reset jacobian | ||
| 188 | ✗ | J = 0; | |
| 189 | |||
| 190 | // calculate the defect (deviation of the current fugacities | ||
| 191 | // from the target fugacities) | ||
| 192 | ✗ | for (int i = 0; i < numComponents; ++ i) { | |
| 193 | ✗ | Scalar phi = FluidSystem::fugacityCoefficient(fluidState, | |
| 194 | paramCache, | ||
| 195 | phaseIdx, | ||
| 196 | i); | ||
| 197 | ✗ | Scalar f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i); | |
| 198 | ✗ | fluidState.setFugacityCoefficient(phaseIdx, i, phi); | |
| 199 | |||
| 200 | ✗ | defect[i] = targetFug[i] - f; | |
| 201 | } | ||
| 202 | |||
| 203 | // assemble jacobian matrix of the constraints for the composition | ||
| 204 | ✗ | for (int i = 0; i < numComponents; ++ i) { | |
| 205 | ✗ | const Scalar eps = 1e-11; //std::max(1e-16, std::abs(x_i)*1e-9); | |
| 206 | |||
| 207 | //////// | ||
| 208 | // approximately calculate partial derivatives of the | ||
| 209 | // fugacity defect of all components in regard to the mole | ||
| 210 | // fraction of the i-th component. This is done via | ||
| 211 | // forward differences | ||
| 212 | |||
| 213 | // deviate the mole fraction of the i-th component | ||
| 214 | ✗ | Scalar x_i = fluidState.moleFraction(phaseIdx, i); | |
| 215 | ✗ | fluidState.setMoleFraction(phaseIdx, i, x_i + eps); | |
| 216 | ✗ | paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i); | |
| 217 | |||
| 218 | // compute new defect and derivative for all component | ||
| 219 | // fugacities | ||
| 220 | ✗ | for (int j = 0; j < numComponents; ++j) { | |
| 221 | // compute the j-th component's fugacity coefficient ... | ||
| 222 | ✗ | Scalar phi = FluidSystem::fugacityCoefficient(fluidState, | |
| 223 | paramCache, | ||
| 224 | phaseIdx, | ||
| 225 | j); | ||
| 226 | // ... and its fugacity ... | ||
| 227 | ✗ | Scalar f = | |
| 228 | ✗ | phi * | |
| 229 | ✗ | fluidState.pressure(phaseIdx) * | |
| 230 | fluidState.moleFraction(phaseIdx, j); | ||
| 231 | // as well as the defect for this fugacity | ||
| 232 | ✗ | Scalar defJPlusEps = targetFug[j] - f; | |
| 233 | |||
| 234 | // use forward differences to calculate the defect's | ||
| 235 | // derivative | ||
| 236 | ✗ | J[j][i] = (defJPlusEps - defect[j])/eps; | |
| 237 | } | ||
| 238 | |||
| 239 | // reset composition to original value | ||
| 240 | ✗ | fluidState.setMoleFraction(phaseIdx, i, x_i); | |
| 241 | ✗ | paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i); | |
| 242 | |||
| 243 | // end forward differences | ||
| 244 | //////// | ||
| 245 | } | ||
| 246 | ✗ | } | |
| 247 | |||
| 248 | template <class FluidState> | ||
| 249 | ✗ | static Scalar update_(FluidState &fluidState, | |
| 250 | ParameterCache ¶mCache, | ||
| 251 | Dune::FieldVector<Scalar, numComponents> &x, | ||
| 252 | Dune::FieldVector<Scalar, numComponents> &b, | ||
| 253 | int phaseIdx, | ||
| 254 | const Dune::FieldVector<Scalar, numComponents> &targetFug) | ||
| 255 | { | ||
| 256 | // store original composition and calculate relative error | ||
| 257 | ✗ | Dune::FieldVector<Scalar, numComponents> origComp; | |
| 258 | ✗ | Scalar relError = 0; | |
| 259 | ✗ | Scalar sumDelta = 0; | |
| 260 | using std::max; | ||
| 261 | using std::abs; | ||
| 262 | using std::min; | ||
| 263 | ✗ | for (unsigned int i = 0; i < numComponents; ++i) | |
| 264 | { | ||
| 265 | ✗ | origComp[i] = fluidState.moleFraction(phaseIdx, i); | |
| 266 | ✗ | relError = max(relError, abs(x[i])); | |
| 267 | ✗ | sumDelta += abs(x[i]); | |
| 268 | } | ||
| 269 | |||
| 270 | // chop update to at most 20% change in composition | ||
| 271 | ✗ | const Scalar maxDelta = 0.2; | |
| 272 | ✗ | if (sumDelta > maxDelta) | |
| 273 | ✗ | x /= (sumDelta/maxDelta); | |
| 274 | |||
| 275 | // change composition | ||
| 276 | ✗ | for (unsigned int i = 0; i < numComponents; ++i) | |
| 277 | { | ||
| 278 | ✗ | Scalar newx = origComp[i] - x[i]; | |
| 279 | |||
| 280 | // only allow negative mole fractions if the target fugacity is negative | ||
| 281 | ✗ | if (targetFug[i] > 0) | |
| 282 | ✗ | newx = max(0.0, newx); | |
| 283 | // only allow positive mole fractions if the target fugacity is positive | ||
| 284 | ✗ | else if (targetFug[i] < 0) | |
| 285 | ✗ | newx = min(0.0, newx); | |
| 286 | // if the target fugacity is zero, the mole fraction must also be zero | ||
| 287 | else | ||
| 288 | ✗ | newx = 0; | |
| 289 | |||
| 290 | ✗ | fluidState.setMoleFraction(phaseIdx, i, newx); | |
| 291 | //sumMoleFrac += abs(newx); | ||
| 292 | } | ||
| 293 | |||
| 294 | ✗ | paramCache.updateComposition(fluidState, phaseIdx); | |
| 295 | |||
| 296 | ✗ | return relError; | |
| 297 | } | ||
| 298 | |||
| 299 | template <class FluidState> | ||
| 300 | static Scalar calculateDefect_(const FluidState ¶ms, | ||
| 301 | int phaseIdx, | ||
| 302 | const ComponentVector &targetFug) | ||
| 303 | { | ||
| 304 | Scalar result = 0.0; | ||
| 305 | using std::abs; | ||
| 306 | for (int i = 0; i < numComponents; ++i) { | ||
| 307 | // sum of the fugacity defect weighted by the inverse | ||
| 308 | // fugacity coefficient | ||
| 309 | result += abs( | ||
| 310 | (targetFug[i] - params.fugacity(phaseIdx, i)) | ||
| 311 | / | ||
| 312 | params.fugacityCoeff(phaseIdx, i) ); | ||
| 313 | } | ||
| 314 | return result; | ||
| 315 | } | ||
| 316 | }; | ||
| 317 | |||
| 318 | } // end namespace Dumux | ||
| 319 | |||
| 320 | #endif | ||
| 321 |