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