OpenCAEPoro  0.2.0 Sep/22/2022
A simulator for multicomponent porous media flow
Bulk.hpp
Go to the documentation of this file.
1 
12 #ifndef __BULK_HEADER__
13 #define __BULK_HEADER__
14 
15 // Standard header files
16 #include <iostream>
17 #include <vector>
18 
19 // OpenCAEPoro header files
20 #include "DenseMat.hpp"
21 #include "FlowUnit.hpp"
22 #include "Grid.hpp"
23 #include "LinearSystem.hpp"
24 #include "Mixture.hpp"
25 #include "MixtureBO.hpp"
26 #include "MixtureComp.hpp"
27 #include "OCPConst.hpp"
28 #include "OCPStructure.hpp"
29 #include "ParamReservoir.hpp"
30 
31 using namespace std;
32 
34 // Note: This class includes reference depth and pressure at it, depth of contacts
35 // between phases, and capillary pressure at phase contact surfaces.
36 class ParamEQUIL
37 {
38  friend class Bulk;
39 
40 private:
41  OCP_DBL Dref;
42  OCP_DBL Pref;
43  OCP_DBL DOWC;
44  OCP_DBL DGOC;
45  OCP_DBL PcOW;
46  OCP_DBL PcGO;
47 
48  OCPTable PBVD;
49 };
50 
52 // Note: Bulk contains main physical infomation of active grids. It describes the
53 // actural geometric domain for simulating. Variables are stored bulk by bulk, and then
54 // phase by phase, then component by component. The bulks are ordered in the alphabetic
55 // order, i.e. the X-axis indices first, followed by the Y-axis and Z-axis indices.
56 // Operations on each bulk are also defined here.
57 class Bulk
58 {
59  friend class BulkConn;
60  friend class Well;
61  friend class DetailInfo;
62 
63  // temp
64  friend class OCP_IMPEC;
65  friend class OCP_FIM;
66  friend class Reservoir;
67  friend class OCP_AIMt;
68  friend class OCP_AIMc;
69 
71  // For general usage
73 
74 public:
75  Bulk() = default;
76 
78  void InputParam(ParamReservoir& rs_param);
80  void Setup(const Grid& myGrid);
82  void InitSjPcBo(const USI& tabrow);
84  void InitSjPcComp(const USI& tabrow, const Grid& myGrid);
86  void InitFlash(const bool& flag = false);
88  void InitFlashDer();
89  void InitFlashDer_n();
91  void Flash();
93  void FlashBLKOIL();
95  void FlashCOMP();
97  void FlashDeriv();
98  void FlashDeriv_n();
100  void FlashDerivBLKOIL();
101  void FlashDerivBLKOIL_n();
103  void FlashDerivCOMP();
104  void FlashDerivCOMP_n();
106  USI CalFlashType(const OCP_USI& n) const;
108  void PassFlashValue(const OCP_USI& n);
109  void PassFlashValueAIMc(const OCP_USI& n);
111  void PassFlashValueDeriv(const OCP_USI& n);
112  void PassFlashValueDeriv_n(const OCP_USI& n);
114  void ResetFlash();
115 
117  void CalKrPc();
119  void CalKrPcDeriv();
121  void CalVpore();
123  OCP_DBL CalFPR() const;
125  void CalMaxChange();
126 
128  void CheckSetup() const;
130  void CheckInitVpore() const;
132  void CheckVpore() const;
134  bool CheckP() const;
136  bool CheckNi();
138  bool CheckVe(const OCP_DBL& Vlim) const;
140  void CheckSat() const;
142  void CheckDiff();
143 
145  OCP_USI GetBulkNum() const { return numBulk; }
147  USI GetComNum() const { return numCom; }
149  USI GetPhaseNum() const { return numPhase; }
151  USI GetMixMode() const;
153  const vector<Mixture*>& GetMixture() const { return flashCal; }
155  OCP_DBL GetP(const OCP_USI& n) const { return P[n]; }
157  OCP_DBL GetSOIL(const OCP_USI& n) const
158  {
159  return S[n * numPhase + phase2Index[OIL]];
160  }
162  OCP_DBL GetSGAS(const OCP_USI& n) const
163  {
164  return S[n * numPhase + phase2Index[GAS]];
165  }
167  OCP_DBL GetSWAT(const OCP_USI& n) const
168  {
169  return S[n * numPhase + phase2Index[WATER]];
170  }
172  OCP_DBL GetdPmax() const { return dPmax; }
174  OCP_DBL GetdNmax() const { return dNmax; }
176  OCP_DBL GetdSmax() const { return dSmax; }
178  OCP_DBL GetdVmax() const { return dVmax; }
179 
180  // Reset phaseNum to the ones of the last time step.
181  void ResetphaseNum() { phaseNum = lphaseNum; }
182  // Reset minEigenSkip to the ones of the last time step.
183  void ResetminEigenSkip() { minEigenSkip = lminEigenSkip; }
184  // Reset flagSkip to the ones of the last time step.
185  void ResetflagSkip() { flagSkip = lflagSkip; }
186  // Reset flagSkip to the ones of the last time step.
187  void ResetziSkip() { ziSkip = lziSkip; }
188  // Reset flagSkip to the ones of the last time step.
189  void ResetPSkip() { PSkip = lPSkip; }
190  // Reser Ks to the ones of the last time step.
191  void ResetKs() { Ks = lKs; }
192 
194  void ResetNt() { Nt = lNt; }
195 
197  void ResetP() { P = lP; }
199  void ResetPj() { Pj = lPj; }
201  void ResetNi() { Ni = lNi; }
203  void ResetKr() { kr = lkr; }
205  void ResetVp() { rockVp = lrockVp; }
206  void CalSomeInfo(const Grid& myGrid) const;
207 
209  void AllocateWellBulkId(const USI& n) { wellBulkId.reserve(n); }
210  void ClearWellBulkId() { wellBulkId.clear(); }
211  OCP_DBL CalNT() { NT = Dnorm1(numBulk, &Nt[0]); return NT; }
212 
213 
214 private:
216  // General variables
218  OCP_USI numBulk;
219  USI numPhase;
220  USI numCom;
221  USI numCom_1;
222 
223 
224  // Initial proportion of each component for EoS : numCom - 1, water is excluded.
225  vector<OCP_DBL> initZi;
226  vector<OCPTable> initZi_T;
227 
228  vector<OCP_DBL> SwatInit;
229  bool SwatInitExist{false};
230  vector<OCP_DBL> ScaleValuePcow;
231  bool ScalePcow{false};
232 
233 
234  USI PVTmode;
235  vector<USI> PVTNUM;
236  USI NTPVT;
237  vector<Mixture*> flashCal;
238  USI SATmode;
239  vector<USI> SATNUM;
240  USI NTSFUN;
241  vector<FlowUnit*> flow;
242  vector<vector<OCP_DBL>> satcm;
243 
244  // Skip stability analysis
245  vector<USI> phaseNum;
246  vector<USI> lphaseNum;
247  vector<USI> NRphaseNum;
249  vector<OCP_SIN> minEigenSkip;
250  vector<bool> flagSkip;
251  vector<OCP_DBL> ziSkip;
252  vector<OCP_DBL> PSkip;
253  vector<OCP_SIN> lminEigenSkip;
254  vector<bool> lflagSkip;
255  vector<OCP_DBL> lziSkip;
256  vector<OCP_DBL> lPSkip;
257 
258  // phase split calculation
259  // IMPORTANT!!!
260  // Ks will change as long as nums of hydroncarbon phase equals 2, and it will has an important effect
261  // on phase spliting calculation as a intial value. So you should not expect to obtain
262  // the exact same result with identical P, T, Ni if the final mixture contains 2 hydroncarbon phase.
263  vector<OCP_DBL> Ks;
264  vector<OCP_DBL> lKs;
265 
267  // Basic model information
269  ParamEQUIL EQUIL;
270  bool blackOil;
271  bool comps;
272  bool oil;
273  bool gas;
274  bool water;
275  bool disGas;
276  bool miscible;
277 
279  // Basic physical variables
281  OCP_DBL T;
282  vector<OCP_DBL> Pb;
283  vector<OCP_DBL> P;
284  vector<OCP_DBL> Pj;
285  vector<OCP_DBL> Pc;
286  vector<bool> phaseExist;
287  vector<OCP_DBL> S;
288  vector<OCP_DBL> nj;
289  vector<OCP_DBL> rho;
290  vector<OCP_DBL> xi;
291  vector<OCP_DBL> xij;
292  vector<OCP_DBL> Ni;
293  vector<OCP_DBL> mu;
294  vector<OCP_DBL> kr;
295  vector<OCP_DBL> vj;
296  vector<OCP_DBL> vf;
297  vector<OCP_DBL> Nt;
298  OCP_DBL NT;
299  // Note: Nij is the moles of component i in phase j, Nj is the moles of phase j.
300 
302  // Derivatives
304  vector<OCP_DBL> vfi;
305  vector<OCP_DBL> vfp;
306 
308  // Properties at the last time step, determined by methods.
310  vector<OCP_DBL> lP;
311  vector<OCP_DBL> lPj;
312  vector<OCP_DBL> lPc;
313  vector<bool> lphaseExist;
314  vector<OCP_DBL> lS;
315  vector<OCP_DBL> lnj;
316  vector<OCP_DBL> lrho;
317  vector<OCP_DBL> lxi;
318  vector<OCP_DBL> lxij;
319  vector<OCP_DBL> lNi;
320  vector<OCP_DBL> lmu;
321  vector<OCP_DBL> lkr;
322  vector<OCP_DBL> lvj;
323  vector<OCP_DBL> lvf;
324  vector<OCP_DBL> lNt;
325  vector<OCP_DBL> lvfi;
326  vector<OCP_DBL> lvfp;
327  vector<OCP_DBL> lrockVp;
328 
329  vector<OCP_DBL> surTen;
330  vector<OCP_DBL> Fk;
331  vector<OCP_DBL> Fp;
332 
333  vector<OCP_DBL> lsurTen;
334 
336  // Reservoir rock infomation of each bulk (size = numBulk)
338  vector<OCP_DBL> dx;
339  vector<OCP_DBL> dy;
340  vector<OCP_DBL> dz;
341  vector<OCP_DBL> depth;
342  vector<OCP_DBL> ntg;
343  vector<OCP_DBL> rockVpInit;
344  vector<OCP_DBL> rockVp;
345  OCP_DBL rockPref;
346  OCP_DBL rockC1;
347  OCP_DBL rockC2;
348  vector<OCP_DBL> rockKxInit;
349  vector<OCP_DBL> rockKx;
350  vector<OCP_DBL> rockKyInit;
351  vector<OCP_DBL> rockKy;
352  vector<OCP_DBL> rockKzInit;
353  vector<OCP_DBL> rockKz;
354 
356  // Auxiliary variables
358  vector<USI> index2Phase;
359  // Note: For example, phase 0 is Oil.
360  vector<USI> phase2Index;
361  // Note: For example, `Oil' is at the i-th location.
362 
364  // Max range for some variables at the current time step
366  OCP_DBL dPmax;
367  OCP_DBL dSmax;
368  OCP_DBL dNmax;
369  OCP_DBL dVmax;
370 
372  // Error
374 
375  vector<OCP_DBL> ePEC;
376  mutable vector<OCP_DBL> eN;
377  mutable vector<OCP_DBL> eV;
378 
380  // For IMPEC
382 
383 public:
385  void AllocateAuxIMPEC();
387  void GetSolIMPEC(const vector<OCP_DBL>& u);
389  void SetCFL2Zero() const { fill(cfl.begin(), cfl.end(), 0); }
391  OCP_DBL CalCFL() const;
393  void UpdateLastStepIMPEC();
394 
395 private:
396  mutable vector<OCP_DBL> cfl;
397 
399  // For FIM
401 
402 public:
404  void AllocateAuxFIM();
406  void GetSolFIM(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
407  const OCP_DBL& dSmaxlim);
408  void GetSolFIM_n(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
409  const OCP_DBL& dSmaxlim);
411  void GetSol01FIM(const vector<OCP_DBL>& u);
413  void CalRelResFIM(ResFIM& resFIM) const;
414  // Show Res
415  void ShowRes(const vector<OCP_DBL>& res) const;
417  void ResetFIM();
419  void UpdateLastStepFIM();
421  OCP_DBL CalNRdSmax(OCP_USI& index);
422 
424  OCP_DBL GetNRdPmax();
426  OCP_DBL GetNRdSmaxP();
427  OCP_DBL GetNRdNmax();
428  void CorrectNi(const vector<OCP_DBL>& res);
429 
430 
431 private:
432  // Derivatives for FIM
433  vector<OCP_DBL> muP;
434  vector<OCP_DBL> xiP;
435  vector<OCP_DBL> rhoP;
436  vector<OCP_DBL> mux;
437  vector<OCP_DBL> xix;
438  vector<OCP_DBL> rhox;
439  vector<OCP_DBL> dPcj_dS;
440  vector<OCP_DBL> dKr_dS;
441  vector<OCP_DBL> dSec_dPri;
442  vector<OCP_DBL> res_n;
443  vector<OCP_DBL> resPc;
444  USI lendSdP;
445  vector<OCP_USI> dSdPindex;
446  vector<OCP_USI> resIndex;
447  // if phase is inexisting or num of components <= 1, set it to zero
448  // it is an auxiliary variable of dSec_dPri
449  vector<USI> pEnumCom;
450 
451  // vars at last step
452  vector<OCP_DBL> lmuP;
453  vector<OCP_DBL> lxiP;
454  vector<OCP_DBL> lrhoP;
455  vector<OCP_DBL> lmux;
456  vector<OCP_DBL> lxix;
457  vector<OCP_DBL> lrhox;
458  vector<OCP_DBL> ldPcj_dS;
459  vector<OCP_DBL> ldKr_dS;
460  vector<OCP_DBL> ldSec_dPri;
461  vector<OCP_DBL> lres_n;
462  vector<OCP_DBL> lresPc;
463  vector<OCP_USI> ldSdPindex;
464  vector<OCP_USI> lresIndex;
465  vector<USI> lpEnumCom;
466 
467 
468  vector<OCP_DBL> dSNR;
469  vector<OCP_DBL> dSNRP;
470  vector<OCP_DBL> dNNR;
471  vector<OCP_DBL> dPNR;
472 
473  OCP_DBL NRdSSP;
474  OCP_DBL maxNRdSSP;
475  OCP_USI index_maxNRdSSP;
476  OCP_DBL NRdPmax;
477  OCP_DBL NRdNmax;
478  OCP_DBL NRdSmax;
479  OCP_DBL NRdSmaxP;
480 
481 
482  vector<OCP_DBL> NRstep;
483 
484 public:
485  // for debug!
486  void OutputInfo(const OCP_USI& n) const;
487  OCP_ULL GetSSMSTAiters()const { return flashCal[0]->GetSSMSTAiters(); }
488  OCP_ULL GetNRSTAiters()const { return flashCal[0]->GetNRSTAiters(); }
489  OCP_ULL GetSSMSPiters()const { return flashCal[0]->GetSSMSPiters(); }
490  OCP_ULL GetNRSPiters()const { return flashCal[0]->GetNRSPiters(); }
491  OCP_ULL GetRRiters()const { return flashCal[0]->GetRRiters(); }
492  OCP_ULL GetSSMSTAcounts()const { return flashCal[0]->GetSSMSTAcounts(); }
493  OCP_ULL GetNRSTAcounts()const { return flashCal[0]->GetNRSTAcounts(); }
494  OCP_ULL GetSSMSPcounts()const { return flashCal[0]->GetSSMSPcounts(); }
495  OCP_ULL GetNRSPcounts()const { return flashCal[0]->GetNRSPcounts(); }
496  OCP_ULL GetRRcounts()const { return flashCal[0]->GetRRcounts(); }
497 
498 
500  // For AIMs, AIMt
502 
503 public:
505  void AllocateAuxAIM(const OCP_DBL& ratio);
506  OCP_USI GetMaxFIMBulk()const { return maxNumFIMBulk; }
508  void FlashDerivAIM(const bool& IfAIMs);
509  void PassFlashValueDerivAIM(const OCP_USI& n);
511  void CalKrPcDerivAIM(const bool& IfAIMs);
513  void CalRelResAIMt(ResFIM& resFIM) const;
515  void GetSolAIMt(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
516  const OCP_DBL& dSmaxlim);
517 
519  void CalRelResAIMs(ResFIM& resFIM) const;
520  void GetSolAIMs(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
521  const OCP_DBL& dSmaxlim);
522  void UpdateLastStepAIM();
523  void ResetFIMBulk();
524  void ShowFIMBulk(const bool& flag = false) const;
526  bool CheckNiFIMBulk() const;
528  void InFIMNi();
530  void OutFIMNi();
531 
533  // For AIMc
535 
536  void AllocateAuxAIMc();
538  void FlashAIMc();
540  void FlashBLKOILAIMc();
542  void FlashCOMPAIMc();
543 
544  void FlashAIMc01();
545  void FlashBLKOILAIMc01();
546  void FlashCOMPAIMc01();
547 
549  void FlashDerivAIMc();
551  void FlashDerivBLKOILAIMc();
553  void FlashDerivCOMPAIMc();
555  void CalKrPcAIMc();
557  void CalKrPcDerivAIMc();
558  void GetSolAIMc(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
559  const OCP_DBL& dSmaxlim);
560  void GetSolAIMc01(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
561  const OCP_DBL& dSmaxlim);
562  void UpdatePj();
563 
564 private:
565  vector<OCP_USI> wellBulkId;
566  vector<OCP_INT> map_Bulk2FIM;
567  vector<OCP_USI> FIMBulk;
568  OCP_USI numFIMBulk;
569  OCP_USI maxNumFIMBulk;
570  vector<OCP_DBL> FIMNi;
571 };
572 
573 #endif /* end if __BULK_HEADER__ */
574 
575 /*----------------------------------------------------------------------------*/
576 /* Brief Change History of This File */
577 /*----------------------------------------------------------------------------*/
578 /* Author Date Actions */
579 /*----------------------------------------------------------------------------*/
580 /* Shizhe Li Oct/01/2021 Create file */
581 /* Chensong Zhang Oct/15/2021 Format file */
582 /* Chensong Zhang Jan/09/2022 Update Doxygen */
583 /*----------------------------------------------------------------------------*/
DenseMat.hpp
Operations about small dense mat.
Dnorm1
double Dnorm1(const int &N, double *x)
Computes the L1-norm of a vector.
Definition: DenseMat.cpp:50
FlowUnit.hpp
FlowUnit class declaration.
Grid.hpp
Grid class declaration.
LinearSystem.hpp
Linear solver class declaration.
MixtureBO.hpp
MixtureBO class declaration.
MixtureComp.hpp
MixtureComp class declaration.
Mixture.hpp
Mixture class declaration.
OCPConst.hpp
Definition of build-in datatypes and consts.
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
GAS
const USI GAS
Fluid type = gas.
Definition: OCPConst.hpp:83
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
WATER
const USI WATER
Fluid type = water.
Definition: OCPConst.hpp:84
OIL
const USI OIL
Fluid type = oil.
Definition: OCPConst.hpp:82
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:24
OCP_ULL
unsigned long long OCP_ULL
Long long unsigned integer.
Definition: OCPConst.hpp:23
OCPStructure.hpp
Some Structure in OpenCAEPoro.
ParamReservoir.hpp
ParamReservoir class declaration.
BulkConn
Properties and operations on connections between bulks (active grids).
Definition: BulkConn.hpp:58
Bulk
Physical information of each active reservoir bulk.
Definition: Bulk.hpp:58
Bulk::GetPhaseNum
USI GetPhaseNum() const
Return the number of phases.
Definition: Bulk.hpp:149
Bulk::AllocateWellBulkId
void AllocateWellBulkId(const USI &n)
Allocate memory for WellbulkId.
Definition: Bulk.hpp:209
Bulk::ResetPj
void ResetPj()
Reset Pj to the ones of the last time step.
Definition: Bulk.hpp:199
Bulk::ResetNi
void ResetNi()
Reset Ni to the ones of the last time step.
Definition: Bulk.hpp:201
Bulk::ResetVp
void ResetVp()
Reset Vp to the ones of the last time step.
Definition: Bulk.hpp:205
Bulk::GetdVmax
OCP_DBL GetdVmax() const
Return dVmax.
Definition: Bulk.hpp:178
Bulk::ResetNt
void ResetNt()
Reset Nt to the ones of the last time step.
Definition: Bulk.hpp:194
Bulk::SetCFL2Zero
void SetCFL2Zero() const
Initialize the CFL number.
Definition: Bulk.hpp:389
Bulk::GetSOIL
OCP_DBL GetSOIL(const OCP_USI &n) const
Return oil saturation of the n-th bulk.
Definition: Bulk.hpp:157
Bulk::GetMixture
const vector< Mixture * > & GetMixture() const
Return flash results (it has not been used by far).
Definition: Bulk.hpp:153
Bulk::GetSWAT
OCP_DBL GetSWAT(const OCP_USI &n) const
Return water saturation of the n-th bulk.
Definition: Bulk.hpp:167
Bulk::GetSGAS
OCP_DBL GetSGAS(const OCP_USI &n) const
Return gas saturation of the n-th bulk.
Definition: Bulk.hpp:162
Bulk::GetdSmax
OCP_DBL GetdSmax() const
Return dSmax.
Definition: Bulk.hpp:176
Bulk::ResetP
void ResetP()
Reset P to the ones of the last time step.
Definition: Bulk.hpp:197
Bulk::ResetKr
void ResetKr()
Reset Kr to the ones of the last time step.
Definition: Bulk.hpp:203
Bulk::GetdNmax
OCP_DBL GetdNmax() const
Return dNmax.
Definition: Bulk.hpp:174
Bulk::GetdPmax
OCP_DBL GetdPmax() const
Return dPmax.
Definition: Bulk.hpp:172
Bulk::GetBulkNum
OCP_USI GetBulkNum() const
Return the number of bulks.
Definition: Bulk.hpp:145
Bulk::GetComNum
USI GetComNum() const
Return the number of components.
Definition: Bulk.hpp:147
Bulk::GetP
OCP_DBL GetP(const OCP_USI &n) const
Return pressure of the n-th bulk.
Definition: Bulk.hpp:155
DetailInfo
Collect more detailed information of each time step.
Definition: OCPOutput.hpp:165
Grid
Basic information of computational grid, including the rock properties.
Definition: Grid.hpp:76
OCP_AIMc
Definition: OCPFluidMethod.hpp:104
OCP_AIMc::Setup
void Setup(Reservoir &rs, LinearSystem &myLS, const OCPControl &ctrl)
Setup AIMc.
Definition: OCPFluidMethod.cpp:552
OCP_AIMt
perform AIM in time, that is, local FIM will be performed after global IMPEC performs
Definition: OCPFluidMethod.hpp:155
OCP_FIM
OCP_FIM is FIM (Fully Implicit Method).
Definition: OCPFluidMethod.hpp:52
OCP_FIM::resFIM
ResFIM resFIM
Resiual for FIM.
Definition: OCPFluidMethod.hpp:81
OCP_IMPEC
OCP_IMPEC is IMPEC (implict pressure explict saturation) method.
Definition: OCPFluidMethod.hpp:24
ParamEQUIL
Initial reservoir infomation for calculating initial equilibration.
Definition: Bulk.hpp:37
Well
Definition: Well.hpp:105
Generated on Sat Oct 15 2022 00:08:52 for OpenCAEPoro by doxygen 1.9.1