My Project
Loading...
Searching...
No Matches
NcpFlash.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_NCP_FLASH_HPP
28#define OPM_NCP_FLASH_HPP
29
31
39
40#include <dune/common/fvector.hh>
41#include <dune/common/fmatrix.hh>
42
43#include <limits>
44
45namespace Opm {
46
86template <class Scalar, class FluidSystem>
88{
89 enum { numPhases = FluidSystem::numPhases };
90 enum { numComponents = FluidSystem::numComponents };
91
92 enum {
93 p0PvIdx = 0,
94 S0PvIdx = 1,
95 x00PvIdx = S0PvIdx + numPhases - 1
96 };
97
98 static const int numEq = numPhases*(numComponents + 1);
99
100public:
104 template <class FluidState, class Evaluation = typename FluidState::Scalar>
105 static void guessInitial(FluidState& fluidState,
106 const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
107 {
108 // the sum of all molarities
109 Evaluation sumMoles = 0;
110 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
111 sumMoles += globalMolarities[compIdx];
112
113 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
114 // composition
115 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
116 fluidState.setMoleFraction(phaseIdx,
117 compIdx,
118 globalMolarities[compIdx]/sumMoles);
119
120 // pressure. use atmospheric pressure as initial guess
121 fluidState.setPressure(phaseIdx, 1.0135e5);
122
123 // saturation. assume all fluids to be equally distributed
124 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
125 }
126
127 // set the fugacity coefficients of all components in all phases
128 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
129 paramCache.updateAll(fluidState);
130 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
131 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
132 const typename FluidState::Scalar phi =
133 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
134 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
135 }
136 }
137 }
138
145 template <class MaterialLaw, class FluidState>
146 static void solve(FluidState& fluidState,
147 const typename MaterialLaw::Params& matParams,
148 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
149 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
150 Scalar tolerance = -1.0)
151 {
152 typedef typename FluidState::Scalar InputEval;
153
154 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
155 typedef Dune::FieldVector<InputEval, numEq> Vector;
156
157 typedef DenseAd::Evaluation</*Scalar=*/InputEval,
158 /*numDerivs=*/numEq> FlashEval;
159
160 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
161 typedef CompositionalFluidState<FlashEval, FluidSystem, /*energy=*/false> FlashFluidState;
162
163 if (tolerance <= 0)
164 tolerance = std::min<Scalar>(1e-3,
165 1e8*std::numeric_limits<Scalar>::epsilon());
166
167 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
168 flashParamCache.assignPersistentData(paramCache);
169
171 // Newton method
173
174 // Jacobian matrix
175 Matrix J;
176 // solution, i.e. phase composition
177 Vector deltaX;
178 // right hand side
179 Vector b;
180
181 Valgrind::SetUndefined(J);
182 Valgrind::SetUndefined(deltaX);
183 Valgrind::SetUndefined(b);
184
185 FlashFluidState flashFluidState;
186 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
187
188 // copy the global molarities to a vector of evaluations. Remember that the
189 // global molarities are constants. (but we need to copy them to a vector of
190 // FlashEvals anyway in order to avoid getting into hell's kitchen.)
191 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
192 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
193 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
194
195 FlashDefectVector defect;
196 const unsigned nMax = 50; // <- maximum number of newton iterations
197 for (unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
198 // calculate the defect of the flash equations and their derivatives
199 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
200 Valgrind::CheckDefined(defect);
201
202 // create field matrices and vectors out of the evaluation vector to solve
203 // the linear system of equations.
204 for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
205 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
206 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
207
208 b[eqIdx] = defect[eqIdx].value();
209 }
210 Valgrind::CheckDefined(J);
211 Valgrind::CheckDefined(b);
212
213 // Solve J*x = b
214 deltaX = 0.0;
215 try { J.solve(deltaX, b); }
216 catch (const Dune::FMatrixError& e) {
217 throw NumericalProblem(e.what());
218 }
219 Valgrind::CheckDefined(deltaX);
220
221 // update the fluid quantities.
222 Scalar relError = update_<MaterialLaw>(flashFluidState, matParams, flashParamCache, deltaX);
223
224 if (relError < tolerance) {
225 assignOutputFluidState_(flashFluidState, fluidState);
226 return;
227 }
228 }
229
230 auto cast = [](const auto d)
231 {
232#if HAVE_QUAD
233 if constexpr (std::is_same_v<decltype(d), const quad>)
234 return static_cast<double>(d);
235 else
236#endif
237 return d;
238 };
239
240 std::string msg = "NcpFlash solver failed: "
241 "{c_alpha^kappa} = {";
242 for (const auto& v : globalMolarities)
243 msg += " " + std::to_string(cast(getValue(v)));
244 msg += " }, T = ";
245 msg += std::to_string(cast(getValue(fluidState.temperature(/*phaseIdx=*/0))));
246 throw NumericalProblem(msg);
247 }
248
256 template <class FluidState, class ComponentVector>
257 static void solve(FluidState& fluidState,
258 const ComponentVector& globalMolarities,
259 Scalar tolerance = 0.0)
260 {
261 typedef NullMaterialTraits<Scalar, numPhases> MaterialTraits;
262 typedef NullMaterial<MaterialTraits> MaterialLaw;
263 typedef typename MaterialLaw::Params MaterialLawParams;
264
265 MaterialLawParams matParams;
266 solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
267 }
268
269
270protected:
271 template <class MaterialLaw, class InputFluidState, class FlashFluidState>
272 static void assignFlashFluidState_(const InputFluidState& inputFluidState,
273 FlashFluidState& flashFluidState,
274 const typename MaterialLaw::Params& matParams,
275 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
276 {
277 typedef typename FlashFluidState::Scalar FlashEval;
278
279 // copy the temperature: even though the model which uses the flash solver might
280 // be non-isothermal, the flash solver does not consider energy. (it could be
281 // modified to do so relatively easily, but it would come at increased
282 // computational cost and normally temperature instead of "total internal energy
283 // of the fluids" is specified.)
284 flashFluidState.setTemperature(inputFluidState.temperature(/*phaseIdx=*/0));
285
286 // copy the saturations: the first N-1 phases are primary variables, the last one
287 // is one minus the sum of the former.
288 FlashEval Slast = 1.0;
289 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
290 FlashEval S = inputFluidState.saturation(phaseIdx);
291 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
292
293 Slast -= S;
294
295 flashFluidState.setSaturation(phaseIdx, S);
296 }
297 flashFluidState.setSaturation(numPhases - 1, Slast);
298
299 // copy the pressures: the first pressure is the first primary variable, the
300 // remaining ones are given as p_beta = p_alpha + p_calpha,beta
301 FlashEval p0 = inputFluidState.pressure(0);
302 p0.setDerivative(p0PvIdx, 1.0);
303
304 std::array<FlashEval, numPhases> pc;
305 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
306 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
307 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
308
309 // copy the mole fractions: all of them are primary variables
310 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
311 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
312 FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
313 x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
314 flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
315 }
316 }
317
318 flashParamCache.updateAll(flashFluidState);
319
320 // compute the density of each phase and the fugacity coefficient of each
321 // component in each phase.
322 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
323 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
324 flashFluidState.setDensity(phaseIdx, rho);
325
326 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
327 const FlashEval& fugCoeff = FluidSystem::fugacityCoefficient(flashFluidState, flashParamCache, phaseIdx, compIdx);
328 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
329 }
330 }
331 }
332
333 template <class FlashFluidState, class OutputFluidState>
334 static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
335 OutputFluidState& outputFluidState)
336 {
337 outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value());
338
339 // copy the saturations, pressures and densities
340 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
341 const auto& S = flashFluidState.saturation(phaseIdx).value();
342 outputFluidState.setSaturation(phaseIdx, S);
343
344 const auto& p = flashFluidState.pressure(phaseIdx).value();
345 outputFluidState.setPressure(phaseIdx, p);
346
347 const auto& rho = flashFluidState.density(phaseIdx).value();
348 outputFluidState.setDensity(phaseIdx, rho);
349 }
350
351 // copy the mole fractions and fugacity coefficients
352 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
353 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
354 const auto& moleFrac =
355 flashFluidState.moleFraction(phaseIdx, compIdx).value();
356 outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
357
358 const auto& fugCoeff =
359 flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
360 outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
361 }
362 }
363 }
364
365 template <class FlashFluidState, class FlashDefectVector, class FlashComponentVector>
366 static void evalDefect_(FlashDefectVector& b,
367 const FlashFluidState& fluidState,
368 const FlashComponentVector& globalMolarities)
369 {
370 typedef typename FlashFluidState::Scalar FlashEval;
371
372 unsigned eqIdx = 0;
373
374 // fugacity of any component must be equal in all phases
375 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
376 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
377 b[eqIdx] =
378 fluidState.fugacity(/*phaseIdx=*/0, compIdx)
379 - fluidState.fugacity(phaseIdx, compIdx);
380 ++eqIdx;
381 }
382 }
383 assert(eqIdx == numComponents*(numPhases - 1));
384
385 // the fact saturations must sum up to 1 is included implicitly and also,
386 // capillary pressures are treated implicitly!
387
388 // global molarities are given
389 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
390 b[eqIdx] = 0.0;
391 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
392 b[eqIdx] +=
393 fluidState.saturation(phaseIdx)
394 * fluidState.molarity(phaseIdx, compIdx);
395 }
396
397 b[eqIdx] -= globalMolarities[compIdx];
398 ++eqIdx;
399 }
400
401 // model assumptions (-> non-linear complementarity functions) must be adhered
402 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
403 FlashEval oneMinusSumMoleFrac = 1.0;
404 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
405 oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
406
407 if (oneMinusSumMoleFrac > fluidState.saturation(phaseIdx))
408 b[eqIdx] = fluidState.saturation(phaseIdx);
409 else
410 b[eqIdx] = oneMinusSumMoleFrac;
411
412 ++eqIdx;
413 }
414 }
415
416 template <class MaterialLaw, class FlashFluidState, class EvalVector>
417 static Scalar update_(FlashFluidState& fluidState,
418 const typename MaterialLaw::Params& matParams,
419 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
420 const EvalVector& deltaX)
421 {
422 // note that it is possible that FlashEval::Scalar is an Evaluation itself
423 typedef typename FlashFluidState::Scalar FlashEval;
424 typedef typename FlashEval::ValueType InnerEval;
425
426#ifndef NDEBUG
427 // make sure we don't swallow non-finite update vectors
428 assert(deltaX.dimension == numEq);
429 for (unsigned i = 0; i < numEq; ++i)
430 assert(std::isfinite(scalarValue(deltaX[i])));
431#endif
432
433 Scalar relError = 0;
434 for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
435 FlashEval tmp = getQuantity_(fluidState, pvIdx);
436 InnerEval delta = deltaX[pvIdx];
437
438 relError = std::max(relError,
439 std::abs(scalarValue(delta))
440 * quantityWeight_(fluidState, pvIdx));
441
442 if (isSaturationIdx_(pvIdx)) {
443 // dampen to at most 25% change in saturation per iteration
444 delta = min(0.25, max(-0.25, delta));
445 }
446 else if (isMoleFracIdx_(pvIdx)) {
447 // dampen to at most 20% change in mole fraction per iteration
448 delta = min(0.20, max(-0.20, delta));
449 }
450 else if (isPressureIdx_(pvIdx)) {
451 // dampen to at most 50% change in pressure per iteration
452 delta = min(0.5*fluidState.pressure(0).value(),
453 max(-0.5*fluidState.pressure(0).value(),
454 delta));
455 }
456
457 tmp -= delta;
458 setQuantity_(fluidState, pvIdx, tmp);
459 }
460
461 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
462
463 return relError;
464 }
465
466 template <class MaterialLaw, class FlashFluidState>
467 static void completeFluidState_(FlashFluidState& flashFluidState,
468 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
469 const typename MaterialLaw::Params& matParams)
470 {
471 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
472
473 typedef typename FlashFluidState::Scalar FlashEval;
474
475 // calculate the saturation of the last phase as a function of
476 // the other saturations
477 FlashEval sumSat = 0.0;
478 for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
479 sumSat += flashFluidState.saturation(phaseIdx);
480 flashFluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);
481
482 // update the pressures using the material law (saturations
483 // and first pressure are already set because it is implicitly
484 // solved for.)
485 std::array<FlashEval, numPhases> pC;
486 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
487 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
488 flashFluidState.setPressure(phaseIdx,
489 flashFluidState.pressure(0)
490 + (pC[phaseIdx] - pC[0]));
491
492 // update the parameter cache
493 paramCache.updateAll(flashFluidState, /*except=*/ParamCache::Temperature);
494
495 // update all densities and fugacity coefficients
496 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
497 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
498 flashFluidState.setDensity(phaseIdx, rho);
499
500 for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
501 const FlashEval& phi = FluidSystem::fugacityCoefficient(flashFluidState, paramCache, phaseIdx, compIdx);
502 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
503 }
504 }
505 }
506
507 static bool isPressureIdx_(unsigned pvIdx)
508 { return pvIdx == 0; }
509
510 static bool isSaturationIdx_(unsigned pvIdx)
511 { return 1 <= pvIdx && pvIdx < numPhases; }
512
513 static bool isMoleFracIdx_(unsigned pvIdx)
514 { return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
515
516 // retrieves a quantity from the fluid state
517 template <class FluidState>
518 static const typename FluidState::Scalar& getQuantity_(const FluidState& fluidState, unsigned pvIdx)
519 {
520 assert(pvIdx < numEq);
521
522 // first pressure
523 if (pvIdx < 1) {
524 unsigned phaseIdx = 0;
525 return fluidState.pressure(phaseIdx);
526 }
527 // first M - 1 saturations
528 else if (pvIdx < numPhases) {
529 unsigned phaseIdx = pvIdx - 1;
530 return fluidState.saturation(phaseIdx);
531 }
532 // mole fractions
533 else // if (pvIdx < numPhases + numPhases*numComponents)
534 {
535 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
536 unsigned compIdx = (pvIdx - numPhases)%numComponents;
537 return fluidState.moleFraction(phaseIdx, compIdx);
538 }
539 }
540
541 // set a quantity in the fluid state
542 template <class FluidState>
543 static void setQuantity_(FluidState& fluidState,
544 unsigned pvIdx,
545 const typename FluidState::Scalar& value)
546 {
547 assert(pvIdx < numEq);
548
549 Valgrind::CheckDefined(value);
550 // first pressure
551 if (pvIdx < 1) {
552 unsigned phaseIdx = 0;
553 fluidState.setPressure(phaseIdx, value);
554 }
555 // first M - 1 saturations
556 else if (pvIdx < numPhases) {
557 unsigned phaseIdx = pvIdx - 1;
558 fluidState.setSaturation(phaseIdx, value);
559 }
560 // mole fractions
561 else {
562 assert(pvIdx < numPhases + numPhases*numComponents);
563 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
564 unsigned compIdx = (pvIdx - numPhases)%numComponents;
565 fluidState.setMoleFraction(phaseIdx, compIdx, value);
566 }
567 }
568
569 template <class FluidState>
570 static Scalar quantityWeight_(const FluidState& /*fluidState*/, unsigned pvIdx)
571 {
572 // first pressure
573 if (pvIdx < 1)
574 return 1e-6;
575 // first M - 1 saturations
576 else if (pvIdx < numPhases)
577 return 1.0;
578 // mole fractions
579 else // if (pvIdx < numPhases + numPhases*numComponents)
580 return 1.0;
581 }
582};
583
584} // namespace Opm
585
586#endif
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Representation of an evaluation of a function and its derivatives w.r.t.
Provides the OPM specific exception classes.
This file contains helper classes for the material laws.
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...
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Some templates to wrap the valgrind client request macros.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition CompositionalFluidState.hpp:46
Represents a function evaluation and its derivatives w.r.t.
Definition Evaluation.hpp:57
Determines the phase compositions, pressures and saturations given the total mass of all components.
Definition NcpFlash.hpp:88
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &globalMolarities)
Guess initial values for all quantities.
Definition NcpFlash.hpp:105
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.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition NcpFlash.hpp:146
static void solve(FluidState &fluidState, const ComponentVector &globalMolarities, Scalar tolerance=0.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition NcpFlash.hpp:257
A generic traits class which does not provide any indices.
Definition MaterialTraits.hpp:45
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Definition NullMaterial.hpp:46
Definition Exceptions.hpp:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30