OpenCAEPoro  0.2.0 Sep/22/2022
A simulator for multicomponent porous media flow
Public Member Functions | Friends | List of all members
Bulk Class Reference

Physical information of each active reservoir bulk. More...

#include <Bulk.hpp>

Public Member Functions

void InputParam (ParamReservoir &rs_param)
 Input param from internal data structure ParamReservoir. More...
 
void Setup (const Grid &myGrid)
 Allocate memory for bulk data of grid. More...
 
void InitSjPcBo (const USI &tabrow)
 Calculate initial equilibrium for blkoil model according to EQUIL. More...
 
void InitSjPcComp (const USI &tabrow, const Grid &myGrid)
 Calculate initial equilibrium for compositional model according to EQUIL. More...
 
void InitFlash (const bool &flag=false)
 Perform flash calculation with saturations. More...
 
void InitFlashDer ()
 Perform flash calculation with saturations and calculate derivatives.
 
void InitFlashDer_n ()
 
void Flash ()
 Perform flash calculation with Ni. More...
 
void FlashBLKOIL ()
 Perform flash calculation with Ni in Black Oil Model.
 
void FlashCOMP ()
 Perform flash calculation with Ni in Compositional Model.
 
void FlashDeriv ()
 Perform flash calculation with Ni and calculate derivatives. More...
 
void FlashDeriv_n ()
 
void FlashDerivBLKOIL ()
 Perform flash calculation with Ni in Black Oil Model.
 
void FlashDerivBLKOIL_n ()
 
void FlashDerivCOMP ()
 Perform flash calculation with Ni in Compositional Model.
 
void FlashDerivCOMP_n ()
 
USI CalFlashType (const OCP_USI &n) const
 determine which flash type will be used
 
void PassFlashValue (const OCP_USI &n)
 Pass values from Flash to Bulk after Flash calculation.
 
void PassFlashValueAIMc (const OCP_USI &n)
 
void PassFlashValueDeriv (const OCP_USI &n)
 Pass derivative values from Flash to Bulk after Flash calculation.
 
void PassFlashValueDeriv_n (const OCP_USI &n)
 
void ResetFlash ()
 Reset variables in flash calculations.
 
void CalKrPc ()
 Calculate relative permeability and capillary pressure with saturation. More...
 
void CalKrPcDeriv ()
 Calculate relative permeability and capillary pressure and their derivatives.
 
void CalVpore ()
 Calculate volume of pore with pressure.
 
OCP_DBL CalFPR () const
 Calculate average pressure in reservoir.
 
void CalMaxChange ()
 Calculate max change of some variables.
 
void CheckSetup () const
 Check if error occurs in Setup.
 
void CheckInitVpore () const
 Check initial pore volume. More...
 
void CheckVpore () const
 Check pore volume.
 
bool CheckP () const
 Check if negative P occurs, return false if so. More...
 
bool CheckNi ()
 Check if negative Ni occurs, return false if so. More...
 
bool CheckVe (const OCP_DBL &Vlim) const
 Check if relative volume error is out of range, return false if so. More...
 
void CheckSat () const
 Check if the sum of saturations is one.
 
void CheckDiff ()
 Check difference from last time step, for Debug and Test.
 
OCP_USI GetBulkNum () const
 Return the number of bulks.
 
USI GetComNum () const
 Return the number of components.
 
USI GetPhaseNum () const
 Return the number of phases.
 
USI GetMixMode () const
 Return the mixture mode.
 
const vector< Mixture * > & GetMixture () const
 Return flash results (it has not been used by far).
 
OCP_DBL GetP (const OCP_USI &n) const
 Return pressure of the n-th bulk.
 
OCP_DBL GetSOIL (const OCP_USI &n) const
 Return oil saturation of the n-th bulk.
 
OCP_DBL GetSGAS (const OCP_USI &n) const
 Return gas saturation of the n-th bulk.
 
OCP_DBL GetSWAT (const OCP_USI &n) const
 Return water saturation of the n-th bulk.
 
OCP_DBL GetdPmax () const
 Return dPmax.
 
OCP_DBL GetdNmax () const
 Return dNmax.
 
OCP_DBL GetdSmax () const
 Return dSmax.
 
OCP_DBL GetdVmax () const
 Return dVmax.
 
void ResetphaseNum ()
 
void ResetminEigenSkip ()
 
void ResetflagSkip ()
 
void ResetziSkip ()
 
void ResetPSkip ()
 
void ResetKs ()
 
void ResetNt ()
 Reset Nt to the ones of the last time step.
 
void ResetP ()
 Reset P to the ones of the last time step.
 
void ResetPj ()
 Reset Pj to the ones of the last time step.
 
void ResetNi ()
 Reset Ni to the ones of the last time step.
 
void ResetKr ()
 Reset Kr to the ones of the last time step.
 
void ResetVp ()
 Reset Vp to the ones of the last time step.
 
void CalSomeInfo (const Grid &myGrid) const
 
void AllocateWellBulkId (const USI &n)
 Allocate memory for WellbulkId.
 
void ClearWellBulkId ()
 
OCP_DBL CalNT ()
 
void AllocateAuxIMPEC ()
 Allocate memory for auxiliary variables used for IMPEC.
 
void GetSolIMPEC (const vector< OCP_DBL > &u)
 Update P and Pj after linear system is solved.
 
void SetCFL2Zero () const
 Initialize the CFL number.
 
OCP_DBL CalCFL () const
 Calculate the CFL number.
 
void UpdateLastStepIMPEC ()
 Update value of last step for IMPEC.
 
void AllocateAuxFIM ()
 Allocate memory for auxiliary variables used for FIM.
 
