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