OpenCAEPoro  0.2.0 Sep/22/2022
A simulator for multicomponent porous media flow
MixtureComp.hpp
Go to the documentation of this file.
1 
12 #ifndef __MIXTURECOMP_HEADER__
13 #define __MIXTURECOMP_HEADER__
14 
15 // Standard header files
16 #include <algorithm>
17 #include <math.h>
18 #include <vector>
19 
20 // OpenCAEPoro header files
21 #include "DenseMat.hpp"
22 #include "Mixture.hpp"
23 #include "OCPTable.hpp"
24 
25 using namespace std;
26 
29 {
30 public:
33  OCP_DBL Ktol{ 1E-4 };
34  OCP_DBL dYtol{ 1E-6 };
35  OCP_DBL eYt{ 1E-8 };
38  bool conflag;
39  // test
41  OCP_DBL curSk;
42 };
43 
46 {
47 public:
52  bool conflag;
53  // test
55 };
56 
57 
60 {
61 public:
66  bool conflag;
67  // test
69 };
70 
71 
73 class NRparamSP
74 {
75 public:
80  bool conflag;
81  // test
83 };
84 
86 class RRparam
87 {
88 public:
92  // test
94 };
95 
97 {
98  friend class MixtureComp;
99 
100 private:
101  SSMparamSTA SSMsta;
102  NRparamSTA NRsta;
103  SSMparamSP SSMsp;
104  NRparamSP NRsp;
105  RRparam RR;
106 };
107 
108 class COMP
109 {
110 public:
111  COMP() = default;
112  COMP(const vector<string>& comp);
113 
114 public:
115  string name;
125 };
126 
127 class MixtureComp : public Mixture
128 {
129 
130 public:
131  OCP_DBL GetErrorPEC() override { return ePEC; }
132  OCP_ULL GetSSMSTAiters() override { return SSMSTAiters; }
133  OCP_ULL GetNRSTAiters() override { return NRSTAiters; }
134  OCP_ULL GetSSMSPiters() override { return SSMSPiters; }
135  OCP_ULL GetNRSPiters() override { return NRSPiters; }
136  OCP_ULL GetRRiters() override { return RRiters; }
137  OCP_ULL GetSSMSTAcounts() override { return SSMSTAcounts; }
138  OCP_ULL GetNRSTAcounts() override { return NRSTAcounts; }
139  OCP_ULL GetSSMSPcounts() override { return SSMSPcounts; }
140  OCP_ULL GetNRSPcounts() override { return NRSPcounts; }
141  OCP_ULL GetRRcounts() override { return RRcounts; }
142 
143 private:
144  // total iters
145  OCP_ULL SSMSTAiters{ 0 };
146  OCP_ULL NRSTAiters{ 0 };
147  OCP_ULL SSMSPiters{ 0 };
148  OCP_ULL NRSPiters{ 0 };
149  OCP_ULL RRiters{ 0 };
150  // total counts, one count may contain many iters
151  OCP_ULL SSMSTAcounts{ 0 };
152  OCP_ULL NRSTAcounts{ 0 };
153  OCP_ULL SSMSPcounts{ 0 };
154  OCP_ULL NRSPcounts{ 0 };
155  OCP_ULL RRcounts{ 0 };
156  // phase equilibrium calculation error
157  // if NP = 1, it's from phase stable analysis, if skiped, it's 0
158  // if NP > 1, it's from phase spliting calculation
159  OCP_DBL ePEC;
160 
161 public:
162  MixtureComp() = default;
163  MixtureComp(const ParamReservoir& rs_param, const USI& i) :MixtureComp(rs_param.EoSp, i) {
164  mixtureType = EOS_PVTW;
165  if (rs_param.PVTW_T.data.size() != 0) {
166  PVTW.Setup(rs_param.PVTW_T.data[i]);
167  if (rs_param.gravity.activity)
168  std_RhoW = RHOW_STD * rs_param.gravity.data[1];
169  if (rs_param.density.activity)
170  std_RhoW = RHOW_STD;
171  std_GammaW = GRAVITY_FACTOR * std_RhoW;
172  data.resize(5);
173  cdata.resize(5);
174  }
175  };
176  MixtureComp(const EoSparam& param, const USI& i);
177  void InitFlash(const OCP_DBL& Pin, const OCP_DBL& Pbbin, const OCP_DBL& Tin,
178  const OCP_DBL* Sjin, const OCP_DBL& Vpore,
179  const OCP_DBL* Ziin) override;
180  void InitFlashDer(const OCP_DBL& Pin, const OCP_DBL& Pbbin, const OCP_DBL& Tin,
181  const OCP_DBL* Sjin, const OCP_DBL& Vpore,
182  const OCP_DBL* Ziin) override;
183  void InitFlashDer_n(const OCP_DBL& Pin, const OCP_DBL& Pbbin, const OCP_DBL& Tin,
184  const OCP_DBL* Sjin, const OCP_DBL& Vpore,
185  const OCP_DBL* Ziin) override;
186  // ftype = 0, flash from single phase
187  // ftype = 1, skip phase stablity analysis and num of phase = 1
188  // ftype = 1, skip phase stablity analysis and num of phase = 2
189  void Flash(const OCP_DBL& Pin, const OCP_DBL& Tin, const OCP_DBL* Niin, const USI& ftype, const USI& lastNP,
190  const OCP_DBL* lastKs) override;
191  void CalFlash(const OCP_DBL& Pin, const OCP_DBL& Tin, const OCP_DBL* Niin);
192  void FlashDeriv(const OCP_DBL& Pin, const OCP_DBL& Tin,
193  const OCP_DBL* Niin, const USI& ftype, const USI& lastNP,
194  const OCP_DBL* lastKs) override;
195  void FlashDeriv_n(const OCP_DBL& Pin, const OCP_DBL& Tin,
196  const OCP_DBL* Niin, const OCP_DBL* Sjin, const OCP_DBL* xijin,
197  const OCP_DBL* njin, const USI& ftype, const USI* phaseExistin,
198  const USI& lastNP, const OCP_DBL* lastKs)override;
199  OCP_DBL XiPhase(const OCP_DBL& Pin, const OCP_DBL& Tin, const OCP_DBL* Ziin) override;
200  OCP_DBL RhoPhase(const OCP_DBL& Pin, const OCP_DBL& Tin,
201  const OCP_DBL* Ziin) override;
202  OCP_DBL GammaPhaseO(const OCP_DBL& Pin, const OCP_DBL& Pbbin) override { OCP_ABORT("Should not be used in Compositional mode!"); };
203  OCP_DBL GammaPhaseG(const OCP_DBL& Pin) override { OCP_ABORT("Should not be used in Compositional mode!"); };
204  OCP_DBL GammaPhaseW(const OCP_DBL& Pin) override;
205  OCP_DBL GammaPhaseOG(const OCP_DBL& Pin, const OCP_DBL& Tin,
206  const OCP_DBL* Ziin) override ;
207  void setPT(const OCP_DBL& p, const OCP_DBL& t) { P = p; T = t; }
208  void setZi(const OCP_DBL* Ziin) { Dcopy(NC, &zi[0], Ziin); }
209  void setZi() { for (USI i = 0; i < NC; i++) zi[i] = Ni[i] / Nh; }
210  void setNi(const OCP_DBL* Niin) { Dcopy(numCom, &Ni[0], Niin); }
211  void CallId();
212  USI GetFtype() override { return ftype; }
213  void CalSurfaceTension();
214  OCP_DBL GetSurTen() override { return surTen; }
215 
216 private:
217 
218  // Basic Components Informations
219  vector<string> Cname;
220  vector<OCP_DBL> Tc;
221  vector<OCP_DBL> Pc;
222  vector<OCP_DBL> Vc;
223  vector<OCP_DBL> MWC;
224  vector<OCP_DBL> Acf;
225  vector<OCP_DBL> OmegaA;
226  vector<OCP_DBL> OmegaB;
227  vector<OCP_DBL> Vshift;
228  vector<OCP_DBL> Parachor;
229  vector<OCP_DBL> Zc;
230  bool ParachorAct;
231  // for viscosity calculation
232  vector<OCP_DBL> Vcvis;
233  vector<OCP_DBL> Zcvis;
234  vector<OCP_DBL> LBCcoef;
235  vector<OCP_DBL> BIC;
236 
237  // Model information
238  bool miscible;
239  OCP_DBL surTen;
240 
241  // Initial properties
242  USI NC;
243  USI NPmax;
244  OCP_DBL P;
245  OCP_DBL T;
246  vector<OCP_DBL> zi;
247  vector<COMP> comp;
248  USI lId;
249  EoScontrol EoSctrl;
250  USI ftype{ 0 };
251  USI tmpFtype;
252 
253  vector<OCP_DBL> Plist;
254  vector<OCP_DBL> Tlist;
255  vector<OCP_DBL> Ytlist;
256 
257 private:
258  OCPTable PVTW;
259  OCP_DBL std_RhoW;
260  OCP_DBL std_GammaW;
261  vector<OCP_DBL> data;
263  vector<OCP_DBL> cdata;
265 
266 public:
267  // EoS Function
268  // Allocate memory for EoS
269  void AllocateEoS();
270  // Setup and Solve EoS for specified A,B
271  void SolEoS(OCP_DBL& ZjT, const OCP_DBL& AjT, const OCP_DBL& BjT) const;
272  // Calculate Ai and Bi
273  void CalAiBi();
274  // Calculate Aj and Bj with specified xj
275  void CalAjBj(OCP_DBL& AjT, OCP_DBL& BjT, const vector<OCP_DBL>& xj) const;
276  void CalAjBj(OCP_DBL& AjT, OCP_DBL& BjT, const OCP_DBL* xj) const;
278  USI CubicRoot(const OCP_DBL& a, const OCP_DBL& b, const OCP_DBL& c,
279  const bool& NTflag = false) const;
281  void PrintZtmp()
282  {
283  for (auto z : Ztmp) cout << z << "\t";
284  }
285 
286 private:
287  // EoS Variables
288  vector<OCP_DBL> Ai;
289  vector<OCP_DBL> Bi;
290  vector<OCP_DBL> Aj;
291  vector<OCP_DBL> Bj;
292  vector<OCP_DBL> Zj;
293  mutable vector<OCP_DBL> Ztmp;
294 
295  // PR default
296  OCP_DBL delta1 = 2.41421356237;
297  OCP_DBL delta2 = -0.41421356237;
298  OCP_DBL delta1P2;
299  OCP_DBL delta1M2;
300  OCP_DBL delta1T2;
301 
302 public:
303  // Phase Function
304  // Allocate memoery for phase variables
305  void AllocatePhase();
306  void CalFugPhi(vector<OCP_DBL>& phiT, vector<OCP_DBL>& fugT,
307  const vector<OCP_DBL>& xj);
308  void CalFugPhi(OCP_DBL* phiT, OCP_DBL* fugT, const OCP_DBL* xj);
309  void CalFugPhi(OCP_DBL* fugT, const OCP_DBL* xj);
310  void CalFugPhiAll();
311  void CalMW();
312  void CalVfXiRho();
313  void CalSaturation();
314  USI FindMWmax();
315  void x2n();
316  void PrintX();
317 
318 private:
319  // Phase Variables
320  USI lNP{ 0 };
321  USI NP;
322  USI inputNP;
323  OCP_DBL Nh;
324  vector<OCP_DBL> vC;
325  vector<OCP_DBL> nu;
326  vector<vector<OCP_DBL>> x;
327  vector<vector<OCP_DBL>> phi;
328  vector<vector<OCP_DBL>> fug;
329  vector<vector<OCP_DBL>> n;
330  vector<vector<OCP_DBL>> ln;
331  vector<OCP_DBL> xiC;
332  vector<OCP_DBL> rhoC;
333  vector<OCP_DBL> MW;
334  vector<USI> phaseLabel;
335 
336 public:
337  // Method Function
338  // Allocate memoery for Method variables
339  void AllocateMethod();
340  void PhaseEquilibrium();
341  void CalKwilson();
342  bool PhaseStable();
343  bool StableSSM(const USI& Id);
344  bool StableSSM01(const USI& Id);
345  bool StableNR(const USI& Id);
346  void CalFugXSTA();
347  void AssembleJmatSTA();
348  bool CheckSplit();
349  void PhaseSplit();
350  void SplitSSM(const bool& flag);
351  void SplitSSM2(const bool& flag);
352  void SplitSSM3(const bool& flag);
353  void RachfordRice2();
354  void RachfordRice2P();
355  void RachfordRice3();
356  void UpdateXRR();
357  void SplitBFGS();
358  void SplitNR();
359  void CalResSP();
360  void CalFugNAll(const bool& Znflag = true);
361  void PrintFugN();
362  void AssembleJmatSP();
364  void CalPhiNSTA();
365  void AssembleSkipMatSTA();
366  OCP_DBL CalStepNRsp();
367 
368 
369  OCP_SIN GetMinEigenSkip() override { return eigenSkip[0]; }
370  bool GetFlagSkip() override { return flagSkip; }
371 
372 private:
373  // Method Variables
374  USI testPId;
375  vector<vector<OCP_DBL>> Kw;
376  vector<vector<OCP_DBL>> Ks;
377  vector<OCP_DBL> lKs;
378  OCP_DBL Asta, Bsta, Zsta;
379  vector<OCP_DBL> phiSta;
380  vector<OCP_DBL> fugSta;
381  // SSM in Stability Analysis
382  vector<OCP_DBL> Y;
383  OCP_DBL Yt;
384  vector<OCP_DBL> di;
385  // NR in Stability Analysis
386  vector<OCP_DBL> resSTA;
387  vector<OCP_DBL> JmatSTA;
388  vector<vector<OCP_DBL>> fugX;
389  vector<OCP_DBL> Ax;
390  vector<OCP_DBL> Bx;
391  vector<OCP_DBL> Zx;
392  // Skip Stability Analysis
393  bool flagSkip;
394  vector<OCP_DBL> phiN;
395  vector<OCP_SIN> skipMatSTA;
396  vector<OCP_SIN> eigenSkip;
397  vector<OCP_SIN> eigenWork;
398  OCP_INT leigenWork;
399 
400  // SSM in Phase Split
401  vector<OCP_DBL> resRR;
402  // NR in Phase Split
403  vector<OCP_DBL> lresSP;
404  vector<OCP_DBL> resSP;
405  vector<OCP_DBL> JmatSP;
406  vector<vector<OCP_DBL>>
407  fugN;
408  vector<OCP_DBL> An;
409  vector<OCP_DBL> Bn;
410  vector<vector<OCP_DBL>> Zn;
411  // for linearsolve with lapack
412  vector<OCP_INT> pivot;
413  vector<OCP_DBL> JmatWork;
414  OCP_INT lJmatWork;
415  char uplo{'U'};
416 
417 public:
418  // After Phase Equilibrium Calculation finishs, properties and some auxiliary
419  // variables will be calculated.
420  void AllocateOthers();
421  void IdentifyPhase();
423  void CopyPhase();
424  void CalViscosity();
425  void CalViscoLBC();
426  void CalViscoHZYT();
427  void CalFugXAll();
428  void CalFugPAll(const bool& Zpflag = true);
429 
430 
431  void CalVjpVfpVfx_partial();
432  void CalXiPNX_partial();
433  void CalRhoPX_partial();
434  void CalMuPX_partial();
435  void CalMuPXLBC_partial();
436  void CalXiRhoMuPN_pfullx();
437  void CaldXsdXpAPI04();
438  void CaldXsdXp04();
439 
440  void CalRhoPNX_full();
441 
442  void CalXiPNX_full01();
443  void CalRhoPNX_full01();
444  void CalMuPX_full01();
445  void CalMuPXLBC_full01();
446  void CalVfiVfp_full01();
447  void AssembleMatVfiVfp_full01();
448  void AssembleRhsVfiVfp_full01();
449  void CaldXsdXp01();
450  void CaldXsdXpAPI01();
451 
452  void CalXiPNX_full02();
453  void CalVfiVfp_full02();
454  void AssembleMatVfiVfp_full02();
455  void AssembleRhsVfiVfp_full02();
456  void CaldXsdXpAPI02();
457  void CaldXsdXpAPI02p();
458 
459  void CalVjpVfpVfn_partial();
460  void CalXiPn_partial();
461  void CalRhoPn_partial();
462  void CalMuPn_partial();
463  void CalMuPnLBC_partial();
464  void CalXiRhoMuPN_pfullxn(const bool& xflag = true);
465 
466  void CaldXsdXpAPI03();
467  void CaldXsdXp03();
468  void CalVfiVfp_full03();
469 
470  void CalKeyDerx();
471  void CalKeyDern();
472 
473 private:
474  // Phase properties and auxiliary variables
475  vector<OCP_DBL> muC;
476  vector<vector<OCP_DBL>>
477  muAux;
478  vector<OCP_DBL>
479  muAux1I;
480  vector<OCP_DBL> sqrtMWi;
481  vector<vector<OCP_DBL>> fugP;
482  vector<OCP_DBL> Zp;
483 
484  vector<OCP_DBL> JmatTmp;
485  vector<OCP_DBL> JmatDer;
488  vector<OCP_DBL> rhsDer;
489  vector<OCP_DBL> xixC;
490  vector<OCP_DBL> xiPC;
491  vector<OCP_DBL> xiNC;
492 };
493 
495 OCP_DBL signD(const OCP_DBL& d);
496 
497 OCP_DBL delta(const USI& i, const USI& j);
498 
499 void NTcubicroot(OCP_DBL& root, const OCP_DBL& a, const OCP_DBL& b, const OCP_DBL& c);
500 
501 #endif
Operations about small dense mat.
void Dcopy(const int &N, double *dst, const double *src)
Calculate the minimal eigenvalue for sysmetric matrix with mkl lapack.
Definition: DenseMat.cpp:37
OCP_DBL signD(const OCP_DBL &d)
Return the sign of double di.
Mixture class declaration.
const USI EOS_PVTW
Mixture model = equation-of-state.
Definition: OCPConst.hpp:89
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
const OCP_DBL GRAVITY_FACTOR
0.00694444 ft2 psi / lb
Definition: OCPConst.hpp:52
float OCP_SIN
Single precision.
Definition: OCPConst.hpp:27
const OCP_DBL RHOW_STD
Water density in lb / ft3.
Definition: OCPConst.hpp:53
int OCP_INT
Long integer.
Definition: OCPConst.hpp:25
unsigned long long OCP_ULL
Long long unsigned integer.
Definition: OCPConst.hpp:23
OCPTable class declaration.
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47
OCP_DBL MW
Molecular Weight.
OCP_DBL OmegaA
Param A of Components.
OCP_DBL Vc
VcMW * MW.
OCP_DBL VcMW
Critical Volume / MW.
OCP_DBL acf
Acentric Factor.
OCP_DBL Tc
Critical Temperature.
OCP_DBL OmegaB
Param B of Components.
OCP_DBL Vshift
shift volume
OCP_DBL Pc
Critical Pressure.
string name
Name of components.
EoSParam contains the params for Compositional Model and functions to read them.
OCP_DBL GammaPhaseG(const OCP_DBL &Pin) override
return gamma of gas phase, gamma equals to mass density times gravity factor.
OCP_DBL GammaPhaseO(const OCP_DBL &Pin, const OCP_DBL &Pbbin) override
return gamma of oil phase, gamma equals to mass density times gravity factor.
void PrintZtmp()
test
Params for NR in Phase Split.
Definition: MixtureComp.hpp:74
USI curIt
current Iters
Definition: MixtureComp.hpp:82
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:78
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:76
bool conflag
convergence flag, if converges, conflag = true
Definition: MixtureComp.hpp:80
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:77
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:79
Params for NR in Phase Stability Analysis.
Definition: MixtureComp.hpp:46
bool conflag
convergence flag, if converges, conflag = true
Definition: MixtureComp.hpp:52
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:49
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:51
USI curIt
current Iters
Definition: MixtureComp.hpp:54
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:48
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:50
Type_A_r< OCP_DBL > density
Density of oil, water, gas in standard conditions.
Type_A_r< OCP_DBL > gravity
Gravity of oil, water, gas in standard conditions.
TableSet PVTW_T
Table set of PVTW.
EoSparam EoSp
Initial component composition, used in compositional models.
Param for Solving Rachford-Rice Equations.
Definition: MixtureComp.hpp:87
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:91
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:89
USI curIt
current Iters
Definition: MixtureComp.hpp:93
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:90
Params for SSM in Phase Split.
Definition: MixtureComp.hpp:60
bool conflag
convergence flag, if converges, conflag = true
Definition: MixtureComp.hpp:66
USI curIt
current Iters
Definition: MixtureComp.hpp:68
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:63
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:64
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:65
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:62
Params for SSM in Phase Stability Analysis.
Definition: MixtureComp.hpp:29
OCP_DBL realTol
Real tol.
Definition: MixtureComp.hpp:37
USI curIt
current Iterations
Definition: MixtureComp.hpp:40
OCP_DBL tol2
tol*tol
Definition: MixtureComp.hpp:36
bool conflag
convergence flag, if converges, conflag = true
Definition: MixtureComp.hpp:38
OCP_DBL tol
Tolerance.
Definition: MixtureComp.hpp:32
USI maxIt
Max Iteration.
Definition: MixtureComp.hpp:31
vector< vector< vector< OCP_DBL > > > data
All table with the same name.
vector< T > data
Data of param.
bool activity
If false, this param is not given.