void GetSolFIM (const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
 Get the solution for FIM after a Newton iteration.
 
void GetSolFIM_n (const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
 
void GetSol01FIM (const vector< OCP_DBL > &u)
 Get the solution for FIM after a Newton iteration???
 
void CalRelResFIM (ResFIM &resFIM) const
 Calculate relative resiual for FIM.
 
void ShowRes (const vector< OCP_DBL > &res) const
 
void ResetFIM ()
 Reset FIM.
 
void UpdateLastStepFIM ()
 Update values of last step for FIM.
 
OCP_DBL CalNRdSmax (OCP_USI &index)
 Calculate some auxiliary variable, for example, dSmax.
 
OCP_DBL GetNRdPmax ()
 Return NRdPmax.
 
OCP_DBL GetNRdSmaxP ()
 Return NRdSmaxP.
 
OCP_DBL GetNRdNmax ()
 
void CorrectNi (const vector< OCP_DBL > &res)
 
void OutputInfo (const OCP_USI &n) const
 
OCP_ULL GetSSMSTAiters () const
 
OCP_ULL GetNRSTAiters () const
 
OCP_ULL GetSSMSPiters () const
 
OCP_ULL GetNRSPiters () const
 
OCP_ULL GetRRiters () const
 
OCP_ULL GetSSMSTAcounts () const
 
OCP_ULL GetNRSTAcounts () const
 
OCP_ULL GetSSMSPcounts () const
 
OCP_ULL GetNRSPcounts () const
 
OCP_ULL GetRRcounts () const
 
void AllocateAuxAIM (const OCP_DBL &ratio)
 Allocate memory for auxiliary variables used for AIMt.
 
OCP_USI GetMaxFIMBulk () const
 
void FlashDerivAIM (const bool &IfAIMs)
 Perform flash calculation with Ni and calculate derivatives.
 
void PassFlashValueDerivAIM (const OCP_USI &n)
 
void CalKrPcDerivAIM (const bool &IfAIMs)
 Calculate relative permeability and capillary pressure and their derivatives.
 
void CalRelResAIMt (ResFIM &resFIM) const
 Calculate relative resiual for local FIM.
 
void GetSolAIMt (const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
 Get the solution for local FIM after a Newton iteration.
 
void CalRelResAIMs (ResFIM &resFIM) const
 Calculate relative resiual for AIMs, parts related to FIM are considered.
 
void GetSolAIMs (const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
 
void UpdateLastStepAIM ()
 
void ResetFIMBulk ()
 
void ShowFIMBulk (const bool &flag=false) const
 
bool CheckNiFIMBulk () const
 Check if negative Ni occurs, return false if so.
 
void InFIMNi ()
 Ni in FIM Bulk -> FIMNi.
 
void OutFIMNi ()
 FIMNi -> Ni in FIM Bulk.
 
void AllocateAuxAIMc ()
 
void FlashAIMc ()
 Perform flash calculation with Ni.
 
void FlashBLKOILAIMc ()
 Perform flash calculation with Ni in Black Oil Model.
 
void FlashCOMPAIMc ()
 Perform flash calculation with Ni in Compositional Model.
 
void FlashAIMc01 ()
 
void FlashBLKOILAIMc01 ()
 
void FlashCOMPAIMc01 ()
 
void FlashDerivAIMc ()
 Perform flash calculation with Ni and calculate derivatives.
 
void FlashDerivBLKOILAIMc ()
 Perform flash calculation with Ni in Black Oil Model.
 
void FlashDerivCOMPAIMc ()
 Perform flash calculation with Ni in Compositional Model.
 
void CalKrPcAIMc ()
 Calculate relative permeability and capillary pressure with saturation.
 
void CalKrPcDerivAIMc ()
 Calculate relative permeability and capillary pressure and their derivatives.
 
void GetSolAIMc (const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
 
void GetSolAIMc01 (const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
 
void UpdatePj ()
 

Friends

class BulkConn
 
class Well
 
class DetailInfo
 
class OCP_IMPEC
 
class OCP_FIM
 
class Reservoir
 
class OCP_AIMt
 
class OCP_AIMc
 

Detailed Description

Physical information of each active reservoir bulk.

Definition at line 57 of file Bulk.hpp.

Member Function Documentation

◆ CalKrPc()

void Bulk::CalKrPc ( )

Calculate relative permeability and capillary pressure with saturation.

Relative permeability and capillary pressure.

Definition at line 1897 of file Bulk.cpp.

1898 {
1899  OCP_FUNCNAME;
1900 
1901  if (!miscible) {
1902  OCP_DBL tmp = 0;
1903  for (OCP_USI n = 0; n < numBulk; n++) {
1904  OCP_USI bId = n * numPhase;
1905  flow[SATNUM[n]]->CalKrPc(&S[bId], &kr[bId], &Pc[bId], 0, tmp, tmp);
1906  for (USI j = 0; j < numPhase; j++)
1907  Pj[n * numPhase + j] = P[n] + Pc[n * numPhase + j];
1908  }
1909  } else {
1910  for (OCP_USI n = 0; n < numBulk; n++) {
1911  OCP_USI bId = n * numPhase;
1912  flow[SATNUM[n]]->CalKrPc(&S[bId], &kr[bId], &Pc[bId], surTen[n], Fk[n],
1913  Fp[n]);
1914  for (USI j = 0; j < numPhase; j++)
1915  Pj[n * numPhase + j] = P[n] + Pc[n * numPhase + j];
1916  }
1917  }
1918 
1919  if (ScalePcow) {
1920  // correct
1921  USI Wid = phase2Index[WATER];
1922  for (USI n = 0; n < numBulk; n++) {
1923  Pc[n * numPhase + Wid] *= ScaleValuePcow[n];
1924  Pj[n * numPhase + Wid] = P[n] + Pc[n * numPhase + Wid];
1925  }
1926  }
1927 }
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
WATER
const USI WATER
Fluid type = water.
Definition: OCPConst.hpp:84
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:24
OCP_FUNCNAME
#define OCP_FUNCNAME
Print Function Name.
Definition: UtilError.hpp:73

References OCP_FUNCNAME, and WATER.

◆ CheckInitVpore()

void Bulk::CheckInitVpore ( ) const

Check initial pore volume.

Check init pore volume.

Definition at line 349 of file Bulk.cpp.

350 {
351  for (OCP_USI n = 0; n < numBulk; n++) {
352  if (rockVpInit[n] < TINY) {
353  OCP_ABORT("Bulk volume is too small: bulk[" + std::to_string(n) +
354  "] = " + std::to_string(rockVpInit[n]));
355  }
356  }
357 }
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:36
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47

References OCP_ABORT, and TINY.

◆ CheckNi()

bool Bulk::CheckNi ( )

Check if negative Ni occurs, return false if so.

Return true if no negative Ni and false otherwise.

Definition at line 2114 of file Bulk.cpp.

2115 {
2116  OCP_FUNCNAME;
2117 
2118  OCP_USI len = numBulk * numCom;
2119  for (OCP_USI n = 0; n < len; n++) {
2120  if (Ni[n] < 0) {
2121  OCP_USI bId = n / numCom;
2122  if (Ni[n] > -1E-3 * Nt[bId] && false) {
2123  Ni[n] = 1E-8 * Nt[bId];
2124  // Ni[n] = 0;
2125  } else {
2126  USI cId = n - bId * numCom;
2127  std::ostringstream NiStringSci;
2128  NiStringSci << std::scientific << Ni[n];
2129  OCP_WARNING("Negative Ni: Ni[" + std::to_string(cId) + "] in Bulk[" +
2130  std::to_string(bId) + "] = " + NiStringSci.str() + " " +
2131  "dNi = " + std::to_string(dNNR[n]) + +" lNi[" +
2132  std::to_string(cId) + "] in Bulk[" + std::to_string(bId) +
2133  "] = " + std::to_string(lNi[n]) +
2134  " Nt = " + std::to_string(Nt[bId]) + " " +
2135  std::to_string(phaseNum[bId]) + " " +
2136  std::to_string(NRstep[bId]));
2137  for (USI i = 0; i < numCom; i++) {
2138  cout << Ni[bId * numCom + i] << " ";
2139  }
2140  cout << endl;
2141  for (USI j = 0; j < numPhase - 1; j++) {
2142  cout << setprecision(9);
2143  if (phaseExist[bId * numPhase + j]) {
2144  cout << j << " ";
2145  for (USI i = 0; i < numCom_1; i++) {
2146  cout << xij[bId * numPhase * numCom + j * numCom + i]
2147  << " ";
2148  }
2149  cout << endl;
2150  }
2151  }
2152  return false;
2153  }
2154  }
2155  }
2156  return true;
2157 }
OCP_WARNING
#define OCP_WARNING(msg)
Log warning messages.
Definition: UtilError.hpp:39

References OCP_FUNCNAME, and OCP_WARNING.

◆ CheckP()

bool Bulk::CheckP ( ) const

Check if negative P occurs, return false if so.

Return true if no negative pressure and false otherwise.

Definition at line 2095 of file Bulk.cpp.

2096 {
2097  OCP_FUNCNAME;
2098 
2099  for (OCP_USI n = 0; n < numBulk; n++) {
2100  if (P[n] < 0) {
2101  std::ostringstream PStringSci;
2102  PStringSci << std::scientific << P[n];
2103  OCP_WARNING("Negative pressure: P[" + std::to_string(n) +
2104  "] = " + PStringSci.str());
2105  cout << "P = " << P[n] << endl;
2106  return false;
2107  }
2108  }
2109 
2110  return true;
2111 }

References OCP_FUNCNAME, and OCP_WARNING.

◆ CheckVe()

bool Bulk::CheckVe ( const OCP_DBL &  Vlim) const

Check if relative volume error is out of range, return false if so.

Return true if all Ve < Vlim and false otherwise.

Definition at line 2160 of file Bulk.cpp.

2161 {
2162  OCP_FUNCNAME;
2163 
2164  // bool tmpflag = true;
2165 
2166  OCP_DBL dVe = 0.0;
2167  for (OCP_USI n = 0; n < numBulk; n++) {
2168  dVe = fabs(vf[n] - rockVp[n]) / rockVp[n];
2169  if (dVe > Vlim) {
2170  cout << "Volume error at Bulk[" << n << "] = " << setprecision(6) << dVe
2171  << " is too big!" << endl;
2172  // OutputInfo(n);
2173  return false;
2174  // tmpflag = false;
2175  }
2176  }
2177  // OutputInfo(39);
2178  // if (!tmpflag) return false;
2179  return true;
2180 }

References OCP_FUNCNAME.

◆ Flash()

void Bulk::Flash ( )

Perform flash calculation with Ni.

Use moles of component and pressure both in blackoil and compositional model.

Definition at line 1340 of file Bulk.cpp.

1341 {
1342  OCP_FUNCNAME;
1343 
1344  if (comps) {
1345  FlashCOMP();
1346  } else {
1347  FlashBLKOIL();
1348  }
1349 
1350 #ifdef DEBUG
1351  CheckSat();
1352 #endif // DEBUG
1353 }
Bulk::FlashBLKOIL
void FlashBLKOIL()
Perform flash calculation with Ni in Black Oil Model.
Definition: Bulk.cpp:1355
Bulk::CheckSat
void CheckSat() const
Check if the sum of saturations is one.
Definition: Bulk.cpp:2249
Bulk::FlashCOMP
void FlashCOMP()
Perform flash calculation with Ni in Compositional Model.
Definition: Bulk.cpp:1363

References CheckSat(), FlashBLKOIL(), FlashCOMP(), and OCP_FUNCNAME.

◆ FlashDeriv()

void Bulk::FlashDeriv ( )

Perform flash calculation with Ni and calculate derivatives.

Use moles of component and pressure both in blackoil and compositional model.

Definition at line 1401 of file Bulk.cpp.

1402 {
1403  OCP_FUNCNAME;
1404 
1405  if (comps) {
1406  FlashDerivCOMP();
1407  } else {
1408  FlashDerivBLKOIL();
1409  }
1410 
1411 #ifdef DEBUG
1412  CheckSat();
1413 #endif // DEBUG
1414 }
Bulk::FlashDerivBLKOIL
void FlashDerivBLKOIL()
Perform flash calculation with Ni in Black Oil Model.
Definition: Bulk.cpp:1428
Bulk::FlashDerivCOMP
void FlashDerivCOMP()
Perform flash calculation with Ni in Compositional Model.
Definition: Bulk.cpp:1440

References CheckSat(), FlashDerivBLKOIL(), FlashDerivCOMP(), and OCP_FUNCNAME.

◆ GetSolFIM_n()

void Bulk::GetSolFIM_n ( const vector< OCP_DBL > &  u,
const OCP_DBL &  dPmaxlim,
const OCP_DBL &  dSmaxlim 
)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Definition at line 2704 of file Bulk.cpp.

2706 {
2707  // For saturations changes:
2708  // 1. maximum changes must be less than dSmaxlim,
2709  // 2. if phases become mobile/immobile, then set it to crtical point,
2710  // 3. After saturations are determined, then scale the nij to conserve Volume
2711  // equations.
2712 
2713  OCP_FUNCNAME;
2714 
2715  dSNR = S;
2716  NRphaseNum = phaseNum;
2717  fill(dSNRP.begin(), dSNRP.end(), 0.0);
2718 
2719  NRdSmaxP = 0;
2720  NRdPmax = 0;
2721  NRdNmax = 0;
2722 
2723  USI row = numPhase * (numCom + 1);
2724  USI ncol = numCom + 1;
2725  USI n_np_j;
2726 
2727  vector<OCP_DBL> dtmp(row, 0);
2728  vector<OCP_DBL> tmpNij(numPhase * numCom, 0);
2729 
2730  OCP_DBL dSmax;
2731  OCP_DBL dP;
2732  OCP_DBL chop = 1;
2733 
2734  // vector<OCP_DBL> tmpxi(numPhase, 0);
2735  // vector<OCP_DBL> tmpVf(numPhase, 0);
2736 
2737  for (OCP_USI n = 0; n < numBulk; n++) {
2738  const vector<OCP_DBL>& scm = satcm[SATNUM[n]];
2739 
2740  dP = u[n * ncol];
2741  dPNR[n] = dP;
2742  P[n] += dP; // seems better
2743  if (fabs(NRdPmax) < fabs(dP)) NRdPmax = dP;
2744 
2745  // rockVp[n] = rockVpInit[n] * (1 + rockC1 * (P[n] - rockPref));
2746  // cout << scientific << setprecision(6) << dP << " " << n << endl;
2747 
2748  dSmax = 0;
2749  chop = 1;
2750 
2751  const USI cNp = phaseNum[n] + 1;
2752  const USI len = resIndex[n + 1] - resIndex[n];
2754  /*Dcopy(cNp, &dtmp[0], &res_n[resIndex[n]]);
2755  DaAxpby(cNp, ncol, 1.0, dSec_dPri.data() + dSdPindex[n], u.data() + n * ncol,
2756  1.0, dtmp.data());*/
2757  Dcopy(len, &dtmp[0], &res_n[resIndex[n]]);
2758  DaAxpby(len, ncol, 1.0, dSec_dPri.data() + dSdPindex[n], u.data() + n * ncol,
2759  1.0, dtmp.data());
2760 
2761  // Calculate dSmax
2762  USI js = 0;
2763  for (USI j = 0; j < numPhase; j++) {
2764  n_np_j = n * numPhase + j;
2765  if (phaseExist[n_np_j]) {
2766  dSmax = max(fabs(dtmp[js]), dSmax);
2767  js++;
2768  }
2769  }
2770 
2771  // Calculate chop
2772  if (dSmax > dSmaxlim) {
2773  chop = min(chop, dSmaxlim / dSmax);
2774  dSmax = dSmaxlim;
2775  }
2776  js = 0;
2777  for (USI j = 0; j < numPhase; j++) {
2778  n_np_j = n * numPhase + j;
2779  if (phaseExist[n_np_j]) {
2780  if (fabs(S[n_np_j] - scm[j]) > TINY &&
2781  (S[n_np_j] - scm[j]) / (chop * dtmp[js]) < 0)
2782  chop *= min(1.0, -((S[n_np_j] - scm[j]) / (chop * dtmp[js])));
2783  if (S[n_np_j] + chop * dtmp[js] < 0)
2784  chop *= min(1.0, -S[n_np_j] / (chop * dtmp[js]));
2785  js++;
2786  }
2787  }
2788 
2789  // Calculate S, phaseExist, xij, nj
2790  fill(tmpNij.begin(), tmpNij.end(), 0.0);
2791  // fill(Ni.begin() + n * numCom, Ni.begin() + n * numCom + numCom, 0.0);
2792 
2793  js = 0;
2794  USI jx = cNp;
2795  for (USI j = 0; j < numPhase; j++) {
2796  n_np_j = n * numPhase + j;
2797  if (phaseExist[n_np_j]) {
2798  dSNRP[n_np_j] = chop * dtmp[js];
2799  if (fabs(NRdSmaxP) < fabs(dSNRP[n_np_j])) NRdSmaxP = dSNRP[n_np_j];
2800  js++;
2801 
2802  // S[n_np_j] += chop * dtmp[js];
2803  // if (S[n_np_j] < TINY) {
2804  // phaseExist[n_np_j] = false;
2805  // }
2806  // js++;
2807  Daxpy(numCom, nj[n_np_j], &xij[n_np_j * numCom], &tmpNij[j * numCom]);
2808  for (USI i = 0; i < pEnumCom[j]; i++) {
2809  tmpNij[j * numCom + i] += chop * dtmp[jx + i];
2810  }
2811  jx += pEnumCom[j];
2812  nj[n_np_j] = Dnorm1(numCom, &tmpNij[j * numCom]);
2813  for (USI i = 0; i < numCom; i++) {
2814  xij[n_np_j * numCom + i] = tmpNij[j * numCom + i] / nj[n_np_j];
2815  // Ni[n * numCom + i] += tmpNij[j * numCom + i];
2816  }
2817  }
2818  }
2819  if (phaseNum[n] == 2 && false) {
2820  for (USI i = 0; i < numCom_1; i++) {
2821  Ks[n * numCom_1 + i] = xij[n * numPhase * numCom + i] / xij[n * numPhase * numCom + numCom + i];
2822  }
2823  }
2824 
2825 
2826  // Vf correction
2827  /* OCP_DBL dVmin = 100 * rockVp[n];
2828  USI index = 0;
2829  for (USI j = 0; j < numPhase; j++) {
2830  if (phaseExist[n * numPhase + j]) {
2831  tmpxi[j] = flashCal[PVTNUM[n]]->XiPhase(P[n], T, &xij[(n *
2832  numPhase + j) * numCom]); tmpVf[j] = nj[n * numPhase + j] / (S[n * numPhase +
2833  j] * tmpxi[j]); if (fabs(tmpVf[j] - rockVp[n]) < dVmin) { dVmin =
2834  fabs(tmpVf[j] - rockVp[n]); index = j;
2835  }
2836  }
2837  } */
2838  // for (USI j = 0; j < numPhase; j++) {
2839  // if (phaseExist[n * numPhase + j]) {
2840  // nj[n * numPhase + j] *= (tmpVf[index] / tmpVf[j]);
2841  // for (USI i = 0; i < numCom; i++) {
2842  // Ni[n * numCom + i] += nj[n * numPhase + j] * xij[(n * numPhase +
2843  // j) * numCom + i];
2844  // }
2845  // }
2846  // }
2847 
2848  NRstep[n] = chop;
2849  for (USI i = 0; i < numCom; i++) {
2850  dNNR[n * numCom + i] = u[n * ncol + 1 + i] * chop;
2851  if (fabs(NRdNmax) < fabs(dNNR[n * numCom + i]) / Nt[n])
2852  NRdNmax = dNNR[n * numCom + i] / Nt[n];
2853  Ni[n * numCom + i] += dNNR[n * numCom + i];
2854  }
2855  }
2856 }
DaAxpby
void DaAxpby(const int &m, const int &n, const double &a, const double *A, const double *x, const double &b, double *y)
Computes y = a A x + b y.
Definition: DenseMat.cpp:173
Daxpy
void Daxpy(const int &n, const double &alpha, const double *x, double *y)
Constant times a vector plus a vector.
Definition: DenseMat.cpp:69
Dnorm1
double Dnorm1(const int &N, double *x)
Computes the L1-norm of a vector.
Definition: DenseMat.cpp:50
Dcopy
void Dcopy(const int &N, double *dst, const double *src)
Calculate the minimal eigenvalue for sysmetric matrix with mkl lapack.
Definition: DenseMat.cpp:37

References DaAxpby(), Daxpy(), Dcopy(), Dnorm1(), OCP_FUNCNAME, and TINY.

◆ InitFlash()

void Bulk::InitFlash ( const bool &  flag = false)

Perform flash calculation with saturations.

Use initial saturation in blackoil model and initial Zi in compositional model. It gives initial properties and some derivatives for IMPEC It just gives Ni for FIM

Definition at line 1280 of file Bulk.cpp.

1281 {
1282  OCP_FUNCNAME;
1283 
1284  for (OCP_USI n = 0; n < numBulk; n++) {
1285  flashCal[PVTNUM[n]]->InitFlash(P[n], Pb[n], T, &S[n * numPhase], rockVp[n],
1286  initZi.data() + n * numCom);
1287  for (USI i = 0; i < numCom; i++) {
1288  Ni[n * numCom + i] = flashCal[PVTNUM[n]]->Ni[i];
1289  }
1290  if (flag) {
1291  PassFlashValue(n);
1292  }
1293  }
1294 
1295 #ifdef DEBUG
1296  if (flag) {
1297  CheckSat();
1298  }
1299 #endif // DEBUG
1300 }
Bulk::PassFlashValue
void PassFlashValue(const OCP_USI &n)
Pass values from Flash to Bulk after Flash calculation.
Definition: Bulk.cpp:1529

References CheckSat(), OCP_FUNCNAME, and PassFlashValue().

◆ InitSjPcBo()

void Bulk::InitSjPcBo ( const USI &  tabrow)

Calculate initial equilibrium for blkoil model according to EQUIL.

Here tabrow is maximum number of depth nodes in table of depth vs pressure.

Definition at line 371 of file Bulk.cpp.

372 {
373  OCP_FUNCNAME;
374 
375  OCP_DBL Dref = EQUIL.Dref;
376  OCP_DBL Pref = EQUIL.Pref;
377  OCP_DBL DOWC = EQUIL.DOWC;
378  OCP_DBL PcOW = EQUIL.PcOW;
379  OCP_DBL DOGC = EQUIL.DGOC;
380  OCP_DBL PcGO = EQUIL.PcGO;
381  OCP_DBL Zmin = 1E8;
382  OCP_DBL Zmax = 0;
383 
384  for (OCP_USI n = 0; n < numBulk; n++) {
385  OCP_DBL temp1 = depth[n] - dz[n] / 2;
386  OCP_DBL temp2 = depth[n] + dz[n] / 2;
387  Zmin = Zmin < temp1 ? Zmin : temp1;
388  Zmax = Zmax > temp2 ? Zmax : temp2;
389  }
390  OCP_DBL tabdz = (Zmax - Zmin) / (tabrow - 1);
391 
392  // create table
393  OCPTable DepthP(tabrow, 4);
394  vector<OCP_DBL>& Ztmp = DepthP.GetCol(0);
395  vector<OCP_DBL>& Potmp = DepthP.GetCol(1);
396  vector<OCP_DBL>& Pgtmp = DepthP.GetCol(2);
397  vector<OCP_DBL>& Pwtmp = DepthP.GetCol(3);
398 
399  // cal Tab_Ztmp
400  Ztmp[0] = Zmin;
401  for (USI i = 1; i < tabrow; i++) {
402  Ztmp[i] = Ztmp[i - 1] + tabdz;
403  }
404 
405  // find the RefId
406  USI beginId = 0;
407  if (Dref <= Ztmp[0]) {
408  beginId = 0;
409  } else if (Dref >= Ztmp[tabrow - 1]) {
410  beginId = tabrow - 1;
411  } else {
412  beginId =
413  distance(Ztmp.begin(), find_if(Ztmp.begin(), Ztmp.end(),
414  [s = Dref](auto& t) { return t > s; }));
415  beginId--;
416  }
417 
418  // begin calculating oil pressure
419  OCP_DBL Pbb = Pref;
420  OCP_DBL gammaOtmp, gammaWtmp, gammaGtmp;
421  OCP_DBL Ptmp;
422  USI mynum = 10;
423  OCP_DBL mydz = 0;
424  OCP_DBL Poref, Pgref, Pwref;
425  OCP_DBL Pbegin = 0;
426 
427  if (Dref < DOGC) {
428 
429  // reference pressure is gas pressure
430  if (!gas) OCP_ABORT("SGOF is missing!");
431 
432  Pgref = Pref;
433  gammaGtmp = flashCal[0]->GammaPhaseG(Pgref);
434  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
435  Pgtmp[beginId] = Pbegin;
436 
437  // find the gas pressure
438  for (USI id = beginId; id > 0; id--) {
439  gammaGtmp = flashCal[0]->GammaPhaseG(Pgtmp[id]);
440  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
441  }
442 
443  for (USI id = beginId; id < tabrow - 1; id++) {
444  gammaGtmp = flashCal[0]->GammaPhaseG(Pgtmp[id]);
445  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
446  }
447 
448  // find the oil pressure in Dref by Pgref
449  Poref = 0;
450  Ptmp = Pgref;
451  mydz = (DOGC - Dref) / mynum;
452 
453  for (USI i = 0; i < mynum; i++) {
454  gammaGtmp = flashCal[0]->GammaPhaseG(Ptmp);
455  Ptmp += gammaGtmp * mydz;
456  }
457  Ptmp -= PcGO;
458  for (USI i = 0; i < mynum; i++) {
459  if (!EQUIL.PBVD.IsEmpty()) {
460  Pbb = EQUIL.PBVD.Eval(0, DOGC - i * mydz, 1);
461  }
462  gammaOtmp = flashCal[0]->GammaPhaseO(Ptmp, Pbb);
463  Ptmp -= gammaOtmp * mydz;
464  }
465  Poref = Ptmp;
466 
467  // find the oil pressure in tab
468  if (!EQUIL.PBVD.IsEmpty()) {
469  Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
470  }
471  gammaOtmp = flashCal[0]->GammaPhaseO(Poref, Pbb);
472  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
473  Potmp[beginId] = Pbegin;
474 
475  for (USI id = beginId; id > 0; id--) {
476  if (!EQUIL.PBVD.IsEmpty()) {
477  Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
478  }
479  gammaOtmp = flashCal[0]->GammaPhaseO(Potmp[id], Pbb);
480  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
481  }
482 
483  for (USI id = beginId; id < tabrow - 1; id++) {
484  if (!EQUIL.PBVD.IsEmpty()) {
485  Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
486  }
487  gammaOtmp = flashCal[0]->GammaPhaseO(Potmp[id], Pbb);
488  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
489  }
490 
491  // find the water pressure in Dref by Poref
492  Pwref = 0;
493  Ptmp = Poref;
494  mydz = (DOWC - Dref) / mynum;
495 
496  for (USI i = 0; i < mynum; i++) {
497  if (!EQUIL.PBVD.IsEmpty()) {
498  Pbb = EQUIL.PBVD.Eval(0, Dref + i * mydz, 1);
499  }
500  gammaOtmp = flashCal[0]->GammaPhaseO(Ptmp, Pbb);
501  Ptmp += gammaOtmp * mydz;
502  }
503  Ptmp -= PcOW;
504  for (USI i = 0; i < mynum; i++) {
505  gammaWtmp = flashCal[0]->GammaPhaseW(Ptmp);
506  Ptmp -= gammaWtmp * mydz;
507  }
508  Pwref = Ptmp;
509 
510  // find the water pressure in tab
511  gammaWtmp = flashCal[0]->GammaPhaseW(Pwref);
512  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
513  Pwtmp[beginId] = Pbegin;
514 
515  for (USI id = beginId; id > 0; id--) {
516  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
517  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
518  }
519 
520  for (USI id = beginId; id < tabrow - 1; id++) {
521  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
522  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
523  }
524  } else if (Dref > DOWC) {
525 
526  // reference pressure is water pressure
527  Pwref = Pref;
528  gammaWtmp = flashCal[0]->GammaPhaseW(Pwref);
529  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
530  Pwtmp[beginId] = Pbegin;
531 
532  // find the water pressure
533  for (USI id = beginId; id > 0; id--) {
534  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
535  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
536  }
537  for (USI id = beginId; id < tabrow - 1; id++) {
538  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
539  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
540  }
541 
542  // find the oil pressure in Dref by Pwref
543  Poref = 0;
544  Ptmp = Pwref;
545  mydz = (DOWC - Dref) / mynum;
546 
547  for (USI i = 0; i < mynum; i++) {
548  gammaWtmp = flashCal[0]->GammaPhaseW(Ptmp);
549  Ptmp += gammaWtmp * mydz;
550  }
551  Ptmp += PcOW;
552 
553  for (USI i = 0; i < mynum; i++) {
554  if (!EQUIL.PBVD.IsEmpty()) {
555  Pbb = EQUIL.PBVD.Eval(0, DOWC - i * mydz, 1);
556  }
557  gammaOtmp = flashCal[0]->GammaPhaseO(Ptmp, Pbb);
558  Ptmp -= gammaOtmp * mydz;
559  }
560  Poref = Ptmp;
561 
562  // find the oil pressure in tab
563  if (!EQUIL.PBVD.IsEmpty()) {
564  Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
565  }
566  gammaOtmp = flashCal[0]->GammaPhaseO(Poref, Pbb);
567  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
568  Potmp[beginId] = Pbegin;
569 
570  for (USI id = beginId; id > 0; id--) {
571  if (!EQUIL.PBVD.IsEmpty()) {
572  Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
573  }
574  gammaOtmp = flashCal[0]->GammaPhaseO(Potmp[id], Pbb);
575  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
576  }
577 
578  for (USI id = beginId; id < tabrow - 1; id++) {
579  if (!EQUIL.PBVD.IsEmpty()) {
580  Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
581  }
582  gammaOtmp = flashCal[0]->GammaPhaseO(Potmp[id], Pbb);
583  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
584  }
585 
586  if (gas) {
587  // find the gas pressure in Dref by Poref
588  Pgref = 0;
589  Ptmp = Poref;
590  mydz = (DOGC - Dref) / mynum;
591 
592  for (USI i = 0; i < mynum; i++) {
593  if (!EQUIL.PBVD.IsEmpty()) {
594  Pbb = EQUIL.PBVD.Eval(0, Dref + i * mydz, 1);
595  }
596  gammaOtmp = flashCal[0]->GammaPhaseO(Ptmp, Pbb);
597  Ptmp += gammaOtmp * mydz;
598  }
599  Ptmp += PcGO;
600  for (USI i = 0; i < mynum; i++) {
601  gammaGtmp = flashCal[0]->GammaPhaseG(Ptmp);
602  Ptmp -= gammaGtmp * mydz;
603  }
604  Pgref = Ptmp;
605 
606  // find the gas pressure in tab
607  gammaGtmp = flashCal[0]->GammaPhaseG(Pgref);
608  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
609  Pgtmp[beginId] = Pbegin;
610 
611  for (USI id = beginId; id > 0; id--) {
612  gammaGtmp = flashCal[0]->GammaPhaseG(Pgtmp[id]);
613  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
614  }
615  for (USI id = beginId; id < tabrow - 1; id++) {
616  gammaGtmp = flashCal[0]->GammaPhaseG(Pgtmp[id]);
617  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
618  }
619  }
620  } else {
621 
622  // reference pressure is oil pressure
623  Poref = Pref;
624  if (!EQUIL.PBVD.IsEmpty()) {
625  Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
626  }
627  gammaOtmp = flashCal[0]->GammaPhaseO(Poref, Pbb);
628  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
629  Potmp[beginId] = Pbegin;
630 
631  // find the oil pressure in tab
632  for (USI id = beginId; id > 0; id--) {
633  if (!EQUIL.PBVD.IsEmpty()) {
634  Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
635  }
636  gammaOtmp = flashCal[0]->GammaPhaseO(Potmp[id], Pbb);
637  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
638  }
639  for (USI id = beginId; id < tabrow - 1; id++) {
640  if (!EQUIL.PBVD.IsEmpty()) {
641  Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
642  }
643  gammaOtmp = flashCal[0]->GammaPhaseO(Potmp[id], Pbb);
644  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
645  }
646 
647  if (gas) {
648  // find the gas pressure in Dref by Poref
649  Pgref = 0;
650  Ptmp = Poref;
651  mydz = (DOGC - Dref) / mynum;
652 
653  for (USI i = 0; i < mynum; i++) {
654  if (!EQUIL.PBVD.IsEmpty()) {
655  Pbb = EQUIL.PBVD.Eval(0, Dref + i * mydz, 1);
656  }
657  gammaOtmp = flashCal[0]->GammaPhaseO(Ptmp, Pbb);
658  Ptmp += gammaOtmp * mydz;
659  }
660  Ptmp += PcGO;
661  for (USI i = 0; i < mynum; i++) {
662  gammaGtmp = flashCal[0]->GammaPhaseG(Ptmp);
663  Ptmp -= gammaGtmp * mydz;
664  }
665  Pgref = Ptmp;
666 
667  // find the gas pressure in tab
668  gammaGtmp = flashCal[0]->GammaPhaseG(Pgref);
669  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
670  Pgtmp[beginId] = Pbegin;
671 
672  for (USI id = beginId; id > 0; id--) {
673  gammaGtmp = flashCal[0]->GammaPhaseG(Pgtmp[id]);
674  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
675  }
676 
677  for (USI id = beginId; id < tabrow - 1; id++) {
678  gammaGtmp = flashCal[0]->GammaPhaseG(Pgtmp[id]);
679  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
680  }
681  }
682 
683  // find the water pressure in Dref by Poref
684  Pwref = 0;
685  Ptmp = Poref;
686  mydz = (DOWC - Dref) / mynum;
687 
688  for (USI i = 0; i < mynum; i++) {
689  if (!EQUIL.PBVD.IsEmpty()) {
690  Pbb = EQUIL.PBVD.Eval(0, Dref + i * mydz, 1);
691  }
692  gammaOtmp = flashCal[0]->GammaPhaseO(Ptmp, Pbb);
693  Ptmp += gammaOtmp * mydz;
694  }
695  Ptmp -= PcOW;
696  for (USI i = 0; i < mynum; i++) {
697  gammaWtmp = flashCal[0]->GammaPhaseW(Ptmp);
698  Ptmp -= gammaWtmp * mydz;
699  }
700  Pwref = Ptmp;
701 
702  // find the water pressure in tab
703  gammaWtmp = flashCal[0]->GammaPhaseW(Pwref);
704  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
705  Pwtmp[beginId] = Pbegin;
706 
707  for (USI id = beginId; id > 0; id--) {
708  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
709  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
710  }
711 
712  for (USI id = beginId; id < tabrow - 1; id++) {
713  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
714  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
715  }
716  }
717 
718  DepthP.Display();
719 
720  // calculate Pc from DepthP to calculate Sj
721  std::vector<OCP_DBL> data(4, 0), cdata(4, 0);
722  // if capillary between water and oil is considered
723  vector<bool> FlagPcow(NTSFUN, true);
724  for (USI i = 0; i < NTSFUN; i++) {
725  if (fabs(flow[i]->GetPcowBySw(0.0 - TINY)) < TINY &&
726  fabs(flow[i]->GetPcowBySw(1.0 + TINY) < TINY)) {
727  FlagPcow[i] = false;
728  }
729  }
730 
731  for (OCP_USI n = 0; n < numBulk; n++) {
732  DepthP.Eval_All(0, depth[n], data, cdata);
733  OCP_DBL Po = data[1];
734  OCP_DBL Pg = data[2];
735  OCP_DBL Pw = data[3];
736  OCP_DBL Pcgo = Pg - Po;
737  OCP_DBL Pcow = Po - Pw;
738  OCP_DBL Sw = flow[SATNUM[n]]->GetSwByPcow(Pcow);
739  OCP_DBL Sg = 0;
740  if (gas) {
741  Sg = flow[SATNUM[n]]->GetSgByPcgo(Pcgo);
742  }
743  if (Sw + Sg > 1) {
744  // should me modified
745  OCP_DBL Pcgw = Pcow + Pcgo;
746  Sw = flow[SATNUM[n]]->GetSwByPcgw(Pcgw);
747  Sg = 1 - Sw;
748  }
749 
750  if (1 - Sw < TINY) {
751  // all water
752  Po = Pw + flow[SATNUM[n]]->GetPcowBySw(1.0);
753  } else if (1 - Sg < TINY) {
754  // all gas
755  Po = Pg - flow[SATNUM[n]]->GetPcgoBySg(1.0);
756  } else if (1 - Sw - Sg < TINY) {
757  // water and gas
758  Po = Pg - flow[SATNUM[n]]->GetPcgoBySg(Sg);
759  }
760  P[n] = Po;
761 
762  if (depth[n] < DOGC) {
763  Pbb = Po;
764  } else if (!EQUIL.PBVD.IsEmpty()) {
765  Pbb = EQUIL.PBVD.Eval(0, depth[n], 1);
766  }
767  Pb[n] = Pbb;
768 
769  // cal Sg and Sw
770  Sw = 0;
771  Sg = 0;
772  USI ncut = 10;
773 
774  for (USI k = 0; k < ncut; k++) {
775  OCP_DBL tmpSw = 0;
776  OCP_DBL tmpSg = 0;
777  OCP_DBL dep = depth[n] + dz[n] / ncut * (k - (ncut - 1) / 2.0);
778  DepthP.Eval_All(0, dep, data, cdata);
779  Po = data[1];
780  Pg = data[2];
781  Pw = data[3];
782  Pcow = Po - Pw;
783  Pcgo = Pg - Po;
784  tmpSw = flow[SATNUM[n]]->GetSwByPcow(Pcow);
785  if (gas) {
786  tmpSg = flow[SATNUM[n]]->GetSgByPcgo(Pcgo);
787  }
788  if (tmpSw + tmpSg > 1) {
789  // should me modified
790  OCP_DBL Pcgw = Pcow + Pcgo;
791  tmpSw = flow[SATNUM[n]]->GetSwByPcgw(Pcgw);
792  tmpSg = 1 - tmpSw;
793  }
794  Sw += tmpSw;
795  Sg += tmpSg;
796  }
797  Sw /= ncut;
798  Sg /= ncut;
799  S[n * numPhase + numPhase - 1] = Sw;
800  if (gas) {
801  S[n * numPhase + numPhase - 2] = Sg;
802  }
803 
804  // correct if Pcow is not considered
805  if (!FlagPcow[SATNUM[n]]) {
806  S[n * numPhase + numPhase - 1] = flow[SATNUM[n]]->GetSwco();
807  }
808  }
809 }
OCPTable::IsEmpty
bool IsEmpty() const
judge if table is empty.
Definition: OCPTable.hpp:42
OCPTable::Eval
OCP_DBL Eval(const USI &j, const OCP_DBL &val, const USI &destj)
Definition: OCPTable.cpp:135

References OCPTable::Display(), OCPTable::Eval(), OCPTable::Eval_All(), OCPTable::GetCol(), OCPTable::IsEmpty(), OCP_ABORT, OCP_FUNCNAME, and TINY.

◆ InitSjPcComp()

void Bulk::InitSjPcComp ( const USI &  tabrow,
const Grid &  myGrid 
)

Calculate initial equilibrium for compositional model according to EQUIL.

Here tabrow is maximum number of depth nodes in table of depth vs pressure.

Definition at line 812 of file Bulk.cpp.

813 {
814  OCP_FUNCNAME;
815 
816  OCP_DBL Dref = EQUIL.Dref;
817  OCP_DBL Pref = EQUIL.Pref;
818  OCP_DBL DOWC = EQUIL.DOWC;
819  OCP_DBL PcOW = EQUIL.PcOW;
820  OCP_DBL DOGC = EQUIL.DGOC;
821  OCP_DBL PcGO = EQUIL.PcGO;
822  OCP_DBL Zmin = 1E8;
823  OCP_DBL Zmax = 0;
824 
825  OCP_DBL tmp;
826 
827  for (OCP_USI n = 0; n < numBulk; n++) {
828  OCP_DBL temp1 = depth[n] - dz[n] / 2;
829  OCP_DBL temp2 = depth[n] + dz[n] / 2;
830  Zmin = Zmin < temp1 ? Zmin : temp1;
831  Zmax = Zmax > temp2 ? Zmax : temp2;
832  }
833  OCP_DBL tabdz = (Zmax - Zmin) / (tabrow - 1);
834 
835  // creater table
836  OCPTable DepthP(tabrow, 4);
837  vector<OCP_DBL>& Ztmp = DepthP.GetCol(0);
838  vector<OCP_DBL>& Potmp = DepthP.GetCol(1);
839  vector<OCP_DBL>& Pwtmp = DepthP.GetCol(2);
840  vector<OCP_DBL>& Pgtmp = DepthP.GetCol(3);
841 
842  vector<OCP_DBL> tmpInitZi(numCom, 0);
843 
844  // cal Tab_Ztmp
845  Ztmp[0] = Zmin;
846  for (USI i = 1; i < tabrow; i++) {
847  Ztmp[i] = Ztmp[i - 1] + tabdz;
848  }
849 
850  // find the RefId
851  USI beginId = 0;
852  if (Dref <= Ztmp[0]) {
853  beginId = 0;
854  } else if (Dref >= Ztmp[tabrow - 1]) {
855  beginId = tabrow - 1;
856  } else {
857  beginId =
858  distance(Ztmp.begin(), find_if(Ztmp.begin(), Ztmp.end(),
859  [s = Dref](auto& t) { return t > s; }));
860  beginId--;
861  }
862 
863  // begin calculating oil pressure:
864  OCP_DBL mytemp = T;
865  OCP_DBL gammaOtmp, gammaWtmp, gammaGtmp;
866  OCP_DBL Ptmp;
867  USI mynum = 10;
868  OCP_DBL mydz = 0;
869  OCP_DBL Poref, Pgref, Pwref;
870  OCP_DBL Pbegin = 0;
871 
872  if (Dref < DOGC) {
873  // reference pressure is gas pressure
874  Pgref = Pref;
875  initZi_T[0].Eval_All0(Dref, tmpInitZi);
876  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgref, mytemp, &tmpInitZi[0]);
877  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
878  Pgtmp[beginId] = Pbegin;
879 
880  // find the gas pressure
881  for (USI id = beginId; id > 0; id--) {
882  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
883  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgtmp[id], mytemp, &tmpInitZi[0]);
884  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
885  }
886 
887  for (USI id = beginId; id < tabrow - 1; id++) {
888  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
889  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgtmp[id], mytemp, &tmpInitZi[0]);
890  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
891  }
892 
893  // find the oil pressure in Dref by Pgref
894  Poref = 0;
895  Ptmp = Pgref;
896  mydz = (DOGC - Dref) / mynum;
897  OCP_DBL myz = Dref;
898 
899  for (USI i = 0; i < mynum; i++) {
900  initZi_T[0].Eval_All0(myz, tmpInitZi);
901  myz += mydz;
902  gammaGtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
903  Ptmp += gammaGtmp * mydz;
904  }
905  Ptmp -= PcGO;
906  for (USI i = 0; i < mynum; i++) {
907  initZi_T[0].Eval_All0(myz, tmpInitZi);
908  myz -= mydz;
909  gammaOtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
910  Ptmp -= gammaOtmp * mydz;
911  }
912  Poref = Ptmp;
913 
914  // find the oil pressure in tab
915  initZi_T[0].Eval_All0(Dref, tmpInitZi);
916  gammaOtmp = flashCal[0]->GammaPhaseOG(Poref, Dref, &tmpInitZi[0]);
917  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
918  Potmp[beginId] = Pbegin;
919 
920  for (USI id = beginId; id > 0; id--) {
921  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
922  gammaOtmp = flashCal[0]->GammaPhaseOG(Potmp[id], mytemp, &tmpInitZi[0]);
923  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
924  }
925 
926  for (USI id = beginId; id < tabrow - 1; id++) {
927  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
928  gammaOtmp = flashCal[0]->GammaPhaseOG(Potmp[id], mytemp, &tmpInitZi[0]);
929  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
930  }
931 
932  // find the water pressure in Dref by Poref
933  Pwref = 0;
934  Ptmp = Poref;
935  mydz = (DOWC - Dref) / mynum;
936  myz = Dref;
937 
938  for (USI i = 0; i < mynum; i++) {
939  initZi_T[0].Eval_All0(myz, tmpInitZi);
940  myz += mydz;
941  gammaOtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
942  Ptmp += gammaOtmp * mydz;
943  }
944  Ptmp -= PcOW;
945  for (USI i = 0; i < mynum; i++) {
946  gammaWtmp = flashCal[0]->GammaPhaseW(Ptmp);
947  Ptmp -= gammaWtmp * mydz;
948  }
949  Pwref = Ptmp;
950 
951  // find the water pressure in tab
952  gammaWtmp = flashCal[0]->GammaPhaseW(Pwref);
953  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
954  Pwtmp[beginId] = Pbegin;
955 
956  for (USI id = beginId; id > 0; id--) {
957  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
958  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
959  }
960 
961  for (USI id = beginId; id < tabrow - 1; id++) {
962  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
963  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
964  }
965  } else if (Dref > DOWC) {
966 
967  // reference pressure is water pressure
968  Pwref = Pref;
969  gammaWtmp = flashCal[0]->GammaPhaseW(Pwref);
970  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
971  Pwtmp[beginId] = Pbegin;
972 
973  // find the water pressure
974  for (USI id = beginId; id > 0; id--) {
975  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
976  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
977  }
978  for (USI id = beginId; id < tabrow - 1; id++) {
979  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
980  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
981  }
982 
983  // find the oil pressure in Dref by Pwref
984  Poref = 0;
985  Ptmp = Pwref;
986  mydz = (DOWC - Dref) / mynum;
987  OCP_DBL myz = Dref;
988 
989  for (USI i = 0; i < mynum; i++) {
990  myz += mydz;
991  gammaWtmp = flashCal[0]->GammaPhaseW(Ptmp);
992  Ptmp += gammaWtmp * mydz;
993  }
994  Ptmp += PcOW;
995 
996  for (USI i = 0; i < mynum; i++) {
997  initZi_T[0].Eval_All0(myz, tmpInitZi);
998  myz -= mydz;
999  gammaOtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
1000  Ptmp -= gammaOtmp * mydz;
1001  }
1002  Poref = Ptmp;
1003 
1004  // find the oil pressure in tab
1005  initZi_T[0].Eval_All0(Dref, tmpInitZi);
1006  gammaOtmp = flashCal[0]->GammaPhaseOG(Poref, mytemp, &tmpInitZi[0]);
1007  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
1008  Potmp[beginId] = Pbegin;
1009 
1010  for (USI id = beginId; id > 0; id--) {
1011  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1012  gammaOtmp = flashCal[0]->GammaPhaseOG(Potmp[id], mytemp, &tmpInitZi[0]);
1013  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
1014  }
1015 
1016  for (USI id = beginId; id < tabrow - 1; id++) {
1017  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1018  gammaOtmp = flashCal[0]->GammaPhaseOG(Potmp[id], mytemp, &tmpInitZi[0]);
1019  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
1020  }
1021 
1022  // find the gas pressure in Dref by Poref
1023  Pgref = 0;
1024  Ptmp = Poref;
1025  mydz = (DOGC - Dref) / mynum;
1026  myz = Dref;
1027 
1028  for (USI i = 0; i < mynum; i++) {
1029  initZi_T[0].Eval_All0(myz, tmpInitZi);
1030  myz += mydz;
1031  gammaOtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
1032  Ptmp += gammaOtmp * mydz;
1033  }
1034  Ptmp += PcGO;
1035  for (USI i = 0; i < mynum; i++) {
1036  initZi_T[0].Eval_All0(myz, tmpInitZi);
1037  myz -= mydz;
1038  gammaGtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
1039  Ptmp -= gammaGtmp * mydz;
1040  }
1041  Pgref = Ptmp;
1042 
1043  // find the gas pressure in tab
1044  initZi_T[0].Eval_All0(Dref, tmpInitZi);
1045  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgref, mytemp, &tmpInitZi[0]);
1046  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
1047  Pgtmp[beginId] = Pbegin;
1048 
1049  for (USI id = beginId; id > 0; id--) {
1050  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1051  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgtmp[id], mytemp, &tmpInitZi[0]);
1052  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
1053  }
1054  for (USI id = beginId; id < tabrow - 1; id++) {
1055  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1056  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgtmp[id], mytemp, &tmpInitZi[0]);
1057  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
1058  }
1059  } else {
1060  // reference pressure is oil pressure
1061  Poref = Pref;
1062  initZi_T[0].Eval_All0(Dref, tmpInitZi);
1063  gammaOtmp = flashCal[0]->GammaPhaseOG(Poref, mytemp, &tmpInitZi[0]);
1064  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
1065  Potmp[beginId] = Pbegin;
1066 
1067  // find the oil pressure
1068  for (USI id = beginId; id > 0; id--) {
1069  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1070  gammaOtmp = flashCal[0]->GammaPhaseOG(Potmp[id], mytemp, &tmpInitZi[0]);
1071  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
1072  }
1073  for (USI id = beginId; id < tabrow - 1; id++) {
1074  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1075  gammaOtmp = flashCal[0]->GammaPhaseOG(Potmp[id], mytemp, &tmpInitZi[0]);
1076  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
1077  }
1078 
1079  // find the gas pressure in Dref by Poref
1080  Pgref = 0;
1081  Ptmp = Poref;
1082  mydz = (DOGC - Dref) / mynum;
1083  OCP_DBL myz = Dref;
1084 
1085  for (USI i = 0; i < mynum; i++) {
1086  initZi_T[0].Eval_All0(myz, tmpInitZi);
1087  myz += mydz;
1088  gammaOtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
1089  Ptmp += gammaOtmp * mydz;
1090  }
1091  Ptmp += PcGO;
1092  for (USI i = 0; i < mynum; i++) {
1093  initZi_T[0].Eval_All0(myz, tmpInitZi);
1094  myz -= mydz;
1095  gammaGtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
1096  Ptmp -= gammaGtmp * mydz;
1097  }
1098  Pgref = Ptmp;
1099 
1100  // find the gas pressure in tab
1101  initZi_T[0].Eval_All0(Dref, tmpInitZi);
1102  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgref, mytemp, &tmpInitZi[0]);
1103  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
1104  Pgtmp[beginId] = Pbegin;
1105 
1106  for (USI id = beginId; id > 0; id--) {
1107  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1108  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgtmp[id], mytemp, &tmpInitZi[0]);
1109  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
1110  }
1111 
1112  for (USI id = beginId; id < tabrow - 1; id++) {
1113  initZi_T[0].Eval_All0(Ztmp[id], tmpInitZi);
1114  gammaGtmp = flashCal[0]->GammaPhaseOG(Pgtmp[id], mytemp, &tmpInitZi[0]);
1115  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
1116  }
1117 
1118  // find the water pressure in Dref by Poref
1119  Pwref = 0;
1120  Ptmp = Poref;
1121  mydz = (DOWC - Dref) / mynum;
1122  myz = Dref;
1123 
1124  for (USI i = 0; i < mynum; i++) {
1125  initZi_T[0].Eval_All0(myz, tmpInitZi);
1126  myz += mydz;
1127  gammaOtmp = flashCal[0]->GammaPhaseOG(Ptmp, mytemp, &tmpInitZi[0]);
1128  Ptmp += gammaOtmp * mydz;
1129  }
1130  Ptmp -= PcOW;
1131  for (USI i = 0; i < mynum; i++) {
1132  gammaWtmp = flashCal[0]->GammaPhaseW(Ptmp);
1133  Ptmp -= gammaWtmp * mydz;
1134  }
1135  Pwref = Ptmp;
1136 
1137  // find the water pressure in tab
1138  gammaWtmp = flashCal[0]->GammaPhaseW(Pwref);
1139  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
1140  Pwtmp[beginId] = Pbegin;
1141 
1142  for (USI id = beginId; id > 0; id--) {
1143  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
1144  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
1145  }
1146 
1147  for (USI id = beginId; id < tabrow - 1; id++) {
1148  gammaWtmp = flashCal[0]->GammaPhaseW(Pwtmp[id]);
1149  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
1150  }
1151  }
1152 
1153  cout << "Depth "
1154  << "Poil "
1155  << "Pwat "
1156  << "Pgas" << endl;
1157  DepthP.Display();
1158 
1159  // calculate Pc from DepthP to calculate Sj
1160  std::vector<OCP_DBL> data(4, 0), cdata(4, 0);
1161  // if capillary between water and oil is considered
1162  vector<bool> FlagPcow(NTSFUN, true);
1163  for (USI i = 0; i < NTSFUN; i++) {
1164  if (fabs(flow[i]->GetPcowBySw(0.0 - TINY)) < TINY &&
1165  fabs(flow[i]->GetPcowBySw(1.0 + TINY) < TINY)) {
1166  FlagPcow[i] = false;
1167  }
1168  }
1169 
1170  for (OCP_USI n = 0; n < numBulk; n++) {
1171  initZi_T[0].Eval_All0(depth[n], tmpInitZi);
1172  for (USI i = 0; i < numCom_1; i++) {
1173  initZi[n * numCom + i] = tmpInitZi[i];
1174  }
1175 
1176  DepthP.Eval_All(0, depth[n], data, cdata);
1177  OCP_DBL Po = data[1];
1178  OCP_DBL Pw = data[2];
1179  OCP_DBL Pg = data[3];
1180  OCP_DBL Pcgo = Pg - Po;
1181  OCP_DBL Pcow = Po - Pw;
1182  OCP_DBL Sw = flow[SATNUM[n]]->GetSwByPcow(Pcow);
1183  OCP_DBL Sg = 0;
1184  if (gas) {
1185  Sg = flow[SATNUM[n]]->GetSgByPcgo(Pcgo);
1186  }
1187  if (Sw + Sg > 1) {
1188  // should me modified
1189  OCP_DBL Pcgw = Pcow + Pcgo;
1190  Sw = flow[SATNUM[n]]->GetSwByPcgw(Pcgw);
1191  Sg = 1 - Sw;
1192  }
1193 
1194  if (1 - Sw < TINY) {
1195  // all water
1196  Po = Pw + flow[SATNUM[n]]->GetPcowBySw(1.0);
1197  } else if (1 - Sg < TINY) {
1198  // all gas
1199  Po = Pg - flow[SATNUM[n]]->GetPcgoBySg(1.0);
1200  } else if (1 - Sw - Sg < TINY) {
1201  // water and gas
1202  Po = Pg - flow[SATNUM[n]]->GetPcgoBySg(Sg);
1203  }
1204  P[n] = Po;
1205 
1206  // cal Sw --- Sg is not needed in initialization of Compositional Model
1207  OCP_DBL swco = flow[SATNUM[n]]->GetSwco();
1208  if (!FlagPcow[SATNUM[n]]) {
1209  S[n * numPhase + numPhase - 1] = swco;
1210  continue;
1211  }
1212 
1213  Sw = 0;
1214  Sg = 0;
1215  USI ncut = 10;
1216  OCP_DBL avePcow = 0;
1217 
1218  for (USI k = 0; k < ncut; k++) {
1219  OCP_DBL tmpSw = 0;
1220  OCP_DBL tmpSg = 0;
1221  OCP_DBL dep = depth[n] + dz[n] / ncut * (k - (ncut - 1) / 2.0);
1222  DepthP.Eval_All(0, dep, data, cdata);
1223  Po = data[1];
1224  Pw = data[2];
1225  Pg = data[3];
1226  Pcow = Po - Pw;
1227  Pcgo = Pg - Po;
1228  avePcow += Pcow;
1229  tmpSw = flow[SATNUM[n]]->GetSwByPcow(Pcow);
1230  if (gas) {
1231  tmpSg = flow[SATNUM[n]]->GetSgByPcgo(Pcgo);
1232  }
1233  if (tmpSw + tmpSg > 1) {
1234  // should be modified
1235  OCP_DBL Pcgw = Pcow + Pcgo;
1236  tmpSw = flow[SATNUM[n]]->GetSwByPcgw(Pcgw);
1237  tmpSg = 1 - tmpSw;
1238  }
1239  Sw += tmpSw;
1240  // Sg += tmpSg;
1241  }
1242  Sw /= ncut;
1243  // Sg /= ncut;
1244  avePcow /= ncut;
1245 
1246  if (SwatInitExist) {
1247  if (ScalePcow) {
1248  ScaleValuePcow[n] = 1.0;
1249  }
1250  if (SwatInit[n] <= swco) {
1251  Sw = swco;
1252  } else {
1253  Sw = SwatInit[n];
1254  if (ScalePcow) {
1255  if (avePcow > 0) {
1256  tmp = flow[SATNUM[n]]->GetPcowBySw(Sw);
1257  if (tmp > 0) {
1258  ScaleValuePcow[n] = avePcow / tmp;
1259  /*USI I, J, K;
1260  myGrid.GetIJKGrid(I, J, K, myGrid.activeMap_B2G[n]);
1261  cout << I << " " << J << " " << K << " "
1262  << fixed << setprecision(3) << " " <<
1263  ScaleValuePcow[n] * flow[0]->GetPcowBySw(0) << endl;*/
1264  }
1265  }
1266  }
1267  }
1268  }
1269  S[n * numPhase + numPhase - 1] = Sw;
1270  // if (gas)
1271  //{
1272  // S[n * numPhase + numPhase - 2] = Sg;
1273  // }
1274  }
1275 }

