My Project
Loading...
Searching...
No Matches
ImmiscibleFlash.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_IMMISCIBLE_FLASH_HPP
28#define OPM_IMMISCIBLE_FLASH_HPP
29
31
37
38#include <dune/common/fvector.hh>
39#include <dune/common/fmatrix.hh>
40
41#include <limits>
42
43namespace Opm {
44
71template <class Scalar, class FluidSystem>
73{
74 static const int numPhases = FluidSystem::numPhases;
75 static const int numComponents = FluidSystem::numComponents;
76 static_assert(numPhases == numComponents,
77 "Immiscibility assumes that the number of phases"
78 " is equal to the number of components");
79
80 enum {
81 p0PvIdx = 0,
82 S0PvIdx = 1
83 };
84
85 static const int numEq = numPhases;
86
87public:
91 template <class FluidState, class Evaluation = typename FluidState::Scalar>
92 static void guessInitial(FluidState& fluidState,
93 const Dune::FieldVector<Evaluation, numComponents>& /*globalMolarities*/)
94 {
95 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
96 // pressure. use 1 bar as initial guess
97 fluidState.setPressure(phaseIdx, 1e5);
98
99 // saturation. assume all fluids to be equally distributed
100 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
101 }
102 }
103
110 template <class MaterialLaw, class FluidState>
111 static void solve(FluidState& fluidState,
112 const typename MaterialLaw::Params& matParams,
113 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
114 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
115 Scalar tolerance = -1)
116 {
117 typedef typename FluidState::Scalar InputEval;
118
120 // Check if all fluid phases are incompressible
122 bool allIncompressible = true;
123 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
124 if (FluidSystem::isCompressible(phaseIdx)) {
125 allIncompressible = false;
126 break;
127 }
128 }
129
130 if (allIncompressible) {
131 // yes, all fluid phases are incompressible. in this case the flash solver
132 // can only determine the saturations, not the pressures. (but this
133 // determination is much simpler than a full flash calculation.)
134 paramCache.updateAll(fluidState);
135 solveAllIncompressible_(fluidState, paramCache, globalMolarities);
136 return;
137 }
138
139 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
140 typedef Dune::FieldVector<InputEval, numEq> Vector;
141
143 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
144 typedef ImmiscibleFluidState<FlashEval, FluidSystem> FlashFluidState;
145
146 if (tolerance <= 0)
147 tolerance = std::min<Scalar>(1e-5,
148 1e8*std::numeric_limits<Scalar>::epsilon());
149
150 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
151 flashParamCache.assignPersistentData(paramCache);
152
154 // Newton method
156
157 // Jacobian matrix
158 Matrix J;
159 // solution, i.e. phase composition
160 Vector deltaX;
161 // right hand side
162 Vector b;
163
164 Valgrind::SetUndefined(J);
165 Valgrind::SetUndefined(deltaX);
166 Valgrind::SetUndefined(b);
167
168 FlashFluidState flashFluidState;
169 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
170
171 // copy the global molarities to a vector of evaluations. Remember that the
172 // global molarities are constants. (but we need to copy them to a vector of
173 // FlashEvals anyway in order to avoid getting into hell's kitchen.)
174 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
175 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
176 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
177
178 FlashDefectVector defect;
179 const unsigned nMax = 50; // <- maximum number of newton iterations
180 for (unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
181 // calculate Jacobian matrix and right hand side
182 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
183 Valgrind::CheckDefined(defect);
184
185 // create field matrices and vectors out of the evaluation vector to solve
186 // the linear system of equations.
187 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
188 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
189 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
190
191 b[eqIdx] = defect[eqIdx].value();
192 }
193 Valgrind::CheckDefined(J);
194 Valgrind::CheckDefined(b);
195
196 // Solve J*x = b
197 deltaX = 0;
198
199 try { J.solve(deltaX, b); }
200 catch (const Dune::FMatrixError& e) {
201 throw NumericalProblem(e.what());
202 }
203 Valgrind::CheckDefined(deltaX);
204
205 // update the fluid quantities.
206 Scalar relError = update_<MaterialLaw>(flashFluidState, flashParamCache, matParams, deltaX);
207
208 if (relError < tolerance) {
209 assignOutputFluidState_(flashFluidState, fluidState);
210 return;
211 }
212 }
213
214 std::string msg = "ImmiscibleFlash solver failed: "
215 "{c_alpha^kappa} = {";
216 for (const auto& v : globalMolarities)
217 msg += " " + std::to_string(v);
218 msg += " }, T = ";
219 msg += std::to_string(fluidState.temperature(/*phaseIdx=*/0));
220 throw NumericalProblem(msg);
221 }
222
223
224protected:
225 template <class MaterialLaw, class InputFluidState, class FlashFluidState>
226 static void assignFlashFluidState_(const InputFluidState& inputFluidState,
227 FlashFluidState& flashFluidState,
228 const typename MaterialLaw::Params& matParams,
229 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
230 {
231 typedef typename FlashFluidState::Scalar FlashEval;
232
233 // copy the temperature: even though the model which uses the flash solver might
234 // be non-isothermal, the flash solver does not consider energy. (it could be
235 // modified to do so relatively easily, but it would come at increased
236 // computational cost and normally temperature instead of "total internal energy
237 // of the fluids" is specified.)
238 flashFluidState.setTemperature(inputFluidState.temperature(/*phaseIdx=*/0));
239
240 // copy the saturations: the first N-1 phases are primary variables, the last one
241 // is one minus the sum of the former.
242 FlashEval Slast = 1.0;
243 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
244 FlashEval S = inputFluidState.saturation(phaseIdx);
245 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
246
247 Slast -= S;
248
249 flashFluidState.setSaturation(phaseIdx, S);
250 }
251 flashFluidState.setSaturation(numPhases - 1, Slast);
252
253 // copy the pressures: the first pressure is the first primary variable, the
254 // remaining ones are given as p_beta = p_alpha + p_calpha,beta
255 FlashEval p0 = inputFluidState.pressure(0);
256 p0.setDerivative(p0PvIdx, 1.0);
257
258 std::array<FlashEval, numPhases> pc;
259 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
260 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
261 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
262
263 flashParamCache.updateAll(flashFluidState);
264
265 // compute the density of each phase.
266 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
267 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
268 flashFluidState.setDensity(phaseIdx, rho);
269 }
270 }
271
272 template <class FlashFluidState, class OutputFluidState>
273 static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
274 OutputFluidState& outputFluidState)
275 {
276 outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value());
277
278 // copy the saturations, pressures and densities
279 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
280 const auto& S = flashFluidState.saturation(phaseIdx).value();
281 outputFluidState.setSaturation(phaseIdx, S);
282
283 const auto& p = flashFluidState.pressure(phaseIdx).value();
284 outputFluidState.setPressure(phaseIdx, p);
285
286 const auto& rho = flashFluidState.density(phaseIdx).value();
287 outputFluidState.setDensity(phaseIdx, rho);
288 }
289 }
290
291 template <class FluidState>
292 static void solveAllIncompressible_(FluidState& fluidState,
293 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
294 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities)
295 {
296 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
297 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
298 fluidState.setDensity(phaseIdx, rho);
299
300 Scalar saturation =
301 globalMolarities[/*compIdx=*/phaseIdx]
302 / fluidState.molarDensity(phaseIdx);
303 fluidState.setSaturation(phaseIdx, saturation);
304 }
305 }
306
307 template <class FluidState, class FlashDefectVector, class FlashComponentVector>
308 static void evalDefect_(FlashDefectVector& b,
309 const FluidState& fluidState,
310 const FlashComponentVector& globalMolarities)
311 {
312 // global molarities are given
313 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
314 b[phaseIdx] =
315 fluidState.saturation(phaseIdx)
316 * fluidState.molarity(phaseIdx, /*compIdx=*/phaseIdx);
317 b[phaseIdx] -= globalMolarities[/*compIdx=*/phaseIdx];
318 }
319 }
320
321 template <class MaterialLaw, class FlashFluidState, class EvalVector>
322 static Scalar update_(FlashFluidState& fluidState,
323 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
324 const typename MaterialLaw::Params& matParams,
325 const EvalVector& deltaX)
326 {
327 // note that it is possible that FlashEval::Scalar is an Evaluation itself
328 typedef typename FlashFluidState::Scalar FlashEval;
329 typedef MathToolbox<FlashEval> FlashEvalToolbox;
330
331 typedef typename FlashEvalToolbox::ValueType InnerEval;
332
333#ifndef NDEBUG
334 // make sure we don't swallow non-finite update vectors
335 assert(deltaX.dimension == numEq);
336 for (unsigned i = 0; i < numEq; ++i)
337 assert(std::isfinite(scalarValue(deltaX[i])));
338#endif
339
340 Scalar relError = 0;
341 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
342 FlashEval tmp = getQuantity_(fluidState, pvIdx);
343 InnerEval delta = deltaX[pvIdx];
344
345 relError = std::max(relError,
346 std::abs(scalarValue(delta))
347 * quantityWeight_(fluidState, pvIdx));
348
349 if (isSaturationIdx_(pvIdx)) {
350 // dampen to at most 20% change in saturation per
351 // iteration
352 delta = min(0.25, max(-0.25, delta));
353 }
354 else if (isPressureIdx_(pvIdx)) {
355 // dampen to at most 30% change in pressure per
356 // iteration
357 delta = min(0.5*fluidState.pressure(0).value(),
358 max(-0.50*fluidState.pressure(0).value(), delta));
359 }
360
361 tmp -= delta;
362 setQuantity_(fluidState, pvIdx, tmp);
363 }
364
365 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
366
367 return relError;
368 }
369
370 template <class MaterialLaw, class FlashFluidState>
371 static void completeFluidState_(FlashFluidState& flashFluidState,
372 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
373 const typename MaterialLaw::Params& matParams)
374 {
375 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
376
377 typedef typename FlashFluidState::Scalar FlashEval;
378
379 // calculate the saturation of the last phase as a function of
380 // the other saturations
381 FlashEval sumSat = 0.0;
382 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
383 sumSat += flashFluidState.saturation(phaseIdx);
384 flashFluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
385
386 // update the pressures using the material law (saturations
387 // and first pressure are already set because it is implicitly
388 // solved for.)
389 std::array<FlashEval, numPhases> pC;
390 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
391 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
392 flashFluidState.setPressure(phaseIdx,
393 flashFluidState.pressure(0)
394 + (pC[phaseIdx] - pC[0]));
395
396 // update the parameter cache
397 paramCache.updateAll(flashFluidState, /*except=*/ParamCache::Temperature|ParamCache::Composition);
398
399 // update all densities
400 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
401 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
402 flashFluidState.setDensity(phaseIdx, rho);
403 }
404 }
405
406 static bool isPressureIdx_(unsigned pvIdx)
407 { return pvIdx == 0; }
408
409 static bool isSaturationIdx_(unsigned pvIdx)
410 { return 1 <= pvIdx && pvIdx < numPhases; }
411
412 // retrieves a quantity from the fluid state
413 template <class FluidState>
414 static const typename FluidState::Scalar& getQuantity_(const FluidState& fluidState, unsigned pvIdx)
415 {
416 assert(pvIdx < numEq);
417
418 // first pressure
419 if (pvIdx < 1) {
420 unsigned phaseIdx = 0;
421 return fluidState.pressure(phaseIdx);
422 }
423 // saturations
424 else {
425 assert(pvIdx < numPhases);
426 unsigned phaseIdx = pvIdx - 1;
427 return fluidState.saturation(phaseIdx);
428 }
429 }
430
431 // set a quantity in the fluid state
432 template <class FluidState>
433 static void setQuantity_(FluidState& fluidState,
434 unsigned pvIdx,
435 const typename FluidState::Scalar& value)
436 {
437 assert(pvIdx < numEq);
438
439 // first pressure
440 if (pvIdx < 1) {
441 unsigned phaseIdx = 0;
442 fluidState.setPressure(phaseIdx, value);
443 }
444 // saturations
445 else {
446 assert(pvIdx < numPhases);
447 unsigned phaseIdx = pvIdx - 1;
448 fluidState.setSaturation(phaseIdx, value);
449 }
450 }
451
452 template <class FluidState>
453 static Scalar quantityWeight_(const FluidState& /*fluidState*/, unsigned pvIdx)
454 {
455 // first pressure
456 if (pvIdx < 1)
457 return 1e-8;
458 // first M - 1 saturations
459 else {
460 assert(pvIdx < numPhases);
461 return 1.0;
462 }
463 }
464};
465
466} // namespace Opm
467
468#endif
Representation of an evaluation of a function and its derivatives w.r.t.
Provides the OPM specific exception classes.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Some templates to wrap the valgrind client request macros.
Represents a function evaluation and its derivatives w.r.t.
Definition Evaluation.hpp:57
Determines the pressures and saturations of all fluid phases given the total mass of all components.
Definition ImmiscibleFlash.hpp:73
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &)
Guess initial values for all quantities.
Definition ImmiscibleFlash.hpp:92
static void solve(FluidState &fluidState, const typename MaterialLaw::Params &matParams, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache, const Dune::FieldVector< typename FluidState::Scalar, numComponents > &globalMolarities, Scalar tolerance=-1)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition ImmiscibleFlash.hpp:111
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition ImmiscibleFluidState.hpp:49
Definition Exceptions.hpp:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30