My Project
Loading...
Searching...
No Matches
CompositionFromFugacities.hpp
Go to the documentation of this file.
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 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
27#ifndef OPM_COMPOSITION_FROM_FUGACITIES_HPP
28#define OPM_COMPOSITION_FROM_FUGACITIES_HPP
29
31
34
35#include <dune/common/fvector.hh>
36#include <dune/common/fmatrix.hh>
37
38#include <cassert>
39#include <limits>
40#include <string_view>
41
42namespace Opm {
43
48template <class Scalar, class FluidSystem, class Evaluation = Scalar>
50{
51 enum { numComponents = FluidSystem::numComponents };
52
53public:
54 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
55
59 template <class FluidState>
60 static void guessInitial(FluidState& fluidState,
61 unsigned phaseIdx,
62 const ComponentVector& /*fugVec*/)
63 {
64 if (FluidSystem::isIdealMixture(phaseIdx))
65 return;
66
67 // Pure component fugacities
68 for (unsigned i = 0; i < numComponents; ++ i) {
69 //std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
70 fluidState.setMoleFraction(phaseIdx,
71 i,
72 1.0/numComponents);
73 }
74 }
75
82 template <class FluidState>
83 static void solve(FluidState& fluidState,
84 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
85 unsigned phaseIdx,
86 const ComponentVector& targetFug)
87 {
88 assert (phaseIdx < static_cast<unsigned int>(FluidSystem::numPhases));
89
90 // use a much more efficient method in case the phase is an
91 // ideal mixture
92 if (FluidSystem::isIdealMixture(phaseIdx)) {
93 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
94 return;
95 }
96
97 // save initial composition in case something goes wrong
98 Dune::FieldVector<Evaluation, numComponents> xInit;
99 for (unsigned i = 0; i < numComponents; ++i) {
100 xInit[i] = fluidState.moleFraction(phaseIdx, i);
101 }
102
104 // Newton method
106
107 // Jacobian matrix
108 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
109 // solution, i.e. phase composition
110 Dune::FieldVector<Evaluation, numComponents> x;
111 // right hand side
112 Dune::FieldVector<Evaluation, numComponents> b;
113
114 paramCache.updatePhase(fluidState, phaseIdx);
115
116 // maximum number of iterations
117 const int nMax = 25;
118 for (int nIdx = 0; nIdx < nMax; ++nIdx) {
119 // calculate Jacobian matrix and right hand side
120 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
121 Valgrind::CheckDefined(J);
122 Valgrind::CheckDefined(b);
123
124 /*
125 std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
126 for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
127 std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
128 std::cout << "\n";
129 std::cout << FluidSystem::phaseName(phaseIdx) << "Phase phi: ";
130 for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
131 std::cout << fluidState.fugacityCoefficient(phaseIdx, i) << " ";
132 std::cout << "\n";
133 */
134
135 // Solve J*x = b
136 x = 0.0;
137 try { J.solve(x, b); }
138 catch (const Dune::FMatrixError& e)
139 { throw NumericalProblem(e.what()); }
140
141 //std::cout << "original delta: " << x << "\n";
142
143 Valgrind::CheckDefined(x);
144
145 /*
146 std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
147 for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
148 std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
149 std::cout << "\n";
150 std::cout << "J: " << J << "\n";
151 std::cout << "rho: " << fluidState.density(phaseIdx) << "\n";
152 std::cout << "delta: " << x << "\n";
153 std::cout << "defect: " << b << "\n";
154
155 std::cout << "J: " << J << "\n";
156
157 std::cout << "---------------------------\n";
158 */
159
160 // update the fluid composition. b is also used to store
161 // the defect for the next iteration.
162 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
163
164 if (relError < 1e-9) {
165 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
166 fluidState.setDensity(phaseIdx, rho);
167
168 //std::cout << "num iterations: " << nIdx << "\n";
169 return;
170 }
171 }
172
173 auto cast = [](const auto d)
174 {
175#if HAVE_QUAD
176 if constexpr (std::is_same_v<decltype(d), const quad>)
177 return static_cast<double>(d);
178 else
179#endif
180 return d;
181 };
182
183 std::string msg =
184 std::string("Calculating the ")
185 + FluidSystem::phaseName(phaseIdx).data()
186 + "Phase composition failed. Initial {x} = {";
187 for (const auto& v : xInit)
188 msg += " " + std::to_string(cast(getValue(v)));
189 msg += " }, {fug_t} = {";
190 for (const auto& v : targetFug)
191 msg += " " + std::to_string(cast(getValue(v)));
192 msg += " }, p = " + std::to_string(cast(getValue(fluidState.pressure(phaseIdx))))
193 + ", T = " + std::to_string(cast(getValue(fluidState.temperature(phaseIdx))));
194 throw NumericalProblem(msg);
195 }
196
197
198protected:
199 // update the phase composition in case the phase is an ideal
200 // mixture, i.e. the component's fugacity coefficients are
201 // independent of the phase's composition.
202 template <class FluidState>
203 static void solveIdealMix_(FluidState& fluidState,
204 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
205 unsigned phaseIdx,
206 const ComponentVector& fugacities)
207 {
208 for (unsigned i = 0; i < numComponents; ++ i) {
209 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
210 paramCache,
211 phaseIdx,
212 i);
213 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
214 Valgrind::CheckDefined(phi);
215 Valgrind::CheckDefined(gamma);
216 Valgrind::CheckDefined(fugacities[i]);
217 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
218 fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
219 };
220
221 paramCache.updatePhase(fluidState, phaseIdx);
222
223 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
224 fluidState.setDensity(phaseIdx, rho);
225 return;
226 }
227
228 template <class FluidState>
229 static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents>& J,
230 Dune::FieldVector<Evaluation, numComponents>& defect,
231 FluidState& fluidState,
232 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
233 unsigned phaseIdx,
234 const ComponentVector& targetFug)
235 {
236 // reset jacobian
237 J = 0;
238
239 Scalar absError = 0;
240 // calculate the defect (deviation of the current fugacities
241 // from the target fugacities)
242 for (unsigned i = 0; i < numComponents; ++ i) {
243 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
244 paramCache,
245 phaseIdx,
246 i);
247 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
248 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
249
250 defect[i] = targetFug[i] - f;
251 absError = std::max(absError, std::abs(scalarValue(defect[i])));
252 }
253
254 // assemble jacobian matrix of the constraints for the composition
255 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
256 for (unsigned i = 0; i < numComponents; ++ i) {
258 // approximately calculate partial derivatives of the
259 // fugacity defect of all components in regard to the mole
260 // fraction of the i-th component. This is done via
261 // forward differences
262
263 // deviate the mole fraction of the i-th component
264 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
265 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
266 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
267
268 // compute new defect and derivative for all component
269 // fugacities
270 for (unsigned j = 0; j < numComponents; ++j) {
271 // compute the j-th component's fugacity coefficient ...
272 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
273 paramCache,
274 phaseIdx,
275 j);
276 // ... and its fugacity ...
277 const Evaluation& f =
278 phi *
279 fluidState.pressure(phaseIdx) *
280 fluidState.moleFraction(phaseIdx, j);
281 // as well as the defect for this fugacity
282 const Evaluation& defJPlusEps = targetFug[j] - f;
283
284 // use forward differences to calculate the defect's
285 // derivative
286 J[j][i] = (defJPlusEps - defect[j])/eps;
287 }
288
289 // reset composition to original value
290 fluidState.setMoleFraction(phaseIdx, i, xI);
291 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
292
293 // end forward differences
295 }
296
297 return absError;
298 }
299
300 template <class FluidState>
301 static Scalar update_(FluidState& fluidState,
302 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
303 Dune::FieldVector<Evaluation, numComponents>& x,
304 Dune::FieldVector<Evaluation, numComponents>& /*b*/,
305 unsigned phaseIdx,
306 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
307 {
308 // store original composition and calculate relative error
309 Dune::FieldVector<Evaluation, numComponents> origComp;
310 Scalar relError = 0;
311 Evaluation sumDelta = 0.0;
312 Evaluation sumx = 0.0;
313 for (unsigned i = 0; i < numComponents; ++i) {
314 origComp[i] = fluidState.moleFraction(phaseIdx, i);
315 relError = std::max(relError, std::abs(scalarValue(x[i])));
316
317 sumx += abs(fluidState.moleFraction(phaseIdx, i));
318 sumDelta += abs(x[i]);
319 }
320
321 // chop update to at most 20% change in composition
322 const Scalar maxDelta = 0.2;
323 if (sumDelta > maxDelta)
324 x /= (sumDelta/maxDelta);
325
326 // change composition
327 for (unsigned i = 0; i < numComponents; ++i) {
328 Evaluation newx = origComp[i] - x[i];
329 // only allow negative mole fractions if the target fugacity is negative
330 if (targetFug[i] > 0)
331 newx = max(0.0, newx);
332 // only allow positive mole fractions if the target fugacity is positive
333 else if (targetFug[i] < 0)
334 newx = min(0.0, newx);
335 // if the target fugacity is zero, the mole fraction must also be zero
336 else
337 newx = 0;
338
339 fluidState.setMoleFraction(phaseIdx, i, newx);
340 }
341
342 paramCache.updateComposition(fluidState, phaseIdx);
343
344 return relError;
345 }
346
347 template <class FluidState>
348 static Scalar calculateDefect_(const FluidState& params,
349 unsigned phaseIdx,
350 const ComponentVector& targetFug)
351 {
352 Scalar result = 0.0;
353 for (unsigned i = 0; i < numComponents; ++i) {
354 // sum of the fugacity defect weighted by the inverse
355 // fugacity coefficient
356 result += std::abs(
357 (targetFug[i] - params.fugacity(phaseIdx, i))
358 /
359 params.fugacityCoefficient(phaseIdx, i) );
360 };
361 return result;
362 }
363}; // namespace Opm
364
365} // end namespace Opm
366
367#endif
Provides the OPM specific exception classes.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Some templates to wrap the valgrind client request macros.
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition CompositionFromFugacities.hpp:50
static void solve(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache, unsigned phaseIdx, const ComponentVector &targetFug)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition CompositionFromFugacities.hpp:83
static void guessInitial(FluidState &fluidState, unsigned phaseIdx, const ComponentVector &)
Guess an initial value for the composition of the phase.
Definition CompositionFromFugacities.hpp:60
Definition Exceptions.hpp:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30