References OCPTable::Display(), OCPTable::Eval_All(), OCPTable::GetCol(), OCP_FUNCNAME, and TINY.

◆ InputParam()

void Bulk::InputParam ( ParamReservoir &  rs_param)

Input param from internal data structure ParamReservoir.

Read parameters from rs_param data structure.

Definition at line 24 of file Bulk.cpp.

25 {
26  OCP_FUNCNAME;
27 
28  rockPref = rs_param.rock.Pref;
29  rockC1 = rs_param.rock.Cr;
30  rockC2 = rockC1;
31  T = rs_param.rsTemp + 460;
32  ScalePcow = rs_param.ScalePcow;
33  blackOil = rs_param.blackOil;
34  comps = rs_param.comps;
35  oil = rs_param.oil;
36  gas = rs_param.gas;
37  water = rs_param.water;
38  disGas = rs_param.disGas;
39  miscible = rs_param.EoSp.miscible;
40  EQUIL.Dref = rs_param.EQUIL[0];
41  EQUIL.Pref = rs_param.EQUIL[1];
42  NTPVT = rs_param.NTPVT;
43  NTSFUN = rs_param.NTSFUN;
44 
45  if (rs_param.PBVD_T.data.size() > 0) EQUIL.PBVD.Setup(rs_param.PBVD_T.data[0]);
46 
47  if (blackOil) {
48 
49  if (water && !oil && !gas) {
50  // water
51  numPhase = 1;
52  numCom = 1;
53  SATmode = PHASE_W;
54  PVTmode = PHASE_W;
55  } else if (water && oil && !gas) {
56  // water, dead oil
57  numPhase = 2;
58  numCom = 2;
59  EQUIL.DOWC = rs_param.EQUIL[2];
60  EQUIL.PcOW = rs_param.EQUIL[3];
61  SATmode = PHASE_OW;
62  PVTmode = PHASE_OW;
63  } else if (water && oil && gas && !disGas) {
64  // water, dead oil, dry gas
65  numPhase = 3;
66  numCom = 3;
67  EQUIL.DOWC = rs_param.EQUIL[2];
68  EQUIL.PcOW = rs_param.EQUIL[3];
69  EQUIL.DGOC = rs_param.EQUIL[4];
70  EQUIL.PcGO = rs_param.EQUIL[5];
71  SATmode = PHASE_DOGW;
72  PVTmode = PHASE_DOGW; // maybe it should be added later
73  } else if (water && oil && gas && disGas) {
74  // water, live oil, dry gas
75  numPhase = 3;
76  numCom = 3;
77  EQUIL.DOWC = rs_param.EQUIL[2];
78  EQUIL.PcOW = rs_param.EQUIL[3];
79  EQUIL.DGOC = rs_param.EQUIL[4];
80  EQUIL.PcGO = rs_param.EQUIL[5];
81  PVTmode = PHASE_ODGW;
82 
83  if (rs_param.SOF3_T.data.size() > 0) {
84  SATmode = PHASE_ODGW02;
85  } else {
86  SATmode = PHASE_ODGW01;
87  }
88  }
89  numCom_1 = numCom - 1;
90  rs_param.numPhase = numPhase;
91  rs_param.numCom = numCom;
92 
93  switch (PVTmode) {
94  case PHASE_W:
95  OCP_ABORT("Wrong Type!");
96  break;
97  case PHASE_OW:
98  for (USI i = 0; i < NTPVT; i++)
99  flashCal.push_back(new BOMixture_OW(rs_param, i));
100  break;
101  case PHASE_DOGW:
102  OCP_ABORT("Wrong Type!");
103  break;
104  case PHASE_ODGW:
105  for (USI i = 0; i < NTPVT; i++)
106  flashCal.push_back(new BOMixture_ODGW(rs_param, i));
107  break;
108  default:
109  OCP_ABORT("Wrong Type!");
110  break;
111  }
112 
113  cout << "Bulk::InputParam --- BLACKOIL" << endl;
114  } else if (comps) {
115 
116  // Water exists and is excluded in EoS model NOW!
117  oil = true;
118  gas = true;
119  water = true;
120 
121  numPhase = rs_param.EoSp.numPhase + 1;
122  numCom = rs_param.EoSp.numCom + 1;
123  numCom_1 = numCom - 1;
124  EQUIL.DOWC = rs_param.EQUIL[2];
125  EQUIL.PcOW = rs_param.EQUIL[3];
126  EQUIL.DGOC = rs_param.EQUIL[4];
127  EQUIL.PcGO = rs_param.EQUIL[5];
128 
129  for (auto& v : rs_param.ZMFVD_T.data) {
130  initZi_T.push_back(OCPTable(v));
131  }
132 
133  if (rs_param.SOF3_T.data.size() > 0) {
134  SATmode = PHASE_ODGW02;
135  } else {
136  SATmode = PHASE_ODGW01;
137  }
138 
139  for (USI i = 0; i < NTPVT; i++)
140  flashCal.push_back(new MixtureComp(rs_param, i));
141 
142  cout << "Bulk::InputParam --- COMPOSITIONAL" << endl;
143  }
144 
145  if (SATmode == PHASE_ODGW01 && miscible) {
146  SATmode = PHASE_ODGW01_MISCIBLE;
147  }
148 
149  satcm.resize(NTSFUN);
150 
151  switch (SATmode) {
152  case PHASE_W:
153  for (USI i = 0; i < NTSFUN; i++)
154  flow.push_back(new FlowUnit_W(rs_param, i));
155  break;
156  case PHASE_OW:
157  for (USI i = 0; i < NTSFUN; i++)
158  flow.push_back(new FlowUnit_OW(rs_param, i));
159  break;
160  case PHASE_ODGW01:
161  for (USI i = 0; i < NTSFUN; i++) {
162  flow.push_back(new FlowUnit_ODGW01(rs_param, i));
163  satcm[i] = flow[i]->GetScm();
164  }
165  break;
166  case PHASE_ODGW01_MISCIBLE:
167  for (USI i = 0; i < NTSFUN; i++) {
168  flow.push_back(new FlowUnit_ODGW01_Miscible(rs_param, i));
169  satcm[i] = flow[i]->GetScm();
170  }
171  break;
172  case PHASE_ODGW02:
173  for (USI i = 0; i < NTSFUN; i++)
174  flow.push_back(new FlowUnit_ODGW02(rs_param, i));
175  break;
176  default:
177  OCP_ABORT("Wrong Type!");
178  break;
179  }
180 }
PHASE_ODGW01
const USI PHASE_ODGW01
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:102
PHASE_ODGW
const USI PHASE_ODGW
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:100
PHASE_ODGW02
const USI PHASE_ODGW02
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:103
PHASE_W
const USI PHASE_W
Phase type = water only.
Definition: OCPConst.hpp:96
PHASE_ODGW01_MISCIBLE
const USI PHASE_ODGW01_MISCIBLE
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:104
PHASE_OW
const USI PHASE_OW
Phase type = oil-water.
Definition: OCPConst.hpp:98
PHASE_DOGW
const USI PHASE_DOGW
Phase type = dead oil-gas-water.
Definition: OCPConst.hpp:101
BOMixture_ODGW
Definition: MixtureBO.hpp:191
BOMixture_OW
Definition: MixtureBO.hpp:140
EoSparam::numCom
USI numCom
num of components, water is excluded.
Definition: ParamReservoir.hpp:94
EoSparam::numPhase
USI numPhase
num of phase, water is excluded, constant now.
Definition: ParamReservoir.hpp:95
EoSparam::miscible
bool miscible
Miscible treatment of hydrocarbons, used in compositional Model.
Definition: ParamReservoir.hpp:114
FlowUnit_ODGW01_Miscible
Definition: FlowUnit.hpp:232
FlowUnit_ODGW01
Definition: FlowUnit.hpp:182
FlowUnit_ODGW02
Definition: FlowUnit.hpp:279
FlowUnit_OW
Definition: FlowUnit.hpp:88
FlowUnit_W
Definition: FlowUnit.hpp:58
OCPTable::Setup
void Setup(const vector< vector< OCP_DBL >> &src)
Setup tables from existing data of table.
Definition: OCPTable.cpp:32
ParamReservoir::NTSFUN
USI NTSFUN
Num of SAT regions.
Definition: ParamReservoir.hpp:187
ParamReservoir::disGas
bool disGas
If true, dissolve gas could exist in oil phase.
Definition: ParamReservoir.hpp:180
ParamReservoir::ZMFVD_T
TableSet ZMFVD_T
Table set of ZMFVD.
Definition: ParamReservoir.hpp:201
ParamReservoir::EQUIL
vector< OCP_DBL > EQUIL
See ParamEQUIL.
Definition: ParamReservoir.hpp:202
ParamReservoir::numPhase
USI numPhase
Number of phases.
Definition: ParamReservoir.hpp:205
ParamReservoir::gas
bool gas
If true, gas phase could exist.
Definition: ParamReservoir.hpp:178
ParamReservoir::SOF3_T
TableSet SOF3_T
Table set of SOF3.
Definition: ParamReservoir.hpp:198
ParamReservoir::PBVD_T
TableSet PBVD_T
Table set of PBVD.
Definition: ParamReservoir.hpp:199
ParamReservoir::ScalePcow
bool ScalePcow
whether Pcow should be scaled.
Definition: ParamReservoir.hpp:168
ParamReservoir::numCom
USI numCom
Number of components(hydrocarbon components), used in Compositional Model when input.
Definition: ParamReservoir.hpp:206
ParamReservoir::oil
bool oil
If true, oil phase could exist.
Definition: ParamReservoir.hpp:177
ParamReservoir::blackOil
bool blackOil
If ture, blackoil model will be used.
Definition: ParamReservoir.hpp:175
ParamReservoir::comps
bool comps
If true, compositional model will be used.
Definition: ParamReservoir.hpp:176
ParamReservoir::NTPVT
USI NTPVT
Num of PVT regions.
Definition: ParamReservoir.hpp:188
ParamReservoir::rsTemp
OCP_DBL rsTemp
Temperature for reservoir.
Definition: ParamReservoir.hpp:159
ParamReservoir::water
bool water
If true, water phase could exist.
Definition: ParamReservoir.hpp:179
ParamReservoir::EoSp
EoSparam EoSp
Initial component composition, used in compositional models.
Definition: ParamReservoir.hpp:184
Rock::Cr
OCP_DBL Cr
Compressibility factor of rock in reservoir.
Definition: ParamReservoir.hpp:52
Rock::Pref
OCP_DBL Pref
Reference pressure at initial porosity.
Definition: ParamReservoir.hpp:51
TableSet::data
vector< vector< vector< OCP_DBL > > > data
All table with the same name.
Definition: ParamReservoir.hpp:35

