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