OpenCAEPoro  v0.5.0
A simulator for multicomponent porous media flow
Bulk.cpp
Go to the documentation of this file.
1 
12 #include <algorithm>
13 #include <cmath>
14 #include <ctime>
15 
16 // OpenCAEPoro header files
17 #include "Bulk.hpp"
18 
20 // HeatLoss
22 
23 void HeatLoss::InputParam(const HLoss& loss)
24 {
25  ifHLoss = loss.ifHLoss;
26  if (ifHLoss) {
27  obC = loss.obC;
28  obK = loss.obK;
29  ubC = loss.ubC;
30  ubK = loss.ubK;
31  }
32 }
33 
34 void HeatLoss::Setup(const OCP_USI& nb)
35 {
36  if (ifHLoss) {
37  obD = obK / obC;
38  ubD = ubK / ubC;
39  numBulk = nb;
40  I.resize(numBulk);
41  p.resize(numBulk);
42  pT.resize(numBulk);
43  lI.resize(numBulk);
44  lp.resize(numBulk);
45  lpT.resize(numBulk);
46  }
47 }
48 
49 void HeatLoss::CalHeatLoss(const vector<USI>& location,
50  const vector<OCP_DBL>& T,
51  const vector<OCP_DBL>& lT,
52  const vector<OCP_DBL>& initT,
53  const OCP_DBL& t,
54  const OCP_DBL& dt)
55 {
56  if (ifHLoss) {
57  OCP_DBL lambda, d, dT, theta;
58  OCP_DBL tmp, q;
59  for (OCP_USI n = 0; n < numBulk; n++) {
60  if (location[n] > 0) {
61  // overburden or underburden
62  lambda = location[n] == 1 ? obD : ubD;
63 
64  dT = T[n] - lT[n];
65  theta = T[n] - initT[n];
66  d = sqrt(lambda * t) / 2;
67  tmp = 3 * pow(d, 2) + lambda * dt;
68  p[n] = (theta * (lambda * dt / d) + lI[n] -
69  dT * (pow(d, 3) / (lambda * dt))) /
70  tmp;
71  pT[n] = (lambda * dt / d - pow(d, 3) / (lambda * dt)) / tmp;
72  q = (2 * p[n] * d - theta + pow(d, 2) * dT / (lambda * dt)) /
73  (2 * pow(d, 2));
74  I[n] = theta * d + p[n] * pow(d, 2) + 2 * q * pow(d, 3);
75  }
76  }
77  }
78 }
79 
80 void HeatLoss::ResetToLastTimeStep()
81 {
82  I = lI;
83  p = lp;
84  pT = lpT;
85 }
86 
87 void HeatLoss::UpdateLastTimeStep()
88 {
89  lI = I;
90  lp = p;
91  lpT = pT;
92 }
93 
95 // Input Param and Setup
97 
99 void Bulk::InputParam(const ParamReservoir& rs_param)
100 {
101  OCP_FUNCNAME;
102 
103  // Common input
104  ifBlackOil = rs_param.blackOil;
105  ifComps = rs_param.comps;
106  ifThermal = rs_param.thermal;
107 
108  if (ifThermal) {
109  ifBlackOil = OCP_FALSE;
110  ifComps = OCP_FALSE;
111  }
112 
113  NTPVT = rs_param.NTPVT;
114  NTSFUN = rs_param.NTSFUN;
115  NTROCC = rs_param.NTROOC;
116 
117  if (rs_param.PBVD_T.data.size() > 0) EQUIL.PBVD.Setup(rs_param.PBVD_T.data[0]);
118 
119  if (ifBlackOil) {
120  // Isothermal blackoil model
121  InputParamBLKOIL(rs_param);
122  } else if (ifComps) {
123  // Isothermal compositional model
124  InputParamCOMPS(rs_param);
125  } else if (ifThermal) {
126  // ifThermal model
127  InputParamTHERMAL(rs_param);
128  }
129 
130  InputSatFunc(rs_param);
131 }
132 
133 void Bulk::InputParamBLKOIL(const ParamReservoir& rs_param)
134 {
135  oil = rs_param.oil;
136  gas = rs_param.gas;
137  water = rs_param.water;
138  disGas = rs_param.disGas;
139 
140  EQUIL.Dref = rs_param.EQUIL[0];
141  EQUIL.Pref = rs_param.EQUIL[1];
142 
143  if (water && !oil && !gas) {
144  // water
145  numPhase = 1;
146  numCom = 1;
147  SATmode = PHASE_W;
148  PVTmodeB = PHASE_W;
149  } else if (water && oil && !gas) {
150  // water, dead oil
151  numPhase = 2;
152  numCom = 2;
153  EQUIL.DOWC = rs_param.EQUIL[2];
154  EQUIL.PcOW = rs_param.EQUIL[3];
155  SATmode = PHASE_OW;
156  PVTmodeB = PHASE_OW;
157  } else if (water && oil && gas && !disGas) {
158  // water, dead oil, dry gas
159  numPhase = 3;
160  numCom = 3;
161  EQUIL.DOWC = rs_param.EQUIL[2];
162  EQUIL.PcOW = rs_param.EQUIL[3];
163  EQUIL.DGOC = rs_param.EQUIL[4];
164  EQUIL.PcGO = rs_param.EQUIL[5];
165  SATmode = PHASE_DOGW;
166  PVTmodeB = PHASE_DOGW; // maybe it should be added later
167  } else if (water && oil && gas && disGas) {
168  // water, live oil, dry gas
169  numPhase = 3;
170  numCom = 3;
171 
172  EQUIL.DOWC = rs_param.EQUIL[2];
173  EQUIL.PcOW = rs_param.EQUIL[3];
174  EQUIL.DGOC = rs_param.EQUIL[4];
175  EQUIL.PcGO = rs_param.EQUIL[5];
176  PVTmodeB = PHASE_ODGW;
177 
178  if (rs_param.SOF3_T.data.size() > 0) {
179  SATmode = PHASE_ODGW02;
180  } else {
181  SATmode = PHASE_ODGW01;
182  }
183  }
184  numComH = numCom - 1;
185 
186  // PVT mode
187  switch (PVTmodeB) {
188  case PHASE_W:
189  OCP_ABORT("Wrong Type!");
190  break;
191  case PHASE_OW:
192  for (USI i = 0; i < NTPVT; i++)
193  flashCal.push_back(new BOMixture_OW(rs_param, i));
194  break;
195  case PHASE_DOGW:
196  OCP_ABORT("Wrong Type!");
197  break;
198  case PHASE_ODGW:
199  for (USI i = 0; i < NTPVT; i++)
200  flashCal.push_back(new BOMixture_ODGW(rs_param, i));
201  break;
202  default:
203  OCP_ABORT("Wrong Type!");
204  break;
205  }
206 
207  phase2Index.resize(3);
208  switch (flashCal[0]->GetMixtureType()) {
209  case BLKOIL_W:
210  phase2Index[WATER] = 0;
211  break;
212  case BLKOIL_OW:
213  phase2Index[OIL] = 0;
214  phase2Index[WATER] = 1;
215  break;
216  case BLKOIL_OG:
217  phase2Index[OIL] = 0;
218  phase2Index[GAS] = 1;
219  break;
220  case BLKOIL_DOGW:
221  case BLKOIL_ODGW:
222  phase2Index[OIL] = 0;
223  phase2Index[GAS] = 1;
224  phase2Index[WATER] = 2;
225  break;
226  default:
227  OCP_ABORT("WRONG Mixture Type!");
228  }
229 
230  InputRockFunc(rs_param);
231  cout << endl << "BLACKOIL model is selected" << endl;
232 }
233 
234 void Bulk::InputParamCOMPS(const ParamReservoir& rs_param)
235 {
236 
237  // Water exists and is excluded in EoS model NOW!
238  oil = OCP_TRUE;
239  gas = OCP_TRUE;
240  water = OCP_TRUE;
241  ifUseEoS = OCP_TRUE;
242 
243  numPhase = rs_param.comsParam.numPhase + 1;
244  numCom = rs_param.comsParam.numCom + 1;
245  numComH = numCom - 1;
246  EQUIL.Dref = rs_param.EQUIL[0];
247  EQUIL.Pref = rs_param.EQUIL[1];
248  EQUIL.DOWC = rs_param.EQUIL[2];
249  EQUIL.PcOW = rs_param.EQUIL[3];
250  EQUIL.DGOC = rs_param.EQUIL[4];
251  EQUIL.PcGO = rs_param.EQUIL[5];
252 
253  // Init Zi
254  for (auto& v : rs_param.ZMFVD_T.data) {
255  initZi_Tab.push_back(OCPTable(v));
256  }
257 
258  // Init T
259  // Use RTEMP
260  rsTemp = rs_param.rsTemp;
261  vector<vector<OCP_DBL>> temp;
262  temp.resize(2);
263  // add depth
264  temp[0].push_back(0);
265  temp[0].push_back(1E8);
266  // add temperature
267  temp[1].push_back(rsTemp);
268  temp[1].push_back(rsTemp);
269  initT_Tab.push_back(OCPTable(temp));
270 
271  // Saturation mode
272  if (rs_param.SOF3_T.data.size() > 0) {
273  SATmode = PHASE_ODGW02;
274  } else {
275  SATmode = PHASE_ODGW01;
276  if (rs_param.comsParam.miscible) {
277  SATmode = PHASE_ODGW01_MISCIBLE;
278  }
279  }
280 
281  // PVT mode
282  for (USI i = 0; i < NTPVT; i++) flashCal.push_back(new MixtureComp(rs_param, i));
283 
284  // phase index
285  phase2Index.resize(3);
286  phase2Index[OIL] = 0;
287  phase2Index[GAS] = 1;
288  phase2Index[WATER] = 2;
289 
290  InputRockFunc(rs_param);
291  cout << endl << "COMPOSITIONAL model is selected" << endl;
292 }
293 
294 void Bulk::InputParamTHERMAL(const ParamReservoir& rs_param)
295 {
296  oil = rs_param.oil;
297  gas = rs_param.gas;
298  water = rs_param.water;
299 
300  // Init T
301  rsTemp = rs_param.rsTemp;
302  for (auto& v : rs_param.TEMPVD_T.data) {
303  initT_Tab.push_back(OCPTable(v));
304  }
305  if (initT_Tab.size() == 0) {
306  // Use RTEMP
307  vector<vector<OCP_DBL>> temp;
308  temp.resize(2);
309  // add depth
310  temp[0].push_back(0);
311  temp[0].push_back(1E8);
312  // add temperature
313  temp[1].push_back(rsTemp);
314  temp[1].push_back(rsTemp);
315  initT_Tab.push_back(OCPTable(temp));
316  }
317  // ifThermal conductivity
318  if (oil) {
319  thconp.push_back(rs_param.thcono);
320  }
321  if (gas) {
322  thconp.push_back(rs_param.thcong);
323  }
324  if (water) {
325  thconp.push_back(rs_param.thconw);
326  }
327 
328  // Only Now
329  SATmode = PHASE_OW;
330  numPhase = 2;
331  numCom = 2;
332  EQUIL.Dref = rs_param.EQUIL[0];
333  EQUIL.Pref = rs_param.EQUIL[1];
334  EQUIL.DOWC = rs_param.EQUIL[2];
335  EQUIL.PcOW = rs_param.EQUIL[3];
336  EQUIL.DGOC = rs_param.EQUIL[4];
337  EQUIL.PcGO = rs_param.EQUIL[5];
338 
339  // PVT mode
340  for (USI i = 0; i < NTPVT; i++)
341  flashCal.push_back(new MixtureThermal_K01(rs_param, i));
342 
343  // phase index
344  phase2Index.resize(3);
345  phase2Index[OIL] = 0;
346  phase2Index[WATER] = 1;
347 
348  InputRockFuncT(rs_param);
349 }
350 
351 void Bulk::InputSatFunc(const ParamReservoir& rs_param)
352 {
353  // Setup Saturation function
354  satcm.resize(NTSFUN);
355  switch (SATmode) {
356  case PHASE_W:
357  for (USI i = 0; i < NTSFUN; i++)
358  flow.push_back(new FlowUnit_W(rs_param, i));
359  break;
360  case PHASE_OW:
361  for (USI i = 0; i < NTSFUN; i++)
362  flow.push_back(new FlowUnit_OW(rs_param, i));
363  break;
364  case PHASE_ODGW01:
365  for (USI i = 0; i < NTSFUN; i++) {
366  flow.push_back(new FlowUnit_ODGW01(rs_param, i));
367  satcm[i] = flow[i]->GetScm();
368  }
369  break;
370  case PHASE_ODGW01_MISCIBLE:
371  for (USI i = 0; i < NTSFUN; i++) {
372  flow.push_back(new FlowUnit_ODGW01_Miscible(rs_param, i));
373  satcm[i] = flow[i]->GetScm();
374  }
375  break;
376  case PHASE_ODGW02:
377  for (USI i = 0; i < NTSFUN; i++)
378  flow.push_back(new FlowUnit_ODGW02(rs_param, i));
379  break;
380  default:
381  OCP_ABORT("Wrong Type!");
382  break;
383  }
384 }
385 
386 void Bulk::InputRockFunc(const ParamReservoir& rs_param)
387 {
388  for (USI i = 0; i < NTROCC; i++) {
389  rock.push_back(new Rock_Linear(rs_param.rockSet[i]));
390  }
391 }
392 
393 void Bulk::InputRockFuncT(const ParamReservoir& rs_param)
394 {
395  for (USI i = 0; i < NTROCC; i++) {
396  if (rs_param.rockSet[i].type == "LINEAR") {
397  rock.push_back(new RockT_Linear(rs_param.rockSet[i]));
398  } else {
399  rock.push_back(new RockT_Exp(rs_param.rockSet[i]));
400  }
401  }
402 
403  // Heat Loss
404  hLoss.InputParam(rs_param.hLoss);
405 }
406 
408 void Bulk::SetupIsoT(const Grid& myGrid)
409 {
410  OCP_FUNCNAME;
411 
412  numBulk = myGrid.activeGridNum;
413  AllocateGridRockIsoT(myGrid);
414  AllocateRegion(myGrid);
415  AllocateError();
416 }
417 
419 void Bulk::SetupT(const Grid& myGrid)
420 {
421  numBulk = myGrid.activeGridNum;
422  AllocateGridRockT(myGrid);
423  AllocateRegion(myGrid);
424  SetupBulkType(myGrid);
425  // Setup Heat Loss
426  hLoss.Setup(numBulk);
427 }
428 
429 void Bulk::SetupOptionalFeatures(const Grid& myGrid, OptionalFeatures& optFeatures)
430 {
431  for (USI i = 0; i < NTPVT; i++) {
432  flashCal[i]->SetupOptionalFeatures(optFeatures, numBulk);
433  }
434  for (USI i = 0; i < NTSFUN; i++) {
435  flow[i]->SetupOptionalFeatures(myGrid, optFeatures);
436  }
437 }
438 
440 // Initial Properties
442 
443 void Bulk::InitPTSw(const USI& tabrow)
444 {
445  OCP_FUNCNAME;
446 
447  initT.resize(numBulk);
448 
449  OCP_DBL Dref = EQUIL.Dref;
450  OCP_DBL Pref = EQUIL.Pref;
451  OCP_DBL DOWC = EQUIL.DOWC;
452  OCP_DBL PcOW = EQUIL.PcOW;
453  OCP_DBL DOGC = EQUIL.DGOC;
454  OCP_DBL PcGO = EQUIL.PcGO;
455  OCP_DBL Zmin = 1E8;
456  OCP_DBL Zmax = 0;
457 
458  for (OCP_USI n = 0; n < numBulk; n++) {
459  OCP_DBL temp1 = depth[n] - dz[n] / 2;
460  OCP_DBL temp2 = depth[n] + dz[n] / 2;
461  Zmin = Zmin < temp1 ? Zmin : temp1;
462  Zmax = Zmax > temp2 ? Zmax : temp2;
463  }
464  OCP_DBL tabdz = (Zmax - Zmin) / (tabrow - 1);
465 
466  // create table
467  OCPTable DepthP(tabrow, 4);
468  vector<OCP_DBL>& Ztmp = DepthP.GetCol(0);
469  vector<OCP_DBL>& Potmp = DepthP.GetCol(1);
470  vector<OCP_DBL>& Pgtmp = DepthP.GetCol(2);
471  vector<OCP_DBL>& Pwtmp = DepthP.GetCol(3);
472 
473  vector<OCP_DBL> tmpInitZi(numCom, 0);
474 
475  // cal Tab_Ztmp
476  Ztmp[0] = Zmin;
477  for (USI i = 1; i < tabrow; i++) {
478  Ztmp[i] = Ztmp[i - 1] + tabdz;
479  }
480 
481  OCP_DBL myTemp = rsTemp;
482 
483  // find the RefId
484  USI beginId = 0;
485  if (Dref <= Ztmp[0]) {
486  beginId = 0;
487  } else if (Dref >= Ztmp[tabrow - 1]) {
488  beginId = tabrow - 1;
489  } else {
490  beginId =
491  distance(Ztmp.begin(), find_if(Ztmp.begin(), Ztmp.end(),
492  [s = Dref](auto& t) { return t > s; }));
493  beginId--;
494  }
495 
496  // begin calculating oil pressure:
497  OCP_DBL Pbb = Pref;
498  OCP_DBL gammaOtmp, gammaWtmp, gammaGtmp;
499  OCP_DBL Ptmp;
500  USI mynum = 10;
501  OCP_DBL mydz = 0;
502  OCP_DBL Poref, Pgref, Pwref;
503  OCP_DBL Pbegin = 0;
504 
505  const OCP_BOOL initZi_flag = initZi_Tab.size() > 0 ? OCP_TRUE : OCP_FALSE;
506  const OCP_BOOL initT_flag = initT_Tab.size() > 0 ? OCP_TRUE : OCP_FALSE;
507  const OCP_BOOL PBVD_flag = EQUIL.PBVD.IsEmpty() ? OCP_FALSE : OCP_TRUE;
508 
509  if (Dref < DOGC) {
510  // reference pressure is gas pressure
511  Pgref = Pref;
512  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
513  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
514  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
515 
516  gammaGtmp = GRAVITY_FACTOR *
517  flashCal[0]->RhoPhase(Pgref, Pbb, myTemp, tmpInitZi.data(), GAS);
518  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
519  Pgtmp[beginId] = Pbegin;
520 
521  // find the gas pressure
522  for (USI id = beginId; id > 0; id--) {
523  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
524  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
525  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
526 
527  gammaGtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgtmp[id], Pbb, myTemp,
528  tmpInitZi.data(), GAS);
529  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
530  }
531 
532  for (USI id = beginId; id < tabrow - 1; id++) {
533  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
534  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
535  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
536 
537  gammaGtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgtmp[id], Pbb, myTemp,
538  tmpInitZi.data(), GAS);
539  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
540  }
541 
542  // find the oil pressure in Dref by Pgref
543  Poref = 0;
544  Ptmp = Pgref;
545  mydz = (DOGC - Dref) / mynum;
546  OCP_DBL myz = Dref;
547 
548  for (USI i = 0; i < mynum; i++) {
549  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
550  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
551  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
552 
553  gammaGtmp = GRAVITY_FACTOR *
554  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), GAS);
555  Ptmp += gammaGtmp * mydz;
556  myz += mydz;
557  }
558  Ptmp -= PcGO;
559  for (USI i = 0; i < mynum; i++) {
560  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
561  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
562  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
563 
564  gammaOtmp = GRAVITY_FACTOR *
565  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), OIL);
566  Ptmp -= gammaOtmp * mydz;
567  myz -= mydz;
568  }
569  Poref = Ptmp;
570 
571  // find the oil pressure in tab
572  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
573  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
574  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
575 
576  gammaOtmp = GRAVITY_FACTOR *
577  flashCal[0]->RhoPhase(Poref, Pbb, myTemp, tmpInitZi.data(), OIL);
578  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
579  Potmp[beginId] = Pbegin;
580 
581  for (USI id = beginId; id > 0; id--) {
582  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
583  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
584  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
585 
586  gammaOtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Potmp[id], Pbb, myTemp,
587  tmpInitZi.data(), OIL);
588  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
589  }
590 
591  for (USI id = beginId; id < tabrow - 1; id++) {
592  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
593  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
594  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
595 
596  gammaOtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Potmp[id], Pbb, myTemp,
597  tmpInitZi.data(), OIL);
598  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
599  }
600 
601  // find the water pressure in Dref by Poref
602  Pwref = 0;
603  Ptmp = Poref;
604  mydz = (DOWC - Dref) / mynum;
605  myz = Dref;
606 
607  for (USI i = 0; i < mynum; i++) {
608  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
609  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
610  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
611 
612  gammaOtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Poref, Pbb, myTemp,
613  tmpInitZi.data(), OIL);
614  Ptmp += gammaOtmp * mydz;
615  myz += mydz;
616  }
617  Ptmp -= PcOW;
618  for (USI i = 0; i < mynum; i++) {
619  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
620  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
621 
622  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp,
623  tmpInitZi.data(), WATER);
624  Ptmp -= gammaWtmp * mydz;
625  myz -= mydz;
626  }
627  Pwref = Ptmp;
628 
629  // find the water pressure in tab
630  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
631  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
632 
633  gammaWtmp = GRAVITY_FACTOR *
634  flashCal[0]->RhoPhase(Pwref, Pbb, myTemp, tmpInitZi.data(), WATER);
635  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
636  Pwtmp[beginId] = Pbegin;
637 
638  for (USI id = beginId; id > 0; id--) {
639  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
640  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
641 
642  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pwtmp[id], Pbb, myTemp,
643  tmpInitZi.data(), WATER);
644  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
645  }
646 
647  for (USI id = beginId; id < tabrow - 1; id++) {
648  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
649  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
650 
651  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pwtmp[id], Pbb, myTemp,
652  tmpInitZi.data(), WATER);
653  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
654  }
655  } else if (Dref > DOWC) {
656  OCP_DBL myz;
657  // reference pressure is water pressure
658  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
659  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
660 
661  Pwref = Pref;
662  gammaWtmp = GRAVITY_FACTOR *
663  flashCal[0]->RhoPhase(Pwref, Pbb, myTemp, tmpInitZi.data(), WATER);
664  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
665  Pwtmp[beginId] = Pbegin;
666 
667  // find the water pressure
668  for (USI id = beginId; id > 0; id--) {
669  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
670  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
671 
672  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pwtmp[id], Pbb, myTemp,
673  tmpInitZi.data(), WATER);
674  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
675  }
676  for (USI id = beginId; id < tabrow - 1; id++) {
677  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
678  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
679 
680  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pwtmp[id], Pbb, myTemp,
681  tmpInitZi.data(), WATER);
682  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
683  }
684 
685  // find the oil pressure in Dref by Pwref
686  Poref = 0;
687  Ptmp = Pwref;
688  mydz = (DOWC - Dref) / mynum;
689  myz = Dref;
690 
691  for (USI i = 0; i < mynum; i++) {
692  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
693  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
694 
695  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp,
696  tmpInitZi.data(), WATER);
697  Ptmp += gammaWtmp * mydz;
698  myz += mydz;
699  }
700  Ptmp += PcOW;
701 
702  for (USI i = 0; i < mynum; i++) {
703  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
704  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
705  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
706 
707  gammaOtmp = GRAVITY_FACTOR *
708  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), OIL);
709  Ptmp -= gammaOtmp * mydz;
710  myz -= mydz;
711  }
712  Poref = Ptmp;
713 
714  // find the oil pressure in tab
715  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
716  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
717  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
718 
719  gammaOtmp = GRAVITY_FACTOR *
720  flashCal[0]->RhoPhase(Poref, Pbb, myTemp, tmpInitZi.data(), OIL);
721  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
722  Potmp[beginId] = Pbegin;
723 
724  for (USI id = beginId; id > 0; id--) {
725  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
726  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
727  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
728 
729  gammaOtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Potmp[id], Pbb, myTemp,
730  tmpInitZi.data(), OIL);
731  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
732  }
733 
734  for (USI id = beginId; id < tabrow - 1; id++) {
735  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
736  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
737  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
738 
739  gammaOtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Potmp[id], Pbb, myTemp,
740  tmpInitZi.data(), OIL);
741  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
742  }
743 
744  if (gas) {
745  // find the gas pressure in Dref by Poref
746  Pgref = 0;
747  Ptmp = Poref;
748  mydz = (DOGC - Dref) / mynum;
749  myz = Dref;
750 
751  for (USI i = 0; i < mynum; i++) {
752  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
753  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
754  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
755 
756  gammaOtmp =
757  GRAVITY_FACTOR *
758  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), OIL);
759  Ptmp += gammaOtmp * mydz;
760  myz += mydz;
761  }
762  Ptmp += PcGO;
763  for (USI i = 0; i < mynum; i++) {
764  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
765  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
766  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
767 
768  gammaGtmp =
769  GRAVITY_FACTOR *
770  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), GAS);
771  Ptmp -= gammaGtmp * mydz;
772  myz -= mydz;
773  }
774  Pgref = Ptmp;
775 
776  // find the gas pressure in tab
777  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
778  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
779  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
780 
781  gammaGtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgref, Pbb, myTemp,
782  tmpInitZi.data(), GAS);
783  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
784  Pgtmp[beginId] = Pbegin;
785 
786  for (USI id = beginId; id > 0; id--) {
787  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
788  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
789  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
790 
791  gammaGtmp =
792  GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgtmp[id], Pbb, myTemp,
793  tmpInitZi.data(), GAS);
794  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
795  }
796  for (USI id = beginId; id < tabrow - 1; id++) {
797  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
798  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
799  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
800 
801  gammaGtmp =
802  GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgtmp[id], Pbb, myTemp,
803  tmpInitZi.data(), GAS);
804  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
805  }
806  }
807 
808  } else {
809  OCP_DBL myz;
810  // reference pressure is oil pressure
811  Poref = Pref;
812  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
813  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
814  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
815 
816  gammaOtmp = GRAVITY_FACTOR *
817  flashCal[0]->RhoPhase(Poref, Pbb, myTemp, tmpInitZi.data(), OIL);
818  Pbegin = Poref + gammaOtmp * (Ztmp[beginId] - Dref);
819  Potmp[beginId] = Pbegin;
820 
821  // find the oil pressure
822  for (USI id = beginId; id > 0; id--) {
823  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
824  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
825  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
826 
827  gammaOtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Potmp[id], Pbb, myTemp,
828  tmpInitZi.data(), OIL);
829  Potmp[id - 1] = Potmp[id] - gammaOtmp * (Ztmp[id] - Ztmp[id - 1]);
830  }
831  for (USI id = beginId; id < tabrow - 1; id++) {
832  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
833  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
834  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
835 
836  gammaOtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Potmp[id], Pbb, myTemp,
837  tmpInitZi.data(), OIL);
838  Potmp[id + 1] = Potmp[id] + gammaOtmp * (Ztmp[id + 1] - Ztmp[id]);
839  }
840 
841  if (gas) {
842  // find the gas pressure in Dref by Poref
843  Pgref = 0;
844  Ptmp = Poref;
845  mydz = (DOGC - Dref) / mynum;
846  myz = Dref;
847 
848  for (USI i = 0; i < mynum; i++) {
849  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
850  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
851  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
852 
853  gammaOtmp =
854  GRAVITY_FACTOR *
855  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), OIL);
856  Ptmp += gammaOtmp * mydz;
857  myz += mydz;
858  }
859  Ptmp += PcGO;
860  for (USI i = 0; i < mynum; i++) {
861  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
862  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
863  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
864 
865  gammaGtmp =
866  GRAVITY_FACTOR *
867  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), GAS);
868  Ptmp -= gammaGtmp * mydz;
869  myz -= mydz;
870  }
871  Pgref = Ptmp;
872 
873  // find the gas pressure in tab
874  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
875  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
876  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Dref, 1);
877 
878  gammaGtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgref, Pbb, myTemp,
879  tmpInitZi.data(), GAS);
880  Pbegin = Pgref + gammaGtmp * (Ztmp[beginId] - Dref);
881  Pgtmp[beginId] = Pbegin;
882 
883  for (USI id = beginId; id > 0; id--) {
884  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
885  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
886  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
887 
888  gammaGtmp =
889  GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgtmp[id], Pbb, myTemp,
890  tmpInitZi.data(), GAS);
891  Pgtmp[id - 1] = Pgtmp[id] - gammaGtmp * (Ztmp[id] - Ztmp[id - 1]);
892  }
893 
894  for (USI id = beginId; id < tabrow - 1; id++) {
895  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
896  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
897  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, Ztmp[id], 1);
898 
899  gammaGtmp =
900  GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pgtmp[id], Pbb, myTemp,
901  tmpInitZi.data(), GAS);
902  Pgtmp[id + 1] = Pgtmp[id] + gammaGtmp * (Ztmp[id + 1] - Ztmp[id]);
903  }
904  }
905 
906  // find the water pressure in Dref by Poref
907  Pwref = 0;
908  Ptmp = Poref;
909  mydz = (DOWC - Dref) / mynum;
910  myz = Dref;
911 
912  for (USI i = 0; i < mynum; i++) {
913  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
914  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
915  if (PBVD_flag) Pbb = EQUIL.PBVD.Eval(0, myz, 1);
916 
917  gammaOtmp = GRAVITY_FACTOR *
918  flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp, tmpInitZi.data(), OIL);
919  Ptmp += gammaOtmp * mydz;
920  myz += mydz;
921  }
922  Ptmp -= PcOW;
923  for (USI i = 0; i < mynum; i++) {
924  if (initZi_flag) initZi_Tab[0].Eval_All0(myz, tmpInitZi);
925  if (initT_flag) myTemp = initT_Tab[0].Eval(0, myz, 1);
926 
927  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Ptmp, Pbb, myTemp,
928  tmpInitZi.data(), WATER);
929  Ptmp -= gammaWtmp * mydz;
930  myz -= mydz;
931  }
932  Pwref = Ptmp;
933 
934  // find the water pressure in tab
935  if (initZi_flag) initZi_Tab[0].Eval_All0(Dref, tmpInitZi);
936  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Dref, 1);
937 
938  gammaWtmp = GRAVITY_FACTOR *
939  flashCal[0]->RhoPhase(Pwref, Pbb, myTemp, tmpInitZi.data(), WATER);
940  Pbegin = Pwref + gammaWtmp * (Ztmp[beginId] - Dref);
941  Pwtmp[beginId] = Pbegin;
942 
943  for (USI id = beginId; id > 0; id--) {
944  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
945  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
946 
947  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pwtmp[id], Pbb, myTemp,
948  tmpInitZi.data(), WATER);
949  Pwtmp[id - 1] = Pwtmp[id] - gammaWtmp * (Ztmp[id] - Ztmp[id - 1]);
950  }
951 
952  for (USI id = beginId; id < tabrow - 1; id++) {
953  if (initZi_flag) initZi_Tab[0].Eval_All0(Ztmp[id], tmpInitZi);
954  if (initT_flag) myTemp = initT_Tab[0].Eval(0, Ztmp[id], 1);
955 
956  gammaWtmp = GRAVITY_FACTOR * flashCal[0]->RhoPhase(Pwtmp[id], Pbb, myTemp,
957  tmpInitZi.data(), WATER);
958  Pwtmp[id + 1] = Pwtmp[id] + gammaWtmp * (Ztmp[id + 1] - Ztmp[id]);
959  }
960  }
961 
962  DepthP.Display();
963 
964  // calculate Pc from DepthP to calculate Sj
965  std::vector<OCP_DBL> data(4, 0), cdata(4, 0);
966  // whether capillary between water and oil is considered
967  vector<OCP_BOOL> FlagPcow(NTSFUN, OCP_TRUE);
968  for (USI i = 0; i < NTSFUN; i++) {
969  if (fabs(flow[i]->GetPcowBySw(0.0 - TINY)) < TINY &&
970  fabs(flow[i]->GetPcowBySw(1.0 + TINY) < TINY)) {
971  FlagPcow[i] = OCP_FALSE;
972  }
973  }
974 
975  for (OCP_USI n = 0; n < numBulk; n++) {
976  if (initZi_flag) {
977  initZi_Tab[0].Eval_All0(depth[n], tmpInitZi);
978  for (USI i = 0; i < numComH; i++) {
979  Ni[n * numCom + i] = tmpInitZi[i];
980  }
981  }
982  if (initT_flag) {
983  myTemp = initT_Tab[0].Eval(0, depth[n], 1);
984  initT[n] = myTemp;
985  T[n] = myTemp;
986  }
987 
988  DepthP.Eval_All(0, depth[n], data, cdata);
989  OCP_DBL Po = data[1];
990  OCP_DBL Pg = data[2];
991  OCP_DBL Pw = data[3];
992  OCP_DBL Pcgo = Pg - Po;
993  OCP_DBL Pcow = Po - Pw;
994  OCP_DBL Sw = flow[SATNUM[n]]->GetSwByPcow(Pcow);
995  OCP_DBL Sg = 0;
996  if (gas) {
997  Sg = flow[SATNUM[n]]->GetSgByPcgo(Pcgo);
998  }
999  if (Sw + Sg > 1) {
1000  // should me modified
1001  OCP_DBL Pcgw = Pcow + Pcgo;
1002  Sw = flow[SATNUM[n]]->GetSwByPcgw(Pcgw);
1003  Sg = 1 - Sw;
1004  }
1005 
1006  if (1 - Sw < TINY) {
1007  // all water
1008  Po = Pw + flow[SATNUM[n]]->GetPcowBySw(1.0);
1009  } else if (1 - Sg < TINY) {
1010  // all gas
1011  Po = Pg - flow[SATNUM[n]]->GetPcgoBySg(1.0);
1012  } else if (1 - Sw - Sg < TINY) {
1013  // water and gas
1014  Po = Pg - flow[SATNUM[n]]->GetPcgoBySg(Sg);
1015  }
1016  P[n] = Po;
1017 
1018  if (depth[n] < DOGC) {
1019  Pbb = Po;
1020  } else if (PBVD_flag) {
1021  Pbb = EQUIL.PBVD.Eval(0, depth[n], 1);
1022  }
1023  Pb[n] = Pbb;
1024 
1025  // cal Sw
1026  OCP_DBL swco = flow[SATNUM[n]]->GetSwco();
1027  if (!FlagPcow[SATNUM[n]]) {
1028  S[n * numPhase + numPhase - 1] = swco;
1029  continue;
1030  }
1031 
1032  Sw = 0;
1033  Sg = 0;
1034  USI ncut = 10;
1035  OCP_DBL avePcow = 0;
1036 
1037  for (USI k = 0; k < ncut; k++) {
1038  OCP_DBL tmpSw = 0;
1039  OCP_DBL tmpSg = 0;
1040  OCP_DBL dep = depth[n] + dz[n] / ncut * (k - (ncut - 1) / 2.0);
1041  DepthP.Eval_All(0, dep, data, cdata);
1042  Po = data[1];
1043  Pg = data[2];
1044  Pw = data[3];
1045  Pcow = Po - Pw;
1046  Pcgo = Pg - Po;
1047  avePcow += Pcow;
1048  tmpSw = flow[SATNUM[n]]->GetSwByPcow(Pcow);
1049  if (gas) {
1050  tmpSg = flow[SATNUM[n]]->GetSgByPcgo(Pcgo);
1051  }
1052  if (tmpSw + tmpSg > 1) {
1053  // should be modified
1054  OCP_DBL Pcgw = Pcow + Pcgo;
1055  tmpSw = flow[SATNUM[n]]->GetSwByPcgw(Pcgw);
1056  tmpSg = 1 - tmpSw;
1057  }
1058  Sw += tmpSw;
1059  // Sg += tmpSg;
1060  }
1061  Sw /= ncut;
1062  // Sg /= ncut;
1063  avePcow /= ncut;
1064 
1065  flow[SATNUM[n]]->SetupScale(n, Sw, avePcow);
1066  S[n * numPhase + numPhase - 1] = Sw;
1067  }
1068 }
1069 
1071 // Optional Features
1073 
1075 // Region
1077 
1078 void Bulk::AllocateRegion(const Grid& myGrid)
1079 {
1080  SATNUM.resize(numBulk, 0);
1081  PVTNUM.resize(numBulk, 0);
1082  ROCKNUM.resize(numBulk, 0);
1083 
1084  // Pass initial grid value
1085  for (OCP_USI bIda = 0; bIda < numBulk; bIda++) {
1086  OCP_USI bId = myGrid.map_Act2All[bIda];
1087 
1088  SATNUM[bIda] = myGrid.SATNUM[bId];
1089  PVTNUM[bIda] = myGrid.PVTNUM[bId];
1090  ROCKNUM[bIda] = myGrid.ROCKNUM[bId];
1091  }
1092 }
1093 
1094 void Bulk::SetupBulkType(const Grid& myGrid)
1095 {
1096  bType.resize(numBulk, 0);
1097 
1098  OCP_USI count = 0;
1099  for (OCP_USI n = 0; n < myGrid.numGrid; n++) {
1100  if (myGrid.map_All2Act[n].IsAct()) {
1101  if (myGrid.map_All2Flu[n].IsAct()) {
1102  bType[count]++;
1103  }
1104  count++;
1105  }
1106  }
1107 }
1108 
1110 // Basic PVT Model Information
1112 
1114 // Basic Grid and Basic Rock Information
1116 
1117 void Bulk::AllocateGridRockIsoT(const Grid& myGrid)
1118 {
1119  dx.resize(numBulk, 0);
1120  dy.resize(numBulk, 0);
1121  dz.resize(numBulk, 0);
1122  v.resize(numBulk, 0);
1123  depth.resize(numBulk, 0);
1124 
1125  ntg.resize(numBulk, 0);
1126  poroInit.resize(numBulk, 0);
1127  rockKx.resize(numBulk, 0);
1128  rockKy.resize(numBulk, 0);
1129  rockKz.resize(numBulk, 0);
1130 
1131  for (OCP_USI bIda = 0; bIda < numBulk; bIda++) {
1132  OCP_USI bId = myGrid.map_Act2All[bIda];
1133 
1134  dx[bIda] = myGrid.dx[bId];
1135  dy[bIda] = myGrid.dy[bId];
1136  dz[bIda] = myGrid.dz[bId];
1137  v[bIda] = myGrid.v[bId];
1138  depth[bIda] = myGrid.depth[bId];
1139 
1140  ntg[bIda] = myGrid.ntg[bId];
1141  poroInit[bIda] = myGrid.poro[bId];
1142  rockKx[bIda] = myGrid.kx[bId];
1143  rockKy[bIda] = myGrid.ky[bId];
1144  rockKz[bIda] = myGrid.kz[bId];
1145  }
1146 }
1147 
1148 void Bulk::AllocateGridRockT(const Grid& myGrid)
1149 {
1150 
1151  dx.resize(numBulk, 0);
1152  dy.resize(numBulk, 0);
1153  dz.resize(numBulk, 0);
1154  v.resize(numBulk, 0);
1155  depth.resize(numBulk, 0);
1156 
1157  ntg.resize(numBulk, 0);
1158  poroInit.resize(numBulk, 0);
1159  rockKx.resize(numBulk, 0);
1160  rockKy.resize(numBulk, 0);
1161  rockKz.resize(numBulk, 0);
1162  thconr.resize(numBulk, 0);
1163 
1164  bLocation.resize(numBulk, 0);
1165 
1166  for (OCP_USI bIda = 0; bIda < numBulk; bIda++) {
1167  OCP_USI bId = myGrid.map_Act2All[bIda];
1168 
1169  dx[bIda] = myGrid.dx[bId];
1170  dy[bIda] = myGrid.dy[bId];
1171  dz[bIda] = myGrid.dz[bId];
1172  v[bIda] = myGrid.v[bId];
1173  depth[bIda] = myGrid.depth[bId];
1174 
1175  ntg[bIda] = myGrid.ntg[bId];
1176  poroInit[bIda] = myGrid.poro[bId];
1177  rockKx[bIda] = myGrid.kx[bId];
1178  rockKy[bIda] = myGrid.ky[bId];
1179  rockKz[bIda] = myGrid.kz[bId];
1180  thconr[bIda] = myGrid.thconr[bId];
1181 
1182  bLocation[bIda] = myGrid.gLocation[bId];
1183  }
1184 }
1185 
1187 // Basic Fluid Information
1189 
1190 OCP_DBL Bulk::CalFPR() const
1191 {
1192  OCP_FUNCNAME;
1193 
1194  OCP_DBL ptmp = 0;
1195  OCP_DBL vtmp = 0;
1196  OCP_DBL tmp = 0;
1197 
1198  if (numPhase == 3) {
1199  for (OCP_USI n = 0; n < numBulk; n++) {
1200  tmp = rockVp[n] * (1 - S[n * numPhase + 2]);
1201  ptmp += P[n] * tmp;
1202  vtmp += tmp;
1203  }
1204  } else if (numPhase < 3) {
1205  for (OCP_USI n = 0; n < numBulk; n++) {
1206  tmp = rockVp[n] * (S[n * numPhase]);
1207  ptmp += P[n] * tmp;
1208  vtmp += tmp;
1209  }
1210  } else {
1211  OCP_ABORT("Number of phases is out of range!");
1212  }
1213  return ptmp / vtmp;
1214 }
1215 
1216 OCP_DBL Bulk::CalFTR() const
1217 {
1218  OCP_FUNCNAME;
1219 
1220  OCP_DBL Ttmp = 0;
1221  OCP_DBL vtmp = 0;
1222 
1223  for (OCP_USI n = 0; n < numBulk; n++) {
1224  Ttmp += T[n] * v[n];
1225  vtmp += v[n];
1226  }
1227 
1228  return Ttmp / vtmp;
1229 }
1230 
1232 // Important Indicator Variable and Check
1234 
1235 OCP_DBL Bulk::CalNRdSmax(OCP_USI& index)
1236 {
1237  NRdSmax = 0;
1238  OCP_USI len = numBulk * numPhase;
1239  for (USI n = 0; n < len; n++) {
1240  if (fabs(NRdSmax) < fabs(dSNR[n])) {
1241  NRdSmax = dSNR[n];
1242  index = n;
1243  }
1244  }
1245  index /= numPhase;
1246  return NRdSmax;
1247 }
1248 
1250 OCP_INT Bulk::CheckP() const
1251 {
1252  OCP_FUNCNAME;
1253 
1254  for (OCP_USI n = 0; n < numBulk; n++) {
1255  if (P[n] < 0) {
1256  std::ostringstream PStringSci;
1257  PStringSci << std::scientific << P[n];
1258  OCP_WARNING("Negative pressure: P[" + std::to_string(n) +
1259  "] = " + PStringSci.str());
1260  cout << "P = " << P[n] << endl;
1261  return BULK_NEGATIVE_PRESSURE;
1262  }
1263  }
1264 
1265  return BULK_SUCCESS;
1266 }
1267 
1268 OCP_INT Bulk::CheckT() const
1269 {
1270  for (OCP_USI n = 0; n < numBulk; n++) {
1271  if (T[n] < 0) {
1272  std::ostringstream PStringSci;
1273  PStringSci << std::scientific << T[n];
1274  OCP_WARNING("Negative pressure: T[" + std::to_string(n) +
1275  "] = " + PStringSci.str());
1276  cout << "T = " << T[n] << endl;
1277  return BULK_NEGATIVE_TEMPERATURE;
1278  }
1279  }
1280 
1281  return BULK_SUCCESS;
1282 }
1283 
1285 OCP_INT Bulk::CheckNi()
1286 {
1287  OCP_FUNCNAME;
1288 
1289  OCP_USI len = numBulk * numCom;
1290  for (OCP_USI n = 0; n < len; n++) {
1291  if (Ni[n] < 0) {
1292  OCP_USI bId = n / numCom;
1293  if (Ni[n] > -1E-3 * Nt[bId] && OCP_FALSE) {
1294  Ni[n] = 1E-8 * Nt[bId];
1295  } else {
1296  USI cId = n - bId * numCom;
1297  std::ostringstream NiStringSci;
1298  NiStringSci << std::scientific << Ni[n];
1299  OCP_WARNING("Negative Ni: Ni[" + std::to_string(cId) + "] in Bulk[" +
1300  std::to_string(bId) + "] = " + NiStringSci.str() + ", " +
1301  "dNi = " + std::to_string(dNNR[n]));
1302 
1303  return BULK_NEGATIVE_COMPONENTS_MOLES;
1304  }
1305  }
1306  }
1307  return BULK_SUCCESS;
1308 }
1309 
1311 OCP_INT Bulk::CheckVe(const OCP_DBL& Vlim) const
1312 {
1313  OCP_FUNCNAME;
1314 
1315  OCP_DBL dVe = 0.0;
1316  for (OCP_USI n = 0; n < numBulk; n++) {
1317  dVe = fabs(vf[n] - rockVp[n]) / rockVp[n];
1318  if (dVe > Vlim) {
1319  cout << "Volume error at Bulk[" << n << "] = " << setprecision(6) << dVe
1320  << " is too big!" << endl;
1321  return BULK_OUTRANGED_VOLUME_ERROR;
1322  }
1323  }
1324  return BULK_SUCCESS;
1325 }
1326 
1327 OCP_INT Bulk::CheckCFL(const OCP_DBL& cflLim) const
1328 {
1329  if (maxCFL > cflLim)
1330  return BULK_OUTRANGED_CFL;
1331  else
1332  return BULK_SUCCESS;
1333 }
1334 
1335 void Bulk::CalMaxChange()
1336 {
1337  OCP_FUNCNAME;
1338 
1339  dPmax = 0;
1340  dTmax = 0;
1341  dNmax = 0;
1342  dSmax = 0;
1343  eVmax = 0;
1344  OCP_DBL tmp = 0;
1345  OCP_USI id;
1346 
1347  for (OCP_USI n = 0; n < numBulk; n++) {
1348 
1349  // dP
1350  tmp = fabs(P[n] - lP[n]);
1351  if (dPmax < tmp) {
1352  dPmax = tmp;
1353  }
1354 
1355  // dT
1356  tmp = fabs(T[n] - lT[n]);
1357  if (dTmax < tmp) {
1358  dTmax = tmp;
1359  }
1360 
1361  // dS
1362  for (USI j = 0; j < numPhase; j++) {
1363  id = n * numPhase + j;
1364  tmp = fabs(S[id] - lS[id]);
1365  if (dSmax < tmp) {
1366  dSmax = tmp;
1367  }
1368  }
1369 
1370  // dN
1371  for (USI i = 0; i < numCom; i++) {
1372  id = n * numCom + i;
1373  tmp = fabs(max(Ni[id], lNi[id]));
1374  if (tmp > TINY) {
1375  tmp = fabs(Ni[id] - lNi[id]) / tmp;
1376  if (dNmax < tmp) {
1377  dNmax = tmp;
1378  }
1379  }
1380  }
1381 
1382  // Ve
1383  tmp = fabs(vf[n] - rockVp[n]) / rockVp[n];
1384  if (eVmax < tmp) {
1385  eVmax = tmp;
1386  }
1387  }
1388 }
1389 
1391 // Error
1393 
1394 void Bulk::AllocateError()
1395 {
1396  if (ifUseEoS) {
1397  ePEC.resize(numBulk);
1398  }
1399 }
1400 
1402 // For AIMc
1404 
1405 void Bulk::ShowFIMBulk(const OCP_BOOL& flag) const
1406 {
1407 
1408  if (flag) {
1409  USI iter = 0;
1410  for (USI n = 0; n < numBulk; n++) {
1411  if (bulkTypeAIM.IfFIMbulk(n)) {
1412  // FIM bulk
1413  cout << setw(6) << n << " ";
1414  iter++;
1415  }
1416 
1417  if ((iter + 1) % 10 == 0) {
1418  cout << endl;
1419  iter = 0;
1420  }
1421  }
1422  cout << endl;
1423  }
1424 }
1425 
1426 /*----------------------------------------------------------------------------*/
1427 /* Brief Change History of This File */
1428 /*----------------------------------------------------------------------------*/
1429 /* Author Date Actions */
1430 /*----------------------------------------------------------------------------*/
1431 /* Shizhe Li Oct/01/2021 Create file */
1432 /* Chensong Zhang Oct/15/2021 Format file */
1433 /* Chensong Zhang Jan/09/2022 Update Doxygen */
1434 /*----------------------------------------------------------------------------*/
Bulk.hpp
Bulk class declaration.
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:38
PHASE_ODGW01
const USI PHASE_ODGW01
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:117
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:23
GAS
const USI GAS
Fluid type = gas.
Definition: OCPConst.hpp:92
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:27
GRAVITY_FACTOR
const OCP_DBL GRAVITY_FACTOR
0.00694444 ft2 psi / lb
Definition: OCPConst.hpp:56
BLKOIL_DOGW
const USI BLKOIL_DOGW
black oil model with dead oil, dry gas, water
Definition: OCPConst.hpp:101
WATER
const USI WATER
Fluid type = water.
Definition: OCPConst.hpp:93
PHASE_ODGW
const USI PHASE_ODGW
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:115
PHASE_ODGW02
const USI PHASE_ODGW02
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:118
BLKOIL_OW
const USI BLKOIL_OW
black oil model with oil and water
Definition: OCPConst.hpp:99
OIL
const USI OIL
Fluid type = oil.
Definition: OCPConst.hpp:91
BLKOIL_OG
const USI BLKOIL_OG
black oil model with oil and gas
Definition: OCPConst.hpp:100
BLKOIL_ODGW
const USI BLKOIL_ODGW
black oil model with live oil, dry gas, water
Definition: OCPConst.hpp:102
PHASE_W
const USI PHASE_W
Phase type = water only.
Definition: OCPConst.hpp:111
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:25
PHASE_ODGW01_MISCIBLE
const USI PHASE_ODGW01_MISCIBLE
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:119
PHASE_OW
const USI PHASE_OW
Phase type = oil-water.
Definition: OCPConst.hpp:113
PHASE_DOGW
const USI PHASE_DOGW
Phase type = dead oil-gas-water.
Definition: OCPConst.hpp:116
OCP_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:26
OCP_BOOL
unsigned int OCP_BOOL
OCP_BOOL in OCP.
Definition: OCPConst.hpp:29
BLKOIL_W
const USI BLKOIL_W
black oil model only with water
Definition: OCPConst.hpp:98
OCP_FUNCNAME
#define OCP_FUNCNAME
Print Function Name.
Definition: UtilError.hpp:73
OCP_WARNING
#define OCP_WARNING(msg)
Log warning messages.
Definition: UtilError.hpp:39
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47
BOMixture_ODGW
Definition: MixtureBO.hpp:257
BOMixture_OW
Definition: MixtureBO.hpp:183
Bulk::SetupOptionalFeatures
void SetupOptionalFeatures(const Grid &myGrid, OptionalFeatures &optFeatures)
Setup optional features.
Definition: Bulk.cpp:429
Bulk::dNNR
vector< OCP_DBL > dNNR
Ni change between NR steps.
Definition: Bulk.hpp:411
Bulk::lP
vector< OCP_DBL > lP
last P
Definition: Bulk.hpp:332
Bulk::lS
vector< OCP_DBL > lS
last S
Definition: Bulk.hpp:336
Bulk::rock
vector< Rock * > rock
rock model
Definition: Bulk.hpp:203
Bulk::flow
vector< FlowUnit * > flow
Vector for capillary pressure, relative perm.
Definition: Bulk.hpp:197
Bulk::NTSFUN
USI NTSFUN
num of SAT regions
Definition: Bulk.hpp:194
Bulk::numPhase
USI numPhase
Number of phase.
Definition: Bulk.hpp:154
Bulk::T
vector< OCP_DBL > T
Temperature: numBulk.
Definition: Bulk.hpp:308
Bulk::eVmax
OCP_DBL eVmax
Definition: Bulk.hpp:462
Bulk::S
vector< OCP_DBL > S
Saturation of phase: numPhase*numBulk.
Definition: Bulk.hpp:314
Bulk::oil
OCP_BOOL oil
If OCP_TRUE, oil phase could exist.
Definition: Bulk.hpp:222
Bulk::bLocation
vector< USI > bLocation
Location of bulk: top, bottom, side.
Definition: Bulk.hpp:206
Bulk::ifComps
OCP_BOOL ifComps
If OCP_TRUE, compositional model will be used.
Definition: Bulk.hpp:219
Bulk::dSmax
OCP_DBL dSmax
Max change in saturation during the current time step.
Definition: Bulk.hpp:460
Bulk::SetupBulkType
void SetupBulkType(const Grid &myGrid)
Setup Bulk type.
Definition: Bulk.cpp:1094
Bulk::ntg
vector< OCP_DBL > ntg
net to gross of bulk.
Definition: Bulk.hpp:243
Bulk::dy
vector< OCP_DBL > dy
Size of cell in y-direction: activeGridNum.
Definition: Bulk.hpp:238
Bulk::PVTNUM
vector< USI > PVTNUM
Identify PVT region in black-oil model: numBulk.
Definition: Bulk.hpp:191
Bulk::gas
OCP_BOOL gas
If OCP_TRUE, gas phase could exist.
Definition: Bulk.hpp:223
Bulk::maxCFL
OCP_DBL maxCFL
max CFL number
Definition: Bulk.hpp:466
Bulk::NTROCC
USI NTROCC
num of Rock regions
Definition: Bulk.hpp:201
Bulk::SATNUM
vector< USI > SATNUM
Identify SAT region: numBulk.
Definition: Bulk.hpp:196
Bulk::depth
vector< OCP_DBL > depth
Depth of center of grid cells: activeGridNum.
Definition: Bulk.hpp:241
Bulk::rockVp
vector< OCP_DBL > rockVp
pore volume = Vgrid * ntg * poro.
Definition: Bulk.hpp:246
Bulk::Ni
vector< OCP_DBL > Ni
Moles of component: numCom*numBulk.
Definition: Bulk.hpp:306
Bulk::SATmode
USI SATmode
Identify SAT mode.
Definition: Bulk.hpp:195
Bulk::CalMaxChange
void CalMaxChange()
Calculate max change of indicator variables.
Definition: Bulk.cpp:1335
Bulk::ePEC
vector< OCP_DBL > ePEC
error for fugacity balance equations, EoS only now
Definition: Bulk.hpp:476
Bulk::NTPVT
USI NTPVT
num of PVT regions
Definition: Bulk.hpp:189
Bulk::dSNR
vector< OCP_DBL > dSNR
saturation change between NR steps
Definition: Bulk.hpp:406
Bulk::hLoss
HeatLoss hLoss
Heat loss iterm.
Definition: Bulk.hpp:207
Bulk::ifThermal
OCP_BOOL ifThermal
Id OCP_TRUE, ifThermal model will be used.
Definition: Bulk.hpp:220
Bulk::CheckNi
OCP_INT CheckNi()
Check if negative Ni occurs.
Definition: Bulk.cpp:1285
Bulk::Nt
vector< OCP_DBL > Nt
Total moles of components in bulks: numBulk.
Definition: Bulk.hpp:305
Bulk::dz
vector< OCP_DBL > dz
Size of cell in z-direction: activeGridNum.
Definition: Bulk.hpp:239
Bulk::CheckT
OCP_INT CheckT() const
Check if negative T occurs.
Definition: Bulk.cpp:1268
Bulk::thconr
vector< OCP_DBL > thconr
Rock ifThermal conductivity: activeGridNum.
Definition: Bulk.hpp:250
Bulk::AllocateRegion
void AllocateRegion(const Grid &myGrid)
Allocate memory for region num.
Definition: Bulk.cpp:1078
Bulk::lT
vector< OCP_DBL > lT
last T
Definition: Bulk.hpp:331
Bulk::flashCal
vector< Mixture * > flashCal
Flash calculation class.
Definition: Bulk.hpp:192
Bulk::InputParam
void InputParam(const ParamReservoir &rs_param)
Input param from internal data structure ParamReservoir.
Definition: Bulk.cpp:99
Bulk::phase2Index
vector< USI > phase2Index
Location of phase according to its name: numPhase.
Definition: Bulk.hpp:302
Bulk::initZi_Tab
vector< OCPTable > initZi_Tab
Initial mole ratio of components vs. depth, table set.
Definition: Bulk.hpp:167
Bulk::satcm
vector< vector< OCP_DBL > > satcm
critical saturation when phase becomes mobile / immobile.
Definition: Bulk.hpp:199
Bulk::CalNRdSmax
OCP_DBL CalNRdSmax(OCP_USI &index)
Calculate some auxiliary variable, for example, dSmax.
Definition: Bulk.cpp:1235
Bulk::vf
vector< OCP_DBL > vf
Total fluid volume: numBulk.
Definition: Bulk.hpp:307
Bulk::rockKz
vector< OCP_DBL > rockKz
current rock permeability along the z direction.
Definition: Bulk.hpp:249
Bulk::CalFTR
OCP_DBL CalFTR() const
Calculate average Temperature in reservoir.
Definition: Bulk.cpp:1216
Bulk::EQUIL
ParamEQUIL EQUIL
Initial Equilibration.
Definition: Bulk.hpp:170
Bulk::ifUseEoS
OCP_BOOL ifUseEoS
If OCP_TRUE, then EoS model is used.
Definition: Bulk.hpp:221
Bulk::initT_Tab
vector< OCPTable > initT_Tab
Initial temperature vs. depth, table set.
Definition: Bulk.hpp:168
Bulk::dNmax
OCP_DBL dNmax
Max change in moles of component during the current time step.
Definition: Bulk.hpp:461
Bulk::rockKx
vector< OCP_DBL > rockKx
current rock permeability along the x direction.
Definition: Bulk.hpp:247
Bulk::NRdSmax
OCP_DBL NRdSmax
Max saturation difference in an NR step(Real)
Definition: Bulk.hpp:420
Bulk::rsTemp
OCP_DBL rsTemp
Reservoir temperature.
Definition: Bulk.hpp:171
Bulk::bType
vector< USI > bType
Indicate bulk type, 0: rock, 1: rock and fluid.
Definition: Bulk.hpp:205
Bulk::ShowFIMBulk
void ShowFIMBulk(const OCP_BOOL &flag=OCP_FALSE) const
Print Bulk which are implicit.
Definition: Bulk.cpp:1405
Bulk::SetupIsoT
void SetupIsoT(const Grid &myGrid)
Allocate memory for fluid grid for isothermal model.
Definition: Bulk.cpp:408
Bulk::PVTmodeB
USI PVTmodeB
Identify PVT mode in black-oil model.
Definition: Bulk.hpp:190
Bulk::rockKy
vector< OCP_DBL > rockKy
current rock permeability along the y direction.
Definition: Bulk.hpp:248
Bulk::disGas
OCP_BOOL disGas
If OCP_TRUE, dissolve gas in live oil could exist.
Definition: Bulk.hpp:225
Bulk::dTmax
OCP_DBL dTmax
Max change in temperature during the current time step.
Definition: Bulk.hpp:459
Bulk::CheckP
OCP_INT CheckP() const
Check if negative P occurs.
Definition: Bulk.cpp:1250
Bulk::water
OCP_BOOL water
If OCP_TRUE, water phase could exist.
Definition: Bulk.hpp:224
Bulk::Pb
vector< OCP_DBL > Pb
Bubble point pressure: numBulk.
Definition: Bulk.hpp:310
Bulk::SetupT
void SetupT(const Grid &myGrid)
Allocate memory for fluid grid for ifThermal model.
Definition: Bulk.cpp:419
Bulk::v
vector< OCP_DBL > v
Volume of grids: activeGridNum.
Definition: Bulk.hpp:240
Bulk::numComH
USI numComH
Number of HydroCarbon.
Definition: Bulk.hpp:156
Bulk::InitPTSw
void InitPTSw(const USI &tabrow)
Calculate initial equilibrium – hydrostatic equilibration.
Definition: Bulk.cpp:443
Bulk::ifBlackOil
OCP_BOOL ifBlackOil
If OCP_TRUE, black-oil model will be used.
Definition: Bulk.hpp:218
Bulk::CalFPR
OCP_DBL CalFPR() const
Calculate average pressure in reservoir.
Definition: Bulk.cpp:1190
Bulk::numBulk
OCP_USI numBulk
Number of bulks (active grids).
Definition: Bulk.hpp:153
Bulk::CheckVe
OCP_INT CheckVe(const OCP_DBL &Vlim) const
Check if relative volume error is outranged.
Definition: Bulk.cpp:1311
Bulk::poroInit
vector< OCP_DBL > poroInit
initial rock porosity * ntg.
Definition: Bulk.hpp:244
Bulk::thconp
vector< OCP_DBL > thconp
Phase thermal conductivity: numPhase.
Definition: Bulk.hpp:172
Bulk::ROCKNUM
vector< USI > ROCKNUM
index of Rock table for each bulk
Definition: Bulk.hpp:202
Bulk::dx
vector< OCP_DBL > dx
Size of cell in x-direction: activeGridNum.
Definition: Bulk.hpp:237
Bulk::initT
vector< OCP_DBL > initT
Initial temperature of each bulk: numBulk.
Definition: Bulk.hpp:169
Bulk::CheckCFL
OCP_INT CheckCFL(const OCP_DBL &cflLim) const
Check if Cfl is outranged.
Definition: Bulk.cpp:1327
Bulk::P
vector< OCP_DBL > P
Pressure: numBulk.
Definition: Bulk.hpp:309
Bulk::lNi
vector< OCP_DBL > lNi
last Ni
Definition: Bulk.hpp:329
Bulk::AllocateGridRockIsoT
void AllocateGridRockIsoT(const Grid &myGrid)
Allocate memory for Rock properties.
Definition: Bulk.cpp:1117
Bulk::numCom
USI numCom
Number of component.
Definition: Bulk.hpp:155
Bulk::dPmax
OCP_DBL dPmax
Max change in pressure during the current time step.
Definition: Bulk.hpp:458
ComponentParam::numPhase
USI numPhase
num of phase, water is excluded, constant now.
Definition: ParamReservoir.hpp:120
ComponentParam::numCom
USI numCom
num of components, water is excluded.
Definition: ParamReservoir.hpp:119
FlowUnit_ODGW01_Miscible
Definition: FlowUnit.hpp:283
FlowUnit_ODGW01
Definition: FlowUnit.hpp:218
FlowUnit_ODGW02
Definition: FlowUnit.hpp:337
FlowUnit_OW
Definition: FlowUnit.hpp:108
FlowUnit_W
Definition: FlowUnit.hpp:72
Grid
Definition: Grid.hpp:89
Grid::ky
vector< OCP_DBL > ky
Absolute permeability in y-direction: numGrid.
Definition: Grid.hpp:172
Grid::v
vector< OCP_DBL > v
Volume of cells: numGrid.
Definition: Grid.hpp:165
Grid::map_Act2All
vector< OCP_USI > map_Act2All
Mapping from active grid to all grid: activeGridNum.
Definition: Grid.hpp:192
Grid::ROCKNUM
vector< USI > ROCKNUM
index of rock table for each grid: numGrid
Definition: Grid.hpp:181
Grid::activeGridNum
OCP_USI activeGridNum
Num of active grid.
Definition: Grid.hpp:190
Grid::PVTNUM
vector< USI > PVTNUM
Identify PVT region for the blackoil model: numGrid.
Definition: Grid.hpp:178
Grid::depth
vector< OCP_DBL > depth
Depth of center of grid cells: numGrid.
Definition: Grid.hpp:166
Grid::SATNUM
vector< USI > SATNUM
Identify SAT region: numGrid.
Definition: Grid.hpp:177
Grid::kz
vector< OCP_DBL > kz
Absolute permeability in z-direction: numGrid.
Definition: Grid.hpp:173
Grid::gLocation
vector< USI > gLocation
Top face, bottom face, side face, numGrid.
Definition: Grid.hpp:156
Grid::dz
vector< OCP_DBL > dz
Size of cell in z-direction: numGrid.
Definition: Grid.hpp:164
Grid::map_All2Flu
vector< GB_Pair > map_All2Flu
Mapping from all grid to fluid grid: numGrid.
Definition: Grid.hpp:196
Grid::thconr
vector< OCP_DBL > thconr
Rock if Thermal conductivity: numGrid.
Definition: Grid.hpp:174
Grid::poro
vector< OCP_DBL > poro
Initial porosity of rock cells: numGrid.
Definition: Grid.hpp:170
Grid::dy
vector< OCP_DBL > dy
Size of cell in y-direction: numGrid.
Definition: Grid.hpp:163
Grid::map_All2Act
vector< GB_Pair > map_All2Act
Mapping from grid to active all grid: numGrid.
Definition: Grid.hpp:193
Grid::dx
vector< OCP_DBL > dx
Size of cell in x-direction: numGrid.
Definition: Grid.hpp:162
Grid::numGrid
OCP_USI numGrid
Number of all cells.
Definition: Grid.hpp:142
Grid::ntg
vector< OCP_DBL > ntg
Net to gross ratio of cells: numGrid.
Definition: Grid.hpp:169
Grid::kx
vector< OCP_DBL > kx
Absolute permeability in x-direction: numGrid.
Definition: Grid.hpp:171
HLoss::obC
OCP_DBL obC
Volumetric heat capacity of overburden rock.
Definition: ParamReservoir.hpp:53
HLoss::ifHLoss
OCP_BOOL ifHLoss
If use Heat loss.
Definition: ParamReservoir.hpp:52
HLoss::obK
OCP_DBL obK
Thermal conductivity of overburden rock.
Definition: ParamReservoir.hpp:54
HLoss::ubK
OCP_DBL ubK
Thermal conductivity of underburden rock.
Definition: ParamReservoir.hpp:56
HLoss::ubC
OCP_DBL ubC
Volumetric heat capacity of underburden rock.
Definition: ParamReservoir.hpp:55
HeatLoss::pT
vector< OCP_DBL > pT
dP / dT
Definition: Bulk.hpp:90
HeatLoss::ubD
OCP_DBL ubD
Thermal diffusivity of underburden rock.
Definition: Bulk.hpp:87
HeatLoss::ubK
OCP_DBL ubK
Thermal conductivity of underburden rock.
Definition: Bulk.hpp:83
HeatLoss::p
vector< OCP_DBL > p
Auxiliary variable.
Definition: Bulk.hpp:89
HeatLoss::obC
OCP_DBL obC
Volumetric heat capacity of overburden rock.
Definition: Bulk.hpp:80
HeatLoss::obD
OCP_DBL obD
Thermal diffusivity of overburden rock.
Definition: Bulk.hpp:86
HeatLoss::lpT
vector< OCP_DBL > lpT
last pT
Definition: Bulk.hpp:94
HeatLoss::obK
OCP_DBL obK
Thermal conductivity of overburden rock.
Definition: Bulk.hpp:81
HeatLoss::I
vector< OCP_DBL > I
Auxiliary variable.
Definition: Bulk.hpp:88
HeatLoss::lp
vector< OCP_DBL > lp
Auxiliary variable.
Definition: Bulk.hpp:93
HeatLoss::ifHLoss
OCP_BOOL ifHLoss
If use Heat loss.
Definition: Bulk.hpp:79
HeatLoss::lI
vector< OCP_DBL > lI
Auxiliary variable.
Definition: Bulk.hpp:92
HeatLoss::numBulk
OCP_USI numBulk
Num of Bulk.
Definition: Bulk.hpp:85
HeatLoss::ubC
OCP_DBL ubC
Volumetric heat capacity of underburden rock.
Definition: Bulk.hpp:82
MixtureThermal_K01
Definition: MixtureThermal.hpp:66
OCPTable::Eval_All
USI Eval_All(const USI &j, const OCP_DBL &val, vector< OCP_DBL > &outdata, vector< OCP_DBL > &slope)
Careful: the memory outdata and slope have not be allocated before.
Definition: OCPTable.cpp:47
OCPTable::Display
void Display() const
Display the data of table on screen.
Definition: OCPTable.cpp:205
OCPTable::Setup
void Setup(const vector< vector< OCP_DBL >> &src)
Setup tables from existing data of table.
Definition: OCPTable.cpp:27
OCPTable::GetCol
vector< OCP_DBL > & GetCol(const USI &j)
return the jth column in table to modify or use.
Definition: OCPTable.hpp:55
OCPTable::IsEmpty
OCP_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:126
ParamReservoir::rockSet
vector< RockParam > rockSet
a set of rock
Definition: ParamReservoir.hpp:234
ParamReservoir::oil
OCP_BOOL oil
If OCP_TRUE, oil phase could exist.
Definition: ParamReservoir.hpp:255
ParamReservoir::NTSFUN
USI NTSFUN
Num of SAT regions.
Definition: ParamReservoir.hpp:263
ParamReservoir::NTROOC
USI NTROOC
Num of Rock regions.
Definition: ParamReservoir.hpp:265
ParamReservoir::ZMFVD_T
TableSet ZMFVD_T
Table set of ZMFVD.
Definition: ParamReservoir.hpp:279
ParamReservoir::EQUIL
vector< OCP_DBL > EQUIL
See ParamEQUIL.
Definition: ParamReservoir.hpp:281
ParamReservoir::comsParam
ComponentParam comsParam
information for components
Definition: ParamReservoir.hpp:260
ParamReservoir::thconw
OCP_DBL thconw
water ifThermal conductivity
Definition: ParamReservoir.hpp:249
ParamReservoir::gas
OCP_BOOL gas
If OCP_TRUE, gas phase could exist.
Definition: ParamReservoir.hpp:256
ParamReservoir::SOF3_T
TableSet SOF3_T
Table set of SOF3.
Definition: ParamReservoir.hpp:276
ParamReservoir::PBVD_T
TableSet PBVD_T
Table set of PBVD.
Definition: ParamReservoir.hpp:277
ParamReservoir::comps
OCP_BOOL comps
If OCP_TRUE, compositional model will be used.
Definition: ParamReservoir.hpp:253
ParamReservoir::TEMPVD_T
TableSet TEMPVD_T
Table set of TEMPVD.
Definition: ParamReservoir.hpp:280
ParamReservoir::thcong
OCP_DBL thcong
gas ifThermal conductivity
Definition: ParamReservoir.hpp:248
ParamReservoir::blackOil
OCP_BOOL blackOil
If ture, blackoil model will be used.
Definition: ParamReservoir.hpp:252
ParamReservoir::disGas
OCP_BOOL disGas
If OCP_TRUE, dissolve gas could exist in oil phase.
Definition: ParamReservoir.hpp:258
ParamReservoir::hLoss
HLoss hLoss
Heat loss property.
Definition: ParamReservoir.hpp:235
ParamReservoir::thcono
OCP_DBL thcono
oil ifThermal conductivity
Definition: ParamReservoir.hpp:247
ParamReservoir::NTPVT
USI NTPVT
Num of PVT regions.
Definition: ParamReservoir.hpp:264
ParamReservoir::rsTemp
OCP_DBL rsTemp
Temperature for reservoir.
Definition: ParamReservoir.hpp:233
ParamReservoir::water
OCP_BOOL water
If OCP_TRUE, water phase could exist.
Definition: ParamReservoir.hpp:257
ParamReservoir::thermal
OCP_BOOL thermal
If OCP_TRUE, ifThermal model will be used.
Definition: ParamReservoir.hpp:254
RockT_Exp
Definition: Rock.hpp:147
RockT_Linear
Definition: Rock.hpp:129
Rock_Linear
Definition: Rock.hpp:70
TableSet::data
vector< vector< vector< OCP_DBL > > > data
All table with the same name.
Definition: ParamReservoir.hpp:37
Generated on Tue Feb 7 2023 01:49:50 for OpenCAEPoro by doxygen 1.9.1