References ParamReservoir::blackOil, ParamReservoir::comps, Rock::Cr, TableSet::data, ParamReservoir::disGas, ParamReservoir::EoSp, ParamReservoir::EQUIL, ParamReservoir::gas, EoSparam::miscible, ParamReservoir::NTPVT, ParamReservoir::NTSFUN, EoSparam::numCom, ParamReservoir::numCom, EoSparam::numPhase, ParamReservoir::numPhase, OCP_ABORT, OCP_FUNCNAME, ParamReservoir::oil, ParamReservoir::PBVD_T, PHASE_DOGW, PHASE_ODGW, PHASE_ODGW01, PHASE_ODGW01_MISCIBLE, PHASE_ODGW02, PHASE_OW, PHASE_W, Rock::Pref, ParamReservoir::rock, ParamReservoir::rsTemp, ParamReservoir::ScalePcow, OCPTable::Setup(), ParamReservoir::SOF3_T, ParamReservoir::water, and ParamReservoir::ZMFVD_T.

◆ Setup()

void Bulk::Setup ( const Grid &  myGrid)

Allocate memory for bulk data of grid.

Setup bulk information.

Definition at line 183 of file Bulk.cpp.

184 {
185  OCP_FUNCNAME;
186 
187  // Rock/Grid properties
188 
189  numBulk = myGrid.activeGridNum;
190  dx.resize(numBulk, 0);
191  dy.resize(numBulk, 0);
192  dz.resize(numBulk, 0);
193  depth.resize(numBulk, 0);
194  ntg.resize(numBulk, 0);
195  rockVpInit.resize(numBulk, 0);
196  rockVp.resize(numBulk, 0);
197  rockKxInit.resize(numBulk, 0);
198  rockKyInit.resize(numBulk, 0);
199  rockKzInit.resize(numBulk, 0);
200  SATNUM.resize(numBulk, 0);
201  PVTNUM.resize(numBulk, 0);
202 
203  if (myGrid.SwatInit.size() > 0) {
204  SwatInitExist = true;
205  SwatInit.resize(numBulk);
206  }
207  if (ScalePcow) {
208  ScaleValuePcow.resize(numBulk, 0);
209  }
210 
211  for (OCP_USI bIdb = 0; bIdb < numBulk; bIdb++) {
212  OCP_USI bIdg = myGrid.activeMap_B2G[bIdb];
213 
214  dx[bIdb] = myGrid.dx[bIdg];
215  dy[bIdb] = myGrid.dy[bIdg];
216  dz[bIdb] = myGrid.dz[bIdg];
217  depth[bIdb] = myGrid.depth[bIdg];
218  ntg[bIdb] = myGrid.ntg[bIdg];
219 
220  rockVpInit[bIdb] = myGrid.v[bIdg] * myGrid.ntg[bIdg] * myGrid.poro[bIdg];
221  rockKxInit[bIdb] = myGrid.kx[bIdg];
222  rockKyInit[bIdb] = myGrid.ky[bIdg];
223  rockKzInit[bIdb] = myGrid.kz[bIdg];
224 
225  SATNUM[bIdb] = myGrid.SATNUM[bIdg];
226  PVTNUM[bIdb] = myGrid.PVTNUM[bIdg];
227 
228  if (SwatInitExist) {
229  SwatInit[bIdb] = myGrid.SwatInit[bIdg];
230  }
231  }
232 
233  rockVp = rockVpInit;
234  rockKx = rockKxInit;
235  rockKy = rockKyInit;
236  rockKz = rockKzInit;
237 
238  // physical variables
239 
240  P.resize(numBulk);
241  Pb.resize(numBulk);
242  Pj.resize(numBulk * numPhase);
243  Pc.resize(numBulk * numPhase);
244  phaseExist.resize(numBulk * numPhase);
245  S.resize(numBulk * numPhase);
246  rho.resize(numBulk * numPhase);
247  xi.resize(numBulk * numPhase);
248  xij.resize(numBulk * numPhase * numCom);
249  Ni.resize(numBulk * numCom);
250  mu.resize(numBulk * numPhase);
251  kr.resize(numBulk * numPhase);
252  vj.resize(numBulk * numPhase);
253  vf.resize(numBulk);
254  Nt.resize(numBulk);
255 
256  // Phase num
257  phaseNum.resize(numBulk);
258  lphaseNum.resize(numBulk);
259  NRphaseNum.resize(numBulk);
260 
261  phase2Index.resize(3);
262 
263  if (blackOil) {
264  switch (PVTmode) {
265  case PHASE_W:
266  index2Phase.resize(1);
267  index2Phase[0] = WATER;
268  phase2Index[WATER] = 0;
269  break;
270  case PHASE_OW:
271  index2Phase.resize(2);
272  index2Phase[0] = OIL;
273  index2Phase[1] = WATER;
274  phase2Index[OIL] = 0;
275  phase2Index[WATER] = 1;
276  break;
277  case PHASE_OG:
278  index2Phase.resize(2);
279  index2Phase[0] = OIL;
280  index2Phase[1] = GAS;
281  phase2Index[OIL] = 0;
282  phase2Index[GAS] = 1;
283  break;
284  case PHASE_GW:
285  index2Phase.resize(2);
286  index2Phase[0] = GAS;
287  index2Phase[1] = WATER;
288  phase2Index[GAS] = 0;
289  phase2Index[WATER] = 1;
290  break;
291  case PHASE_ODGW:
292  case PHASE_DOGW:
293  index2Phase.resize(3);
294  index2Phase[0] = OIL;
295  index2Phase[1] = GAS;
296  index2Phase[2] = WATER;
297  phase2Index[OIL] = 0;
298  phase2Index[GAS] = 1;
299  phase2Index[WATER] = 2;
300  break;
301  default:
302  OCP_ABORT("Unknown PVT model!");
303  }
304  } else if (comps) {
305  initZi.resize(numBulk * numCom);
306 
307  phase2Index[OIL] = 0;
308  phase2Index[GAS] = 1;
309  phase2Index[WATER] = 2;
310 
311  minEigenSkip.resize(numBulk);
312  flagSkip.resize(numBulk);
313  ziSkip.resize(numBulk * numCom);
314  PSkip.resize(numBulk);
315  lminEigenSkip.resize(numBulk);
316  lflagSkip.resize(numBulk);
317  lziSkip.resize(numBulk * numCom);
318  lPSkip.resize(numBulk);
319  Ks.resize(numBulk * (numCom - 1));
320  lKs.resize(numBulk * (numCom - 1));
321 
322  if (miscible) {
323  surTen.resize(numBulk);
324  Fk.resize(numBulk);
325  Fp.resize(numBulk);
326  lsurTen.resize(numBulk);
327  }
328  }
329 
330  // error
331  ePEC.resize(numBulk);
332  eN.resize(numBulk);
333  eV.resize(numBulk);
334 
335  CalSomeInfo(myGrid);
336 
337 #if DEBUG
338  CheckSetup();
339 #endif
340 }
GAS
const USI GAS
Fluid type = gas.
Definition: OCPConst.hpp:83
OIL
const USI OIL
Fluid type = oil.
Definition: OCPConst.hpp:82
PHASE_OG
const USI PHASE_OG
Phase type = oil-gas.
Definition: OCPConst.hpp:99
PHASE_GW
const USI PHASE_GW
Phase type = gas-water.
Definition: OCPConst.hpp:97
Bulk::CheckSetup
void CheckSetup() const
Check if error occurs in Setup.
Definition: Bulk.cpp:342

References CheckSetup(), GAS, OCP_ABORT, OCP_FUNCNAME, OIL, PHASE_DOGW, PHASE_GW, PHASE_ODGW, PHASE_OG, PHASE_OW, PHASE_W, and WATER.

The documentation for this class was generated from the following files:
Generated on Sat Oct 15 2022 00:08:52 for OpenCAEPoro by doxygen 1.9.1