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