My Project
Loading...
Searching...
No Matches
PTFlash.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 Copyright 2022 NORCE.
5 Copyright 2022 SINTEF Digital, Mathematics and Cybernetics.
6
7 This file is part of the Open Porous Media project (OPM).
8
9 OPM is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 2 of the License, or
12 (at your option) any later version.
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20
21 Consult the COPYING file in the top-level source directory of this
22 module for the precise wording of the license and the list of
23 copyright holders.
24*/
29#ifndef OPM_CHI_FLASH_HPP
30#define OPM_CHI_FLASH_HPP
31
41
42#include <dune/common/fvector.hh>
43#include <dune/common/fmatrix.hh>
44#include <dune/common/classname.hh>
45
46#include <limits>
47#include <iostream>
48#include <iomanip>
49#include <stdexcept>
50#include <type_traits>
51
52#include <fmt/format.h>
53
54namespace Opm {
55
61template <class Scalar, class FluidSystem>
63{
64 enum { numPhases = FluidSystem::numPhases };
65 enum { numComponents = FluidSystem::numComponents };
66 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx};
67 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx};
68 enum { numMiscibleComponents = FluidSystem::numMiscibleComponents};
69 enum { numMisciblePhases = FluidSystem::numMisciblePhases}; //oil, gas
70 enum {
71 numEq =
72 numMisciblePhases+
73 numMisciblePhases*numMiscibleComponents
74 };
75
76public:
81 template <class FluidState>
82 static void solve(FluidState& fluid_state,
83 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& z,
84 std::string twoPhaseMethod,
85 Scalar tolerance = -1.,
86 int verbosity = 0)
87 {
88
89 using InputEval = typename FluidState::Scalar;
90 using ComponentVector = Dune::FieldVector<typename FluidState::Scalar, numComponents>;
91
92 if (tolerance <= 0) {
93 tolerance = std::min<Scalar>(1e-3, 1e8 * std::numeric_limits<Scalar>::epsilon());
94 }
95
96 // K and L from previous timestep (wilson and -1 initially)
97 ComponentVector K;
98 for(int compIdx = 0; compIdx < numComponents; ++compIdx) {
99 K[compIdx] = fluid_state.K(compIdx);
100 }
101 InputEval L;
102 // TODO: L has all the derivatives to be all ZEROs here.
103 L = fluid_state.L();
104
105 // Print header
106 if (verbosity >= 1) {
107 std::cout << "********" << std::endl;
108 std::cout << "Inputs are K = [" << K << "], L = [" << L << "], z = [" << z << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl;
109 }
110
111 using ScalarVector = Dune::FieldVector<Scalar, numComponents>;
112 Scalar L_scalar = Opm::getValue(L);
113 ScalarVector z_scalar;
114 ScalarVector K_scalar;
115 for (unsigned i = 0; i < numComponents; ++i) {
116 z_scalar[i] = Opm::getValue(z[i]);
117 K_scalar[i] = Opm::getValue(K[i]);
118 }
119 using ScalarFluidState = CompositionalFluidState<Scalar, FluidSystem>;
120 ScalarFluidState fluid_state_scalar;
121
122 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
123 fluid_state_scalar.setMoleFraction(oilPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(oilPhaseIdx, compIdx)));
124 fluid_state_scalar.setMoleFraction(gasPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(gasPhaseIdx, compIdx)));
125 fluid_state_scalar.setKvalue(compIdx, Opm::getValue(fluid_state.K(compIdx)));
126 }
127
128 fluid_state_scalar.setLvalue(L_scalar);
129 // other values need to be Scalar, but I guess the fluidstate does not support it yet.
130 fluid_state_scalar.setPressure(FluidSystem::oilPhaseIdx,
131 Opm::getValue(fluid_state.pressure(FluidSystem::oilPhaseIdx)));
132 fluid_state_scalar.setPressure(FluidSystem::gasPhaseIdx,
133 Opm::getValue(fluid_state.pressure(FluidSystem::gasPhaseIdx)));
134
135 fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0)));
136
137 // Do a stability test to check if cell is is_single_phase-phase (do for all cells the first time).
138 bool is_stable = false;
139 if ( L <= 0 || L == 1 ) {
140 if (verbosity >= 1) {
141 std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl;
142 }
143 phaseStabilityTest_(is_stable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
144 }
145 if (verbosity >= 1) {
146 std::cout << "Inputs after stability test are K = [" << K_scalar << "], L = [" << L_scalar << "], z = [" << z_scalar << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl;
147 }
148 // TODO: we do not need two variables is_stable and is_single_hase, while lacking a good name
149 // TODO: from the later code, is good if we knows whether single_phase_gas or single_phase_oil here
150 const bool is_single_phase = is_stable;
151
152 // Update the composition if cell is two-phase
153 if ( !is_single_phase ) {
154 // Rachford Rice equation to get initial L for composition solver
155 L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity);
156 flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity);
157 } else {
158 // Cell is one-phase. Use Li's phase labeling method to see if it's liquid or vapor
159 L_scalar = li_single_phase_label_(fluid_state_scalar, z_scalar, verbosity);
160 }
161 fluid_state_scalar.setLvalue(L_scalar);
162
163 // Print footer
164 if (verbosity >= 1) {
165 std::cout << "********" << std::endl;
166 }
167
168 // the flash solution process were performed in scalar form, after the flash calculation finishes,
169 // ensure that things in fluid_state_scalar is transformed to fluid_state
170 for (int compIdx=0; compIdx<numComponents; ++compIdx){
171 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, compIdx);
172 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x_i);
173 const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, compIdx);
174 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y_i);
175 }
176
177 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
178 fluid_state.setKvalue(compIdx, K_scalar[compIdx]);
179 fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]);
180 }
181 fluid_state.setLvalue(L_scalar);
182 // we update the derivatives in fluid_state
183 updateDerivatives_(fluid_state_scalar, z, fluid_state, is_single_phase);
184 }//end solve
185
193 template <class FluidState, class ComponentVector>
194 static void solve(FluidState& fluid_state,
195 const ComponentVector& globalMolarities,
196 Scalar tolerance = 0.0)
197 {
198 using MaterialTraits = NullMaterialTraits<Scalar, numPhases>;
199 using MaterialLaw = NullMaterial<MaterialTraits>;
200 using MaterialLawParams = typename MaterialLaw::Params;
201
202 MaterialLawParams matParams;
203 solve<MaterialLaw>(fluid_state, matParams, globalMolarities, tolerance);
204 }
205
206 template <class Vector>
207 static typename Vector::field_type solveRachfordRice_g_(const Vector& K, const Vector& z, int verbosity)
208 {
209 // Find min and max K. Have to do a laborious for loop to avoid water component (where K=0)
210 // TODO: Replace loop with Dune::min_value() and Dune::max_value() when water component is properly handled
211 using field_type = typename Vector::field_type;
212 constexpr field_type tol = 1e-12;
213 constexpr int itmax = 10000;
214 field_type Kmin = K[0];
215 field_type Kmax = K[0];
216 for (int compIdx = 1; compIdx < numComponents; ++compIdx){
217 if (K[compIdx] < Kmin)
218 Kmin = K[compIdx];
219 else if (K[compIdx] >= Kmax)
220 Kmax = K[compIdx];
221 }
222 // Lower and upper bound for solution
223 auto Vmin = 1 / (1 - Kmax);
224 auto Vmax = 1 / (1 - Kmin);
225 // Initial guess
226 auto V = (Vmin + Vmax)/2;
227 // Print initial guess and header
228 if (verbosity == 3 || verbosity == 4) {
229 std::cout << "Initial guess " << numComponents << "c : V = " << V << " and [Vmin, Vmax] = [" << Vmin << ", " << Vmax << "]" << std::endl;
230 std::cout << std::setw(10) << "Iteration" << std::setw(16) << "abs(step)" << std::setw(16) << "V" << std::endl;
231 }
232 // Newton-Raphson loop
233 for (int iteration = 1; iteration < itmax; ++iteration) {
234 // Calculate function and derivative values
235 field_type denum = 0.0;
236 field_type r = 0.0;
237 for (int compIdx = 0; compIdx < numComponents; ++compIdx){
238 auto dK = K[compIdx] - 1.0;
239 auto a = z[compIdx] * dK;
240 auto b = (1 + V * dK);
241 r += a/b;
242 denum += z[compIdx] * (dK*dK) / (b*b);
243 }
244 auto delta = r / denum;
245 V += delta;
246
247 // Check if V is within the bounds, and if not, we apply bisection method
248 if (V < Vmin || V > Vmax)
249 {
250 // Print info
251 if (verbosity == 3 || verbosity == 4) {
252 std::cout << "V = " << V << " is not within the the range [Vmin, Vmax], solve using Bisection method!" << std::endl;
253 }
254
255 // Run bisection
256 auto Lmin = 1 - Vmax;
257 auto Lmax = 1 - Vmin;
258 // TODO: This is required for some cases. Not clear why
259 // since the objective function should be monotone with a
260 // single zero between the Lmin/Lmax interval defined by
261 // K-values.
262 Lmax = 1.0;
263 Lmin = 0.0;
264 auto L = bisection_g_(K, Lmin, Lmax, z, verbosity);
265
266 // Print final result
267 if (verbosity >= 1) {
268 std::cout << "Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl;
269 }
270 return L;
271 }
272
273 // Print iteration info
274 if (verbosity == 3 || verbosity == 4) {
275 std::cout << std::setw(10) << iteration << std::setw(16) << Opm::abs(delta) << std::setw(16) << V << std::endl;
276 }
277 // Check for convergence
278 if ( Opm::abs(r) < tol ) {
279 auto L = 1 - V;
280 // Should we make sure the range of L is within (0, 1)?
281
282 // Print final result
283 if (verbosity >= 1) {
284 std::cout << "Rachford-Rice converged to final solution L = " << L << std::endl;
285 }
286 return L;
287 }
288 }
289
290 // Throw error if Rachford-Rice fails
291 throw std::runtime_error(" Rachford-Rice did not converge within maximum number of iterations" );
292 }
293
294 template <class Vector>
295 static typename Vector::field_type bisection_g_(const Vector& K, typename Vector::field_type Lmin,
296 typename Vector::field_type Lmax, const Vector& z, int verbosity)
297 {
298 // Calculate for g(Lmin) for first comparison with gMid = g(L)
299 typename Vector::field_type gLmin = rachfordRice_g_(K, Lmin, z);
300
301 // Print new header
302 if (verbosity >= 3) {
303 std::cout << std::setw(10) << "Iteration" << std::setw(16) << "g(Lmid)" << std::setw(16) << "L" << std::endl;
304 }
305
306 constexpr int max_it = 10000;
307
308 auto closeLmaxLmin = [](double max_v, double min_v) {
309 return Opm::abs(max_v - min_v) / 2. < 1e-10;
310 // what if max_v < min_v?
311 };
312
313 // Bisection loop
314 if (closeLmaxLmin(Lmax, Lmin) ){
315 throw std::runtime_error(fmt::format("Strange bisection with Lmax {} and Lmin {}?", Lmax, Lmin));
316 }
317 for (int iteration = 0; iteration < max_it; ++iteration){
318 // New midpoint
319 auto L = (Lmin + Lmax) / 2;
320 auto gMid = rachfordRice_g_(K, L, z);
321 // std::cout << ">>> Lmin = " << Lmin << "g(Lmin) = " << gLmin << " L = " << L << " g(L) = " << gMid << std::endl;
322 if (verbosity == 3 || verbosity == 4) {
323 std::cout << std::setw(10) << iteration << std::setw(16) << gMid << std::setw(16) << L << std::endl;
324 }
325
326 // Check if midpoint fulfills g=0 or L - Lmin is sufficiently small
327 if (Opm::abs(gMid) < 1e-16 || closeLmaxLmin(Lmax, Lmin)){
328 return L;
329 }
330 // Else we repeat with midpoint being either Lmin og Lmax (depending on the signs).
331 else if (Dune::sign(gMid) != Dune::sign(gLmin)) {
332 // gMid has different sign as gLmin, so we set L as the new Lmax
333 Lmax = L;
334 }
335 else {
336 // gMid and gLmin have same sign so we set L as the new Lmin
337 Lmin = L;
338 gLmin = gMid;
339 }
340 }
341 throw std::runtime_error(" Rachford-Rice bisection failed with " + std::to_string(max_it) + " iterations!");
342 }
343
344 template <class Vector, class FlashFluidState>
345 static typename Vector::field_type li_single_phase_label_(const FlashFluidState& fluid_state, const Vector& z, int verbosity)
346 {
347 // Calculate intermediate sum
348 typename Vector::field_type sumVz = 0.0;
349 for (int compIdx=0; compIdx<numComponents; ++compIdx){
350 // Get component information
351 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
352
353 // Sum calculation
354 sumVz += (V_crit * z[compIdx]);
355 }
356
357 // Calculate approximate (pseudo) critical temperature using Li's method
358 typename Vector::field_type Tc_est = 0.0;
359 for (int compIdx=0; compIdx<numComponents; ++compIdx){
360 // Get component information
361 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
362 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
363
364 // Sum calculation
365 Tc_est += (V_crit * T_crit * z[compIdx] / sumVz);
366 }
367
368 // Get temperature
369 const auto& T = fluid_state.temperature(0);
370
371 // If temperature is below estimated critical temperature --> phase = liquid; else vapor
372 typename Vector::field_type L;
373 if (T < Tc_est) {
374 // Liquid
375 L = 1.0;
376
377 // Print
378 if (verbosity >= 1) {
379 std::cout << "Cell is single-phase, liquid (L = 1.0) due to Li's phase labeling method giving T < Tc_est (" << T << " < " << Tc_est << ")!" << std::endl;
380 }
381 }
382 else {
383 // Vapor
384 L = 0.0;
385
386 // Print
387 if (verbosity >= 1) {
388 std::cout << "Cell is single-phase, vapor (L = 0.0) due to Li's phase labeling method giving T >= Tc_est (" << T << " >= " << Tc_est << ")!" << std::endl;
389 }
390 }
391
392
393 return L;
394 }
395
396 template <class FlashFluidState, class ComponentVector>
397 static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluid_state, const ComponentVector& z, int verbosity)
398 {
399 // Declarations
400 bool isTrivialL, isTrivialV;
401 ComponentVector x, y;
402 typename FlashFluidState::Scalar S_l, S_v;
403 ComponentVector K0 = K;
404 ComponentVector K1 = K;
405
406 // Check for vapour instable phase
407 if (verbosity == 3 || verbosity == 4) {
408 std::cout << "Stability test for vapor phase:" << std::endl;
409 }
410 checkStability_(fluid_state, isTrivialV, K0, y, S_v, z, /*isGas=*/true, verbosity);
411 bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV;
412
413 // Check for liquids stable phase
414 if (verbosity == 3 || verbosity == 4) {
415 std::cout << "Stability test for liquid phase:" << std::endl;
416 }
417 checkStability_(fluid_state, isTrivialL, K1, x, S_l, z, /*isGas=*/false, verbosity);
418 bool L_stable = (S_l < (1.0 + 1e-5)) || isTrivialL;
419
420 // L-stable means success in making liquid, V-unstable means no success in making vapour
421 isStable = L_stable && V_unstable;
422 if (isStable) {
423 // Single phase, i.e. phase composition is equivalent to the global composition
424 // Update fluid_state with mole fraction
425 for (int compIdx=0; compIdx<numComponents; ++compIdx){
426 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, z[compIdx]);
427 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, z[compIdx]);
428 }
429 }
430 // If not stable: use the mole fractions from Michelsen's test to update K
431 else {
432 for (int compIdx = 0; compIdx<numComponents; ++compIdx) {
433 K[compIdx] = y[compIdx] / x[compIdx];
434 }
435 }
436 }
437
438protected:
439
440 template <class FlashFluidState>
441 static typename FlashFluidState::Scalar wilsonK_(const FlashFluidState& fluid_state, int compIdx)
442 {
443 const auto& acf = FluidSystem::acentricFactor(compIdx);
444 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
445 const auto& T = fluid_state.temperature(0);
446 const auto& p_crit = FluidSystem::criticalPressure(compIdx);
447 const auto& p = fluid_state.pressure(0); //for now assume no capillary pressure
448
449 const auto& tmp = Opm::exp(5.3727 * (1+acf) * (1-T_crit/T)) * (p_crit/p);
450 return tmp;
451 }
452
453 template <class Vector>
454 static typename Vector::field_type rachfordRice_g_(const Vector& K, typename Vector::field_type L, const Vector& z)
455 {
456 typename Vector::field_type g=0;
457 for (int compIdx=0; compIdx<numComponents; ++compIdx){
458 g += (z[compIdx]*(K[compIdx]-1))/(K[compIdx]-L*(K[compIdx]-1));
459 }
460 return g;
461 }
462
463 template <class Vector>
464 static typename Vector::field_type rachfordRice_dg_dL_(const Vector& K, const typename Vector::field_type L, const Vector& z)
465 {
466 typename Vector::field_type dg=0;
467 for (int compIdx=0; compIdx<numComponents; ++compIdx){
468 dg += (z[compIdx]*(K[compIdx]-1)*(K[compIdx]-1))/((K[compIdx]-L*(K[compIdx]-1))*(K[compIdx]-L*(K[compIdx]-1)));
469 }
470 return dg;
471 }
472
473 template <class FlashFluidState, class ComponentVector>
474 static void checkStability_(const FlashFluidState& fluid_state, bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc,
475 typename FlashFluidState::Scalar& S_loc, const ComponentVector& z, bool isGas, int verbosity)
476 {
477 using FlashEval = typename FlashFluidState::Scalar;
478 using PengRobinsonMixture = typename Opm::PengRobinsonMixture<Scalar, FluidSystem>;
479
480 // Declarations
481 FlashFluidState fluid_state_fake = fluid_state;
482 FlashFluidState fluid_state_global = fluid_state;
483
484 // Setup output
485 if (verbosity >= 3 || verbosity >= 4) {
486 std::cout << std::setw(10) << "Iteration" << std::setw(16) << "K-Norm" << std::setw(16) << "R-Norm" << std::endl;
487 }
488
489 // Michelsens stability test.
490 // Make two fake phases "inside" one phase and check for positive volume
491 for (int i = 0; i < 20000; ++i) {
492 S_loc = 0.0;
493 if (isGas) {
494 for (int compIdx=0; compIdx<numComponents; ++compIdx){
495 xy_loc[compIdx] = K[compIdx] * z[compIdx];
496 S_loc += xy_loc[compIdx];
497 }
498 for (int compIdx=0; compIdx<numComponents; ++compIdx){
499 xy_loc[compIdx] /= S_loc;
500 fluid_state_fake.setMoleFraction(gasPhaseIdx, compIdx, xy_loc[compIdx]);
501 }
502 }
503 else {
504 for (int compIdx=0; compIdx<numComponents; ++compIdx){
505 xy_loc[compIdx] = z[compIdx]/K[compIdx];
506 S_loc += xy_loc[compIdx];
507 }
508 for (int compIdx=0; compIdx<numComponents; ++compIdx){
509 xy_loc[compIdx] /= S_loc;
510 fluid_state_fake.setMoleFraction(oilPhaseIdx, compIdx, xy_loc[compIdx]);
511 }
512 }
513
514 int phaseIdx = (isGas ? static_cast<int>(gasPhaseIdx) : static_cast<int>(oilPhaseIdx));
515 int phaseIdx2 = (isGas ? static_cast<int>(oilPhaseIdx) : static_cast<int>(gasPhaseIdx));
516 // TODO: not sure the following makes sense
517 for (int compIdx=0; compIdx<numComponents; ++compIdx){
518 fluid_state_global.setMoleFraction(phaseIdx2, compIdx, z[compIdx]);
519 }
520
521 typename FluidSystem::template ParameterCache<FlashEval> paramCache_fake;
522 paramCache_fake.updatePhase(fluid_state_fake, phaseIdx);
523
524 typename FluidSystem::template ParameterCache<FlashEval> paramCache_global;
525 paramCache_global.updatePhase(fluid_state_global, phaseIdx2);
526
527 //fugacity for fake phases each component
528 for (int compIdx=0; compIdx<numComponents; ++compIdx){
529 auto phiFake = PengRobinsonMixture::computeFugacityCoefficient(fluid_state_fake, paramCache_fake, phaseIdx, compIdx);
530 auto phiGlobal = PengRobinsonMixture::computeFugacityCoefficient(fluid_state_global, paramCache_global, phaseIdx2, compIdx);
531
532 fluid_state_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake);
533 fluid_state_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
534 }
535
536
537 ComponentVector R;
538 for (int compIdx=0; compIdx<numComponents; ++compIdx){
539 if (isGas){
540 auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
541 auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
542 auto fug_ratio = fug_global / fug_fake;
543 R[compIdx] = fug_ratio / S_loc;
544 }
545 else{
546 auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
547 auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
548 auto fug_ratio = fug_fake / fug_global;
549 R[compIdx] = fug_ratio * S_loc;
550 }
551 }
552
553 for (int compIdx=0; compIdx<numComponents; ++compIdx){
554 K[compIdx] *= R[compIdx];
555 }
556 Scalar R_norm = 0.0;
557 Scalar K_norm = 0.0;
558 for (int compIdx=0; compIdx<numComponents; ++compIdx){
559 auto a = Opm::getValue(R[compIdx]) - 1.0;
560 auto b = Opm::log(Opm::getValue(K[compIdx]));
561 R_norm += a*a;
562 K_norm += b*b;
563 }
564
565 // Print iteration info
566 if (verbosity >= 3) {
567 std::cout << std::setw(10) << i << std::setw(16) << K_norm << std::setw(16) << R_norm << std::endl;
568 }
569
570 // Check convergence
571 isTrivial = (K_norm < 1e-5);
572 if (isTrivial || R_norm < 1e-10)
573 return;
574 //todo: make sure that no mole fraction is smaller than 1e-8 ?
575 //todo: take care of water!
576 }
577 throw std::runtime_error(" Stability test did not converge");
578 }//end checkStability
579
580 // TODO: basically FlashFluidState and ComponentVector are both depending on the one Scalar type
581 template <class FlashFluidState, class ComponentVector>
582 static void computeLiquidVapor_(FlashFluidState& fluid_state, typename FlashFluidState::Scalar& L, ComponentVector& K, const ComponentVector& z)
583 {
584 // Calculate x and y, and normalize
585 ComponentVector x;
586 ComponentVector y;
587 typename FlashFluidState::Scalar sumx=0;
588 typename FlashFluidState::Scalar sumy=0;
589 for (int compIdx=0; compIdx<numComponents; ++compIdx){
590 x[compIdx] = z[compIdx]/(L + (1-L)*K[compIdx]);
591 sumx += x[compIdx];
592 y[compIdx] = (K[compIdx]*z[compIdx])/(L + (1-L)*K[compIdx]);
593 sumy += y[compIdx];
594 }
595 x /= sumx;
596 y /= sumy;
597
598 for (int compIdx=0; compIdx<numComponents; ++compIdx){
599 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]);
600 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y[compIdx]);
601 }
602 }
603
604 template <class FluidState, class ComponentVector>
605 static void flash_2ph(const ComponentVector& z_scalar,
606 const std::string& flash_2p_method,
607 ComponentVector& K_scalar,
608 typename FluidState::Scalar& L_scalar,
609 FluidState& fluid_state_scalar,
610 int verbosity = 0) {
611 if (verbosity >= 1) {
612 std::cout << "Cell is two-phase! Solve Rachford-Rice with initial K = [" << K_scalar << "]" << std::endl;
613 }
614
615 // Calculate composition using nonlinear solver
616 // Newton
617 if (flash_2p_method == "newton") {
618 if (verbosity >= 1) {
619 std::cout << "Calculate composition using Newton." << std::endl;
620 }
621 newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
622 } else if (flash_2p_method == "ssi") {
623 // Successive substitution
624 if (verbosity >= 1) {
625 std::cout << "Calculate composition using Succcessive Substitution." << std::endl;
626 }
627 successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, false, verbosity);
628 } else if (flash_2p_method == "ssi+newton") {
629 successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, true, verbosity);
630 newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
631 } else {
632 throw std::runtime_error("unknown two phase flash method " + flash_2p_method + " is specified");
633 }
634 }
635
636 template <class FlashFluidState, class ComponentVector>
637 static void newtonComposition_(ComponentVector& K, typename FlashFluidState::Scalar& L,
638 FlashFluidState& fluid_state, const ComponentVector& z,
639 int verbosity)
640 {
641 // Note: due to the need for inverse flash update for derivatives, the following two can be different
642 // Looking for a good way to organize them
643 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
644 constexpr size_t num_primary_variables = numMisciblePhases * numMiscibleComponents + 1;
645 using NewtonVector = Dune::FieldVector<Scalar, num_equations>;
646 using NewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, num_primary_variables>;
647 constexpr Scalar tolerance = 1.e-8;
648
649 NewtonVector soln(0.);
650 NewtonVector res(0.);
651 NewtonMatrix jac (0.);
652
653 // Compute x and y from K, L and Z
654 computeLiquidVapor_(fluid_state, L, K, z);
655 if (verbosity >= 1) {
656 std::cout << " the current L is " << Opm::getValue(L) << std::endl;
657 }
658
659 // Print initial condition
660 if (verbosity >= 1) {
661 std::cout << "Initial guess: x = [";
662 for (int compIdx=0; compIdx<numComponents; ++compIdx){
663 if (compIdx < numComponents - 1)
664 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) << " ";
665 else
666 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
667 }
668 std::cout << "], y = [";
669 for (int compIdx=0; compIdx<numComponents; ++compIdx){
670 if (compIdx < numComponents - 1)
671 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) << " ";
672 else
673 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
674 }
675 std::cout << "], and " << "L = " << L << std::endl;
676 }
677
678 // Print header
679 if (verbosity == 2 || verbosity == 4) {
680 std::cout << std::setw(10) << "Iteration" << std::setw(16) << "Norm2(step)" << std::setw(16) << "Norm2(Residual)" << std::endl;
681 }
682
683 // AD type
684 using Eval = DenseAd::Evaluation<Scalar, num_primary_variables>;
685 // TODO: we might need to use numMiscibleComponents here
686 std::vector<Eval> x(numComponents), y(numComponents);
687 Eval l;
688
689 // TODO: I might not need to set soln anything here.
690 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
691 x[compIdx] = Eval(fluid_state.moleFraction(oilPhaseIdx, compIdx), compIdx);
692 const unsigned idx = compIdx + numComponents;
693 y[compIdx] = Eval(fluid_state.moleFraction(gasPhaseIdx, compIdx), idx);
694 }
695 l = Eval(L, num_primary_variables - 1);
696
697 // it is created for the AD calculation for the flash calculation
698 CompositionalFluidState<Eval, FluidSystem> flash_fluid_state;
699 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
700 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
701 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
702 // TODO: should we use wilsonK_ here?
703 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
704 }
705 flash_fluid_state.setLvalue(l);
706 // other values need to be Scalar, but I guess the fluid_state does not support it yet.
707 flash_fluid_state.setPressure(FluidSystem::oilPhaseIdx,
708 fluid_state.pressure(FluidSystem::oilPhaseIdx));
709 flash_fluid_state.setPressure(FluidSystem::gasPhaseIdx,
710 fluid_state.pressure(FluidSystem::gasPhaseIdx));
711
712 // TODO: not sure whether we need to set the saturations
713 flash_fluid_state.setSaturation(FluidSystem::gasPhaseIdx,
714 fluid_state.saturation(FluidSystem::gasPhaseIdx));
715 flash_fluid_state.setSaturation(FluidSystem::oilPhaseIdx,
716 fluid_state.saturation(FluidSystem::oilPhaseIdx));
717 flash_fluid_state.setTemperature(fluid_state.temperature(0));
718
719 using ParamCache = typename FluidSystem::template ParameterCache<typename CompositionalFluidState<Eval, FluidSystem>::Scalar>;
720 ParamCache paramCache;
721
722 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
723 paramCache.updatePhase(flash_fluid_state, phaseIdx);
724 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
725 // TODO: will phi here carry the correct derivatives?
726 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
727 flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
728 }
729 }
730 bool converged = false;
731 unsigned iter = 0;
732 constexpr unsigned max_iter = 1000;
733 while (iter < max_iter) {
734 assembleNewton_<CompositionalFluidState<Eval, FluidSystem>, ComponentVector, num_primary_variables, num_equations>
735 (flash_fluid_state, z, jac, res);
736 if (verbosity >= 1) {
737 std::cout << " newton residual is " << res.two_norm() << std::endl;
738 }
739 converged = res.two_norm() < tolerance;
740 if (converged) {
741 break;
742 }
743
744 jac.solve(soln, res);
745 constexpr Scalar damping_factor = 1.0;
746 // updating x and y
747 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
748 x[compIdx] -= soln[compIdx] * damping_factor;
749 y[compIdx] -= soln[compIdx + numComponents] * damping_factor;
750 }
751 l -= soln[num_equations - 1] * damping_factor;
752
753 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
754 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
755 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
756 // TODO: should we use wilsonK_ here?
757 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
758 }
759 flash_fluid_state.setLvalue(l);
760
761 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
762 paramCache.updatePhase(flash_fluid_state, phaseIdx);
763 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
764 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
765 flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
766 }
767 }
768 }
769 if (verbosity >= 1) {
770 for (unsigned i = 0; i < num_equations; ++i) {
771 for (unsigned j = 0; j < num_primary_variables; ++j) {
772 std::cout << " " << jac[i][j];
773 }
774 std::cout << std::endl;
775 }
776 std::cout << std::endl;
777 }
778 if (!converged) {
779 throw std::runtime_error(" Newton composition update did not converge within maxIterations " + std::to_string(max_iter));
780 }
781
782 // fluid_state is scalar, we need to update all the values for fluid_state here
783 for (unsigned idx = 0; idx < numComponents; ++idx) {
784 const auto x_i = Opm::getValue(flash_fluid_state.moleFraction(oilPhaseIdx, idx));
785 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
786 const auto y_i = Opm::getValue(flash_fluid_state.moleFraction(gasPhaseIdx, idx));
787 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
788 const auto K_i = Opm::getValue(flash_fluid_state.K(idx));
789 fluid_state.setKvalue(idx, K_i);
790 // TODO: not sure we need K and L here, because they are in the flash_fluid_state anyway.
791 K[idx] = K_i;
792 }
793 L = Opm::getValue(l);
794 fluid_state.setLvalue(L);
795 }
796
797 // TODO: the interface will need to refactor for later usage
798 template<typename FlashFluidState, typename ComponentVector, size_t num_primary, size_t num_equation >
799 static void assembleNewton_(const FlashFluidState& fluid_state,
800 const ComponentVector& global_composition,
801 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
802 Dune::FieldVector<double, num_equation>& res)
803 {
804 using Eval = DenseAd::Evaluation<double, num_primary>;
805 std::vector<Eval> x(numComponents), y(numComponents);
806 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
807 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
808 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
809 }
810 const Eval& l = fluid_state.L();
811 // TODO: clearing zero whether necessary?
812 jac = 0.;
813 res = 0.;
814 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
815 {
816 // z - L*x - (1-L) * y
817 auto local_res = -global_composition[compIdx] + l * x[compIdx] + (1 - l) * y[compIdx];
818 res[compIdx] = Opm::getValue(local_res);
819 for (unsigned i = 0; i < num_primary; ++i) {
820 jac[compIdx][i] = local_res.derivative(i);
821 }
822 }
823
824 {
825 // f_liquid - f_vapor = 0
826 auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
827 fluid_state.fugacity(gasPhaseIdx, compIdx));
828 res[compIdx + numComponents] = Opm::getValue(local_res);
829 for (unsigned i = 0; i < num_primary; ++i) {
830 jac[compIdx + numComponents][i] = local_res.derivative(i);
831 }
832 }
833 }
834 // sum(x) - sum(y) = 0
835 Eval sumx = 0.;
836 Eval sumy = 0.;
837 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
838 sumx += x[compIdx];
839 sumy += y[compIdx];
840 }
841 auto local_res = sumx - sumy;
842 res[num_equation - 1] = Opm::getValue(local_res);
843 for (unsigned i = 0; i < num_primary; ++i) {
844 jac[num_equation - 1][i] = local_res.derivative(i);
845 }
846 }
847
848 // TODO: the interface will need to refactor for later usage
849 template<typename FlashFluidState, typename ComponentVector, size_t num_primary, size_t num_equation >
850 static void assembleNewtonSingle_(const FlashFluidState& fluid_state,
851 const ComponentVector& global_composition,
852 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
853 Dune::FieldVector<double, num_equation>& res)
854 {
855 using Eval = DenseAd::Evaluation<double, num_primary>;
856 std::vector<Eval> x(numComponents), y(numComponents);
857 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
858 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
859 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
860 }
861 const Eval& l = fluid_state.L();
862
863 // TODO: clearing zero whether necessary?
864 jac = 0.;
865 res = 0.;
866 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
867 {
868 // z - L*x - (1-L) * y ---> z - x;
869 auto local_res = -global_composition[compIdx] + x[compIdx];
870 res[compIdx] = Opm::getValue(local_res);
871 for (unsigned i = 0; i < num_primary; ++i) {
872 jac[compIdx][i] = local_res.derivative(i);
873 }
874 }
875
876 {
877 // f_liquid - f_vapor = 0 -->z - y;
878 auto local_res = -global_composition[compIdx] + y[compIdx];
879 res[compIdx + numComponents] = Opm::getValue(local_res);
880 for (unsigned i = 0; i < num_primary; ++i) {
881 jac[compIdx + numComponents][i] = local_res.derivative(i);
882 }
883 }
884 }
885
886 // TODO: better we have isGas or isLiquid here
887 const bool isGas = Opm::abs(l - 1.0) > std::numeric_limits<double>::epsilon();
888
889 // sum(x) - sum(y) = 0
890 auto local_res = l;
891 if(isGas) {
892 local_res = l-1;
893 }
894
895 res[num_equation - 1] = Opm::getValue(local_res);
896 for (unsigned i = 0; i < num_primary; ++i) {
897 jac[num_equation - 1][i] = local_res.derivative(i);
898 }
899 }
900
901 template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
902 static void updateDerivatives_(const FlashFluidStateScalar& fluid_state_scalar,
903 const ComponentVector& z,
904 FluidState& fluid_state,
905 bool is_single_phase)
906 {
907 if(!is_single_phase)
908 updateDerivativesTwoPhase_(fluid_state_scalar, z, fluid_state);
909 else
910 updateDerivativesSinglePhase_(fluid_state_scalar, z, fluid_state);
911
912 }
913
914 template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
915 static void updateDerivativesTwoPhase_(const FlashFluidStateScalar& fluid_state_scalar,
916 const ComponentVector& z,
917 FluidState& fluid_state)
918 {
919 // getting the secondary Jocobian matrix
920 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
921 constexpr size_t secondary_num_pv = numComponents + 1; // pressure, z for all the components
922 using SecondaryEval = Opm::DenseAd::Evaluation<double, secondary_num_pv>; // three z and one pressure
923 using SecondaryComponentVector = Dune::FieldVector<SecondaryEval, numComponents>;
924 using SecondaryFlashFluidState = Opm::CompositionalFluidState<SecondaryEval, FluidSystem>;
925
926 SecondaryFlashFluidState secondary_fluid_state;
927 SecondaryComponentVector secondary_z;
928 // p and z are the primary variables here
929 // pressure
930 const SecondaryEval sec_p = SecondaryEval(fluid_state_scalar.pressure(FluidSystem::oilPhaseIdx), 0);
931 secondary_fluid_state.setPressure(FluidSystem::oilPhaseIdx, sec_p);
932 secondary_fluid_state.setPressure(FluidSystem::gasPhaseIdx, sec_p);
933
934 // set the temperature // TODO: currently we are not considering the temperature derivatives
935 secondary_fluid_state.setTemperature(Opm::getValue(fluid_state_scalar.temperature(0)));
936
937 for (unsigned idx = 0; idx < numComponents; ++idx) {
938 secondary_z[idx] = SecondaryEval(Opm::getValue(z[idx]), idx + 1);
939 }
940 // set up the mole fractions
941 for (unsigned idx = 0; idx < numComponents; ++idx) {
942 // TODO: double checking that fluid_state_scalar returns a scalar here
943 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, idx);
944 secondary_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
945 const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, idx);
946 secondary_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
947 // TODO: double checking make sure those are consistent
948 const auto K_i = fluid_state_scalar.K(idx);
949 secondary_fluid_state.setKvalue(idx, K_i);
950 }
951 const auto L = fluid_state_scalar.L();
952 secondary_fluid_state.setLvalue(L);
953 // TODO: Do we need to update the saturations?
954 // compositions
955 // TODO: we probably can simplify SecondaryFlashFluidState::Scalar
956 using SecondaryParamCache = typename FluidSystem::template ParameterCache<typename SecondaryFlashFluidState::Scalar>;
957 SecondaryParamCache secondary_param_cache;
958 for (unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
959 secondary_param_cache.updatePhase(secondary_fluid_state, phase_idx);
960 for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
961 SecondaryEval phi = FluidSystem::fugacityCoefficient(secondary_fluid_state, secondary_param_cache, phase_idx, comp_idx);
962 secondary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
963 }
964 }
965
966 using SecondaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
967 using SecondaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, secondary_num_pv>;
968 SecondaryNewtonMatrix sec_jac;
969 SecondaryNewtonVector sec_res;
970
971
972 //use the regular equations
973 assembleNewton_<SecondaryFlashFluidState, SecondaryComponentVector, secondary_num_pv, num_equations>
974 (secondary_fluid_state, secondary_z, sec_jac, sec_res);
975
976
977 // assembly the major matrix here
978 // primary variables are x, y and L
979 constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1;
981 using PrimaryComponentVector = Dune::FieldVector<double, numComponents>;
982 using PrimaryFlashFluidState = Opm::CompositionalFluidState<PrimaryEval, FluidSystem>;
983
984 PrimaryFlashFluidState primary_fluid_state;
985 // primary_z is not needed, because we use z will be okay here
986 PrimaryComponentVector primary_z;
987 for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
988 primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
989 }
990 for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
991 const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx);
992 primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii);
993 const unsigned idx = comp_idx + numComponents;
994 const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx);
995 primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii);
996 primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii);
997 }
998 PrimaryEval l;
999 l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1);
1000 primary_fluid_state.setLvalue(l);
1001 primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx));
1002 primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx));
1003 primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0));
1004
1005 // TODO: is PrimaryFlashFluidState::Scalar> PrimaryEval here?
1006 using PrimaryParamCache = typename FluidSystem::template ParameterCache<typename PrimaryFlashFluidState::Scalar>;
1007 PrimaryParamCache primary_param_cache;
1008 for (unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
1009 primary_param_cache.updatePhase(primary_fluid_state, phase_idx);
1010 for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
1011 PrimaryEval phi = FluidSystem::fugacityCoefficient(primary_fluid_state, primary_param_cache, phase_idx, comp_idx);
1012 primary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
1013 }
1014 }
1015
1016 using PrimaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
1017 using PrimaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, primary_num_pv>;
1018 PrimaryNewtonVector pri_res;
1019 PrimaryNewtonMatrix pri_jac;
1020
1021 //use the regular equations
1022 assembleNewton_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
1023 (primary_fluid_state, primary_z, pri_jac, pri_res);
1024
1025 // the following code does not compile with DUNE2.6
1026 // SecondaryNewtonMatrix xx;
1027 // pri_jac.solve(xx, sec_jac);
1028 pri_jac.invert();
1029 sec_jac.template leftmultiply(pri_jac);
1030
1031 using InputEval = typename FluidState::Scalar;
1032 using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>;
1033
1034 ComponentVectorMoleFraction x(numComponents), y(numComponents);
1035 InputEval L_eval = L;
1036
1037 // use the chainrule (and using partial instead of total
1038 // derivatives, DF / Dp = dF / dp + dF / ds * ds/dp.
1039 // where p is the primary variables and s the secondary variables. We then obtain
1040 // ds / dp = -inv(dF / ds)*(DF / Dp)
1041
1042 const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx);
1043 const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx);
1044 std::vector<double> K(numComponents);
1045
1046 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1047 K[compIdx] = fluid_state_scalar.K(compIdx);
1048 x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);//;z[compIdx] * 1. / (L + (1 - L) * K[compIdx]);
1049 y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);//;x[compIdx] * K[compIdx];
1050 }
1051
1052 // then we try to set the derivatives for x, y and L against P and x.
1053 // p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary
1054 // pressure joins
1055
1056 constexpr size_t num_deri = numComponents;
1057 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1058 std::vector<double> deri(num_deri, 0.);
1059 // derivatives from P
1060 for (unsigned idx = 0; idx < num_deri; ++idx) {
1061 deri[idx] = -sec_jac[compIdx][0] * p_l.derivative(idx);
1062 }
1063
1064 for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1065 const double pz = -sec_jac[compIdx][cIdx + 1];
1066 const auto& zi = z[cIdx];
1067 for (unsigned idx = 0; idx < num_deri; ++idx) {
1068 deri[idx] += pz * zi.derivative(idx);
1069 }
1070 }
1071 for (unsigned idx = 0; idx < num_deri; ++idx) {
1072 x[compIdx].setDerivative(idx, deri[idx]);
1073 }
1074 // handling y
1075 for (unsigned idx = 0; idx < num_deri; ++idx) {
1076 deri[idx] = -sec_jac[compIdx + numComponents][0] * p_v.derivative(idx);
1077 }
1078 for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1079 const double pz = -sec_jac[compIdx + numComponents][cIdx + 1];
1080 const auto& zi = z[cIdx];
1081 for (unsigned idx = 0; idx < num_deri; ++idx) {
1082 deri[idx] += pz * zi.derivative(idx);
1083 }
1084 }
1085 for (unsigned idx = 0; idx < num_deri; ++idx) {
1086 y[compIdx].setDerivative(idx, deri[idx]);
1087 }
1088
1089 // handling derivatives of L
1090 std::vector<double> deriL(num_deri, 0.);
1091 for (unsigned idx = 0; idx < num_deri; ++idx) {
1092 deriL[idx] = -sec_jac[2 * numComponents][0] * p_v.derivative(idx);
1093 }
1094 for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1095 const double pz = -sec_jac[2 * numComponents][cIdx + 1];
1096 const auto& zi = z[cIdx];
1097 for (unsigned idx = 0; idx < num_deri; ++idx) {
1098 deriL[idx] += pz * zi.derivative(idx);
1099 }
1100 }
1101
1102 for (unsigned idx = 0; idx < num_deri; ++idx) {
1103 L_eval.setDerivative(idx, deriL[idx]);
1104 }
1105 }
1106
1107 // set up the mole fractions
1108 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1109 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
1110 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
1111 }
1112 fluid_state.setLvalue(L_eval);
1113 } //end updateDerivativesTwoPhase
1114
1115 template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
1116 static void updateDerivativesSinglePhase_(const FlashFluidStateScalar& fluid_state_scalar,
1117 const ComponentVector& z,
1118 FluidState& fluid_state)
1119 {
1120 using InputEval = typename FluidState::Scalar;
1121 // L_eval is converted from a scalar, so all derivatives are zero at this point
1122 InputEval L_eval = fluid_state_scalar.L();;
1123
1124 // for single phase situation, x = y = z;
1125 // and L_eval have all zero derivatives
1126 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1127 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, z[compIdx]);
1128 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, z[compIdx]);
1129 }
1130 fluid_state.setLvalue(L_eval);
1131 } //end updateDerivativesSinglePhase
1132
1133
1134 // TODO: or use typename FlashFluidState::Scalar
1135 template <class FlashFluidState, class ComponentVector>
1136 static void successiveSubstitutionComposition_(ComponentVector& K, typename ComponentVector::field_type& L, FlashFluidState& fluid_state, const ComponentVector& z,
1137 const bool newton_afterwards, const int verbosity)
1138 {
1139 // Determine max. iterations based on if it will be used as a standalone flash or as a pre-process to Newton (or other) method.
1140 const int maxIterations = newton_afterwards ? 3 : 100;
1141
1142 // Store cout format before manipulation
1143 std::ios_base::fmtflags f(std::cout.flags());
1144
1145 // Print initial guess
1146 if (verbosity >= 1)
1147 std::cout << "Initial guess: K = [" << K << "] and L = " << L << std::endl;
1148
1149 if (verbosity == 2 || verbosity == 4) {
1150 // Print header
1151 int fugWidth = (numComponents * 12)/2;
1152 int convWidth = fugWidth + 7;
1153 std::cout << std::setw(10) << "Iteration" << std::setw(fugWidth) << "fL/fV" << std::setw(convWidth) << "norm2(fL/fv-1)" << std::endl;
1154 }
1155 //
1156 // Successive substitution loop
1157 //
1158 for (int i=0; i < maxIterations; ++i){
1159 // Compute (normalized) liquid and vapor mole fractions
1160 computeLiquidVapor_(fluid_state, L, K, z);
1161
1162 // Calculate fugacity coefficient
1163 using ParamCache = typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>;
1164 ParamCache paramCache;
1165 for (int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
1166 paramCache.updatePhase(fluid_state, phaseIdx);
1167 for (int compIdx=0; compIdx<numComponents; ++compIdx){
1168 auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, phaseIdx, compIdx);
1169 fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
1170 }
1171 }
1172
1173 // Calculate fugacity ratio
1174 ComponentVector newFugRatio;
1175 ComponentVector convFugRatio;
1176 for (int compIdx=0; compIdx<numComponents; ++compIdx){
1177 newFugRatio[compIdx] = fluid_state.fugacity(oilPhaseIdx, compIdx)/fluid_state.fugacity(gasPhaseIdx, compIdx);
1178 convFugRatio[compIdx] = newFugRatio[compIdx] - 1.0;
1179 }
1180
1181 // Print iteration info
1182 if (verbosity >= 2) {
1183 int prec = 5;
1184 int fugWidth = (prec + 3);
1185 int convWidth = prec + 9;
1186 std::cout << std::defaultfloat;
1187 std::cout << std::fixed;
1188 std::cout << std::setw(5) << i;
1189 std::cout << std::setw(fugWidth);
1190 std::cout << std::setprecision(prec);
1191 std::cout << newFugRatio;
1192 std::cout << std::scientific;
1193 std::cout << std::setw(convWidth) << convFugRatio.two_norm() << std::endl;
1194 }
1195
1196 // Check convergence
1197 if (convFugRatio.two_norm() < 1e-6) {
1198 // Restore cout format
1199 std::cout.flags(f);
1200
1201 // Print info
1202 if (verbosity >= 1) {
1203 std::cout << "Solution converged to the following result :" << std::endl;
1204 std::cout << "x = [";
1205 for (int compIdx=0; compIdx<numComponents; ++compIdx){
1206 if (compIdx < numComponents - 1)
1207 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) << " ";
1208 else
1209 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
1210 }
1211 std::cout << "]" << std::endl;
1212 std::cout << "y = [";
1213 for (int compIdx=0; compIdx<numComponents; ++compIdx){
1214 if (compIdx < numComponents - 1)
1215 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) << " ";
1216 else
1217 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
1218 }
1219 std::cout << "]" << std::endl;
1220 std::cout << "K = [" << K << "]" << std::endl;
1221 std::cout << "L = " << L << std::endl;
1222 }
1223 // Restore cout format format
1224 return;
1225 }
1226
1227 // If convergence is not met, K is updated in a successive substitution manner
1228 else{
1229 // Update K
1230 for (int compIdx=0; compIdx<numComponents; ++compIdx){
1231 K[compIdx] *= newFugRatio[compIdx];
1232 }
1233
1234 // Solve Rachford-Rice to get L from updated K
1235 L = solveRachfordRice_g_(K, z, 0);
1236 }
1237
1238 }
1239 if (!newton_afterwards) {
1240 throw std::runtime_error(
1241 "Successive substitution composition update did not converge within maxIterations " + std::to_string(maxIterations));
1242 }
1243 }
1244
1245};//end PTFlash
1246
1247} // namespace Opm
1248
1249#endif
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
A central place for various physical constants occuring in some equations.
Representation of an evaluation of a function and its derivatives w.r.t.
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...
Implements the Peng-Robinson equation of state for a mixture.
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
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
Determines the phase compositions, pressures and saturations given the total mass of all components f...
Definition PTFlash.hpp:63
static void solve(FluidState &fluid_state, const ComponentVector &globalMolarities, Scalar tolerance=0.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition PTFlash.hpp:194
static void solve(FluidState &fluid_state, const Dune::FieldVector< typename FluidState::Scalar, numComponents > &z, std::string twoPhaseMethod, Scalar tolerance=-1., int verbosity=0)
Calculates the fluid state from the global mole fractions of the components and the phase pressures.
Definition PTFlash.hpp:82
Implements the Peng-Robinson equation of state for a mixture.
Definition PengRobinsonMixture.hpp:41
static LhsEval computeFugacityCoefficient(const FluidState &fs, const Params &params, unsigned phaseIdx, unsigned compIdx)
Returns the fugacity coefficient of an individual component in the phase.
Definition PengRobinsonMixture.hpp:74
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30