My Project
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 <limits>
39
40namespace Opm {
41
46template <class Scalar, class FluidSystem, class Evaluation = Scalar>
48{
49 enum { numComponents = FluidSystem::numComponents };
50
51public:
52 typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
53
57 template <class FluidState>
58 static void guessInitial(FluidState& fluidState,
59 unsigned phaseIdx,
60 const ComponentVector& /*fugVec*/)
61 {
62 if (FluidSystem::isIdealMixture(phaseIdx))
63 return;
64
65 // Pure component fugacities
66 for (unsigned i = 0; i < numComponents; ++ i) {
67 //std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
68 fluidState.setMoleFraction(phaseIdx,
69 i,
70 1.0/numComponents);
71 }
72 }
73
80 template <class FluidState>
81 static void solve(FluidState& fluidState,
82 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
83 unsigned phaseIdx,
84 const ComponentVector& targetFug)
85 {
86 // use a much more efficient method in case the phase is an
87 // ideal mixture
88 if (FluidSystem::isIdealMixture(phaseIdx)) {
89 solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
90 return;
91 }
92
93 // save initial composition in case something goes wrong
94 Dune::FieldVector<Evaluation, numComponents> xInit;
95 for (unsigned i = 0; i < numComponents; ++i) {
96 xInit[i] = fluidState.moleFraction(phaseIdx, i);
97 }
98
100 // Newton method
102
103 // Jacobian matrix
104 Dune::FieldMatrix<Evaluation, numComponents, numComponents> J;
105 // solution, i.e. phase composition
106 Dune::FieldVector<Evaluation, numComponents> x;
107 // right hand side
108 Dune::FieldVector<Evaluation, numComponents> b;
109
110 paramCache.updatePhase(fluidState, phaseIdx);
111
112 // maximum number of iterations
113 const int nMax = 25;
114 for (int nIdx = 0; nIdx < nMax; ++nIdx) {
115 // calculate Jacobian matrix and right hand side
116 linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug);
117 Valgrind::CheckDefined(J);
118 Valgrind::CheckDefined(b);
119
120 /*
121 std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
122 for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
123 std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
124 std::cout << "\n";
125 std::cout << FluidSystem::phaseName(phaseIdx) << "Phase phi: ";
126 for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
127 std::cout << fluidState.fugacityCoefficient(phaseIdx, i) << " ";
128 std::cout << "\n";
129 */
130
131 // Solve J*x = b
132 x = 0.0;
133 try { J.solve(x, b); }
134 catch (const Dune::FMatrixError& e)
135 { throw NumericalProblem(e.what()); }
136
137 //std::cout << "original delta: " << x << "\n";
138
139 Valgrind::CheckDefined(x);
140
141 /*
142 std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
143 for (unsigned i = 0; i < FluidSystem::numComponents; ++i)
144 std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
145 std::cout << "\n";
146 std::cout << "J: " << J << "\n";
147 std::cout << "rho: " << fluidState.density(phaseIdx) << "\n";
148 std::cout << "delta: " << x << "\n";
149 std::cout << "defect: " << b << "\n";
150
151 std::cout << "J: " << J << "\n";
152
153 std::cout << "---------------------------\n";
154 */
155
156 // update the fluid composition. b is also used to store
157 // the defect for the next iteration.
158 Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
159
160 if (relError < 1e-9) {
161 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
162 fluidState.setDensity(phaseIdx, rho);
163
164 //std::cout << "num iterations: " << nIdx << "\n";
165 return;
166 }
167 }
168
169 auto cast = [](const auto d)
170 {
171#if HAVE_QUAD
172 if constexpr (std::is_same_v<decltype(d), const quad>)
173 return static_cast<double>(d);
174 else
175#endif
176 return d;
177 };
178
179 std::string msg =
180 std::string("Calculating the ")
181 + FluidSystem::phaseName(phaseIdx)
182 + "Phase composition failed. Initial {x} = {";
183 for (const auto& v : xInit)
184 msg += " " + std::to_string(cast(getValue(v)));
185 msg += " }, {fug_t} = {";
186 for (const auto& v : targetFug)
187 msg += " " + std::to_string(cast(getValue(v)));
188 msg += " }, p = " + std::to_string(cast(getValue(fluidState.pressure(phaseIdx))))
189 + ", T = " + std::to_string(cast(getValue(fluidState.temperature(phaseIdx))));
190 throw NumericalProblem(msg);
191 }
192
193
194protected:
195 // update the phase composition in case the phase is an ideal
196 // mixture, i.e. the component's fugacity coefficients are
197 // independent of the phase's composition.
198 template <class FluidState>
199 static void solveIdealMix_(FluidState& fluidState,
200 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
201 unsigned phaseIdx,
202 const ComponentVector& fugacities)
203 {
204 for (unsigned i = 0; i < numComponents; ++ i) {
205 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
206 paramCache,
207 phaseIdx,
208 i);
209 const Evaluation& gamma = phi * fluidState.pressure(phaseIdx);
210 Valgrind::CheckDefined(phi);
211 Valgrind::CheckDefined(gamma);
212 Valgrind::CheckDefined(fugacities[i]);
213 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
214 fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
215 };
216
217 paramCache.updatePhase(fluidState, phaseIdx);
218
219 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
220 fluidState.setDensity(phaseIdx, rho);
221 return;
222 }
223
224 template <class FluidState>
225 static Scalar linearize_(Dune::FieldMatrix<Evaluation, numComponents, numComponents>& J,
226 Dune::FieldVector<Evaluation, numComponents>& defect,
227 FluidState& fluidState,
228 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
229 unsigned phaseIdx,
230 const ComponentVector& targetFug)
231 {
232 // reset jacobian
233 J = 0;
234
235 Scalar absError = 0;
236 // calculate the defect (deviation of the current fugacities
237 // from the target fugacities)
238 for (unsigned i = 0; i < numComponents; ++ i) {
239 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
240 paramCache,
241 phaseIdx,
242 i);
243 const Evaluation& f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
244 fluidState.setFugacityCoefficient(phaseIdx, i, phi);
245
246 defect[i] = targetFug[i] - f;
247 absError = std::max(absError, std::abs(scalarValue(defect[i])));
248 }
249
250 // assemble jacobian matrix of the constraints for the composition
251 static const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e6;
252 for (unsigned i = 0; i < numComponents; ++ i) {
254 // approximately calculate partial derivatives of the
255 // fugacity defect of all components in regard to the mole
256 // fraction of the i-th component. This is done via
257 // forward differences
258
259 // deviate the mole fraction of the i-th component
260 Evaluation xI = fluidState.moleFraction(phaseIdx, i);
261 fluidState.setMoleFraction(phaseIdx, i, xI + eps);
262 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
263
264 // compute new defect and derivative for all component
265 // fugacities
266 for (unsigned j = 0; j < numComponents; ++j) {
267 // compute the j-th component's fugacity coefficient ...
268 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState,
269 paramCache,
270 phaseIdx,
271 j);
272 // ... and its fugacity ...
273 const Evaluation& f =
274 phi *
275 fluidState.pressure(phaseIdx) *
276 fluidState.moleFraction(phaseIdx, j);
277 // as well as the defect for this fugacity
278 const Evaluation& defJPlusEps = targetFug[j] - f;
279
280 // use forward differences to calculate the defect's
281 // derivative
282 J[j][i] = (defJPlusEps - defect[j])/eps;
283 }
284
285 // reset composition to original value
286 fluidState.setMoleFraction(phaseIdx, i, xI);
287 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, i);
288
289 // end forward differences
291 }
292
293 return absError;
294 }
295
296 template <class FluidState>
297 static Scalar update_(FluidState& fluidState,
298 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
299 Dune::FieldVector<Evaluation, numComponents>& x,
300 Dune::FieldVector<Evaluation, numComponents>& /*b*/,
301 unsigned phaseIdx,
302 const Dune::FieldVector<Evaluation, numComponents>& targetFug)
303 {
304 // store original composition and calculate relative error
305 Dune::FieldVector<Evaluation, numComponents> origComp;
306 Scalar relError = 0;
307 Evaluation sumDelta = 0.0;
308 Evaluation sumx = 0.0;
309 for (unsigned i = 0; i < numComponents; ++i) {
310 origComp[i] = fluidState.moleFraction(phaseIdx, i);
311 relError = std::max(relError, std::abs(scalarValue(x[i])));
312
313 sumx += abs(fluidState.moleFraction(phaseIdx, i));
314 sumDelta += abs(x[i]);
315 }
316
317 // chop update to at most 20% change in composition
318 const Scalar maxDelta = 0.2;
319 if (sumDelta > maxDelta)
320 x /= (sumDelta/maxDelta);
321
322 // change composition
323 for (unsigned i = 0; i < numComponents; ++i) {
324 Evaluation newx = origComp[i] - x[i];
325 // only allow negative mole fractions if the target fugacity is negative
326 if (targetFug[i] > 0)
327 newx = max(0.0, newx);
328 // only allow positive mole fractions if the target fugacity is positive
329 else if (targetFug[i] < 0)
330 newx = min(0.0, newx);
331 // if the target fugacity is zero, the mole fraction must also be zero
332 else
333 newx = 0;
334
335 fluidState.setMoleFraction(phaseIdx, i, newx);
336 }
337
338 paramCache.updateComposition(fluidState, phaseIdx);
339
340 return relError;
341 }
342
343 template <class FluidState>
344 static Scalar calculateDefect_(const FluidState& params,
345 unsigned phaseIdx,
346 const ComponentVector& targetFug)
347 {
348 Scalar result = 0.0;
349 for (unsigned i = 0; i < numComponents; ++i) {
350 // sum of the fugacity defect weighted by the inverse
351 // fugacity coefficient
352 result += std::abs(
353 (targetFug[i] - params.fugacity(phaseIdx, i))
354 /
355 params.fugacityCoefficient(phaseIdx, i) );
356 };
357 return result;
358 }
359}; // namespace Opm
360
361} // end namespace Opm
362
363#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:48
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:81
static void guessInitial(FluidState &fluidState, unsigned phaseIdx, const ComponentVector &)
Guess an initial value for the composition of the phase.
Definition: CompositionFromFugacities.hpp:58
Definition: Exceptions.hpp:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition: Exceptions.hpp:30