OpenCAEPoro  0.2.0 Sep/22/2022
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 // General
22 
24 void Bulk::InputParam(ParamReservoir& rs_param)
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 }
181 
183 void Bulk::Setup(const Grid& myGrid)
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 }
341 
342 void Bulk::CheckSetup() const
343 {
344  CheckInitVpore();
345  CheckVpore();
346 }
347 
349 void Bulk::CheckInitVpore() const
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 }
358 
359 // Check pore volume.
360 void Bulk::CheckVpore() const
361 {
362  for (OCP_USI n = 0; n < numBulk; n++) {
363  if (rockVp[n] < TINY) {
364  OCP_ABORT("Bulk volume is too small: bulk[" + std::to_string(n) +
365  "] = " + std::to_string(rockVp[n]));
366  }
367  }
368 }
369 
371 void Bulk::InitSjPcBo(const USI& tabrow)
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 }
810 
812 void Bulk::InitSjPcComp(const USI& tabrow, const Grid& myGrid)
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 }
1276 
1280 void Bulk::InitFlash(const bool& flag)
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 }
1301 
1302 void Bulk::InitFlashDer()
1303 {
1304  OCP_FUNCNAME;
1305 
1306  if (comps) {
1307  for (OCP_USI n = 0; n < numBulk; n++) {
1308  flashCal[PVTNUM[n]]->InitFlashDer(P[n], Pb[n], T, &S[n * numPhase],
1309  rockVp[n], initZi.data() + n * numCom);
1310  for (USI i = 0; i < numCom; i++) {
1311  Ni[n * numCom + i] = flashCal[PVTNUM[n]]->Ni[i];
1312  }
1313  PassFlashValueDeriv(n);
1314  }
1315  } else {
1316  InitFlash(false);
1317  FlashDeriv();
1318  }
1319 }
1320 
1321 void Bulk::InitFlashDer_n()
1322 {
1323  OCP_FUNCNAME;
1324 
1325  if (comps) {
1326  for (OCP_USI n = 0; n < numBulk; n++) {
1327  flashCal[PVTNUM[n]]->InitFlashDer_n(P[n], Pb[n], T, &S[n * numPhase],
1328  rockVp[n], initZi.data() + n * numCom);
1329  for (USI i = 0; i < numCom; i++) {
1330  Ni[n * numCom + i] = flashCal[PVTNUM[n]]->Ni[i];
1331  }
1332  PassFlashValueDeriv_n(n);
1333  }
1334  } else {
1335  OCP_ABORT("Not Completed in BLKOIL MODEL!");
1336  }
1337 }
1338 
1340 void Bulk::Flash()
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 }
1354 
1355 void Bulk::FlashBLKOIL()
1356 {
1357  for (OCP_USI n = 0; n < numBulk; n++) {
1358  flashCal[PVTNUM[n]]->Flash(P[n], T, &Ni[n * numCom], 0, 0, 0);
1359  PassFlashValue(n);
1360  }
1361 }
1362 
1363 void Bulk::FlashCOMP()
1364 {
1365  USI ftype;
1366  OCP_USI bId;
1367  OCP_DBL Ntw;
1368  OCP_DBL minEig;
1369  // cout << endl << "==================================" << endl;
1370  for (OCP_USI n = 0; n < numBulk; n++) {
1371  ftype = 1;
1372  if (flagSkip[n]) {
1373  minEig = minEigenSkip[n];
1374  if (fabs(1 - PSkip[n] / P[n]) >= minEig / 10) {
1375  ftype = 0;
1376  }
1377  // cout << setprecision(2) << scientific << minEig / 10 << " " << fabs(1 -
1378  // lP[n] / P[n]) << " ";
1379  if (ftype == 1) {
1380  bId = n * numCom;
1381  Ntw = Nt[n] - Ni[bId + numCom - 1];
1382  for (USI i = 0; i < numCom - 1; i++) {
1383  // cout << fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) << " ";
1384  if (fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) >= minEig / 10) {
1385  ftype = 0;
1386  break;
1387  }
1388  }
1389  }
1390  // cout << n << endl;
1391  } else {
1392  ftype = 0;
1393  }
1394  flashCal[PVTNUM[n]]->Flash(P[n], T, &Ni[n * numCom], ftype, phaseNum[n],
1395  &Ks[n * numCom_1]);
1396  PassFlashValue(n);
1397  }
1398 }
1399 
1401 void Bulk::FlashDeriv()
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 }
1415 
1416 void Bulk::FlashDeriv_n()
1417 {
1418  OCP_FUNCNAME;
1419 
1420  if (comps) {
1421  FlashDerivCOMP_n();
1422  } else {
1423  FlashDerivBLKOIL_n();
1424  }
1425 }
1426 
1428 void Bulk::FlashDerivBLKOIL()
1429 {
1430  // dSec_dPri.clear();
1431  for (OCP_USI n = 0; n < numBulk; n++) {
1432  flashCal[PVTNUM[n]]->FlashDeriv(P[n], T, &Ni[n * numCom], 0, 0, 0);
1433  PassFlashValueDeriv(n);
1434  }
1435 }
1436 
1437 void Bulk::FlashDerivBLKOIL_n() {}
1438 
1440 void Bulk::FlashDerivCOMP()
1441 {
1442  USI ftype;
1443  NRdSSP = 0;
1444  maxNRdSSP = 0;
1445  index_maxNRdSSP = 0;
1446 
1447  for (OCP_USI n = 0; n < numBulk; n++) {
1448 
1449  ftype = CalFlashType(n);
1450 
1451  flashCal[PVTNUM[n]]->FlashDeriv(P[n], T, &Ni[n * numCom], ftype, phaseNum[n],
1452  &Ks[n * numCom_1]);
1453  PassFlashValueDeriv(n);
1454  }
1455 }
1456 
1457 void Bulk::FlashDerivCOMP_n()
1458 {
1459  USI ftype;
1460  vector<USI> flagB(numPhase, 0);
1461  NRdSSP = 0;
1462  maxNRdSSP = 0;
1463  index_maxNRdSSP = 0;
1464 
1465  for (OCP_USI n = 0; n < numBulk; n++) {
1466 
1467  ftype = CalFlashType(n);
1468 
1469  for (USI j = 0; j < numPhase; j++) flagB[j] = phaseExist[n * numPhase + j];
1470 
1471  flashCal[PVTNUM[n]]->FlashDeriv_n(
1472  P[n], T, &Ni[n * numCom], &S[n * numPhase], &xij[n * numPhase * numCom],
1473  &nj[n * numPhase], ftype, &flagB[0], phaseNum[n], &Ks[n * numCom_1]);
1474 
1475  PassFlashValueDeriv_n(n);
1476  }
1477 }
1478 
1479 
1480 USI Bulk::CalFlashType(const OCP_USI& n) const
1481 {
1482  USI ftype = 1;
1483  OCP_USI bId;
1484  OCP_DBL Ntw;
1485  OCP_DBL minEig;
1486 
1487  if (flagSkip[n]) {
1488  // If only single hydrocarbon phase existed at last NR step and flagSkip
1489  // is available, then test if phase stability could be skipped.
1490  minEig = minEigenSkip[n];
1491  if (fabs(1 - PSkip[n] / P[n]) >= minEig / 10) {
1492  ftype = 0;
1493  }
1494  if (ftype == 1) {
1495  bId = n * numCom;
1496  Ntw = Nt[n] - Ni[bId + numCom_1];
1497  for (USI i = 0; i < numCom_1; i++) {
1498  if (fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) >= minEig / 10) {
1499  ftype = 0;
1500  break;
1501  }
1502  }
1503  }
1504  }
1505  else {
1506  ftype = 2;
1507  if (phaseNum[n] == 2 && true) {
1508  // if num of hydrocarbon phases is 2 at last NR step and predicted saturations
1509  // indicates that the stae will be kept, then skip phase stability analysis,
1510  // carry phase splitting directly
1511  bId = n * numPhase;
1512  for (USI j = 0; j < numPhase; j++) {
1513  if (dSNR[bId + j] + dSNRP[bId + j] < 1E-4) {
1514  // phases change (predicted)
1515  ftype = 0;
1516  break;
1517  }
1518  }
1519  }
1520  else {
1521  ftype = 0;
1522  }
1523  }
1524 
1525  return ftype;
1526 }
1527 
1528 
1529 void Bulk::PassFlashValue(const OCP_USI& n)
1530 {
1531  OCP_FUNCNAME;
1532 
1533  OCP_USI bIdp = n * numPhase;
1534  USI pvtnum = PVTNUM[n];
1535  USI nptmp = 0;
1536  for (USI j = 0; j < numPhase; j++) {
1537  phaseExist[bIdp + j] = flashCal[pvtnum]->phaseExist[j];
1538  // Important! Saturation must be passed no matter if the phase exists. This is
1539  // because it will be used to calculate relative permeability and capillary
1540  // pressure at each time step. Make sure that all saturations are updated at
1541  // each step!
1542  S[bIdp + j] = flashCal[pvtnum]->S[j];
1543  if (phaseExist[bIdp + j]) { // j -> bId + j fix bugs.
1544  nptmp++;
1545  rho[bIdp + j] = flashCal[pvtnum]->rho[j];
1546  xi[bIdp + j] = flashCal[pvtnum]->xi[j];
1547  for (USI i = 0; i < numCom; i++) {
1548  xij[bIdp * numCom + j * numCom + i] =
1549  flashCal[pvtnum]->xij[j * numCom + i];
1550  }
1551  mu[bIdp + j] = flashCal[pvtnum]->mu[j];
1552  vj[bIdp + j] = flashCal[pvtnum]->v[j];
1553  }
1554  }
1555  Nt[n] = flashCal[pvtnum]->Nt;
1556  vf[n] = flashCal[pvtnum]->vf;
1557  vfp[n] = flashCal[pvtnum]->vfp;
1558  OCP_USI bIdc = n * numCom;
1559  for (USI i = 0; i < numCom; i++) {
1560  vfi[bIdc + i] = flashCal[pvtnum]->vfi[i];
1561  }
1562 
1563  phaseNum[n] = nptmp - 1; // water is excluded
1564  if (comps) {
1565  if (nptmp == 3) {
1566  // num of hydrocarbon phase equals 2
1567  // Calculate Ks
1568  OCP_USI bIdc1 = n * numCom_1;
1569  for (USI i = 0; i < numCom_1; i++) {
1570  Ks[bIdc1 + i] =
1571  flashCal[pvtnum]->xij[i] / flashCal[pvtnum]->xij[numCom + i];
1572  }
1573  }
1574 
1575  if (flashCal[pvtnum]->GetFtype() == 0) {
1576  flagSkip[n] = flashCal[pvtnum]->GetFlagSkip();
1577  if (flagSkip[n]) {
1578  minEigenSkip[n] = flashCal[pvtnum]->GetMinEigenSkip();
1579  for (USI j = 0; j < numPhase - 1; j++) {
1580  if (phaseExist[bIdp + j]) {
1581  for (USI i = 0; i < numCom - 1; i++) {
1582  ziSkip[bIdc + i] = flashCal[pvtnum]->xij[j * numCom + i];
1583  }
1584  break;
1585  }
1586  }
1587  PSkip[n] = P[n];
1588  }
1589  }
1590 
1591  if (miscible) {
1592  surTen[n] = flashCal[pvtnum]->GetSurTen();
1593  }
1594  }
1595 }
1596 
1597 void Bulk::PassFlashValueAIMc(const OCP_USI& n)
1598 {
1599  // only var about volume needs, some flash var also
1600  OCP_FUNCNAME;
1601 
1602  OCP_USI bIdp = n * numPhase;
1603  USI pvtnum = PVTNUM[n];
1604 
1605  Nt[n] = flashCal[pvtnum]->Nt;
1606  vf[n] = flashCal[pvtnum]->vf;
1607  vfp[n] = flashCal[pvtnum]->vfp;
1608  OCP_USI bIdc = n * numCom;
1609  for (USI i = 0; i < numCom; i++) {
1610  vfi[bIdc + i] = flashCal[pvtnum]->vfi[i];
1611  }
1612 
1613  USI nptmp = 0;
1614  for (USI j = 0; j < numPhase; j++) {
1615  if (flashCal[pvtnum]->phaseExist[j]) {
1616  nptmp++;
1617  }
1618  }
1619 
1620  phaseNum[n] = nptmp - 1; // water is excluded
1621 
1622  if (comps) {
1623  if (nptmp == 3) {
1624  // num of hydrocarbon phase equals 2
1625  // Calculate Ks
1626  OCP_USI bIdc1 = n * numCom_1;
1627  for (USI i = 0; i < numCom_1; i++) {
1628  Ks[bIdc1 + i] =
1629  flashCal[pvtnum]->xij[i] / flashCal[pvtnum]->xij[numCom + i];
1630  }
1631  }
1632 
1633  if (flashCal[pvtnum]->GetFtype() == 0) {
1634  flagSkip[n] = flashCal[pvtnum]->GetFlagSkip();
1635  if (flagSkip[n]) {
1636  minEigenSkip[n] = flashCal[pvtnum]->GetMinEigenSkip();
1637  for (USI j = 0; j < numPhase - 1; j++) {
1638  if (phaseExist[bIdp + j]) {
1639  for (USI i = 0; i < numCom - 1; i++) {
1640  ziSkip[bIdc + i] = flashCal[pvtnum]->xij[j * numCom + i];
1641  }
1642  break;
1643  }
1644  }
1645  PSkip[n] = P[n];
1646  }
1647  }
1648 
1649  if (miscible) {
1650  surTen[n] = flashCal[pvtnum]->GetSurTen();
1651  }
1652  }
1653 }
1654 
1655 void Bulk::PassFlashValueDeriv(const OCP_USI& n)
1656 {
1657  OCP_FUNCNAME;
1658 
1659  OCP_USI bIdp = n * numPhase;
1660  USI pvtnum = PVTNUM[n];
1661  USI nptmp = 0;
1662  USI len = 0;
1663 
1664  for (USI j = 0; j < numPhase; j++) {
1665  // Important! Saturation must be passed no matter if the phase exists. This is
1666  // because it will be used to calculate relative permeability and capillary
1667  // pressure at each time step. Make sure that all saturations are updated at
1668  // each step!
1669  S[bIdp + j] = flashCal[pvtnum]->S[j];
1670  dSNR[bIdp + j] = S[bIdp + j] - dSNR[bIdp + j];
1671  if (phaseExist[bIdp + j]) {
1672  NRdSSP +=
1673  (dSNR[bIdp + j] - dSNRP[bIdp + j]) * (dSNR[bIdp + j] - dSNRP[bIdp + j]);
1674  if (fabs(maxNRdSSP) < fabs(dSNR[bIdp + j] - dSNRP[bIdp + j])) {
1675  maxNRdSSP = dSNR[bIdp + j] - dSNRP[bIdp + j];
1676  index_maxNRdSSP = n;
1677  }
1678  // cout << n << " " << scientific << setprecision(6) << dSNR[bIdp + j] <<
1679  // " " << dSNRP[bIdp + j] << endl;
1680  }
1681 
1682  phaseExist[bIdp + j] = flashCal[pvtnum]->phaseExist[j];
1683  pEnumCom[bIdp + j] = flashCal[pvtnum]->pEnumCom[j];
1684  len += pEnumCom[bIdp + j];
1685  if (phaseExist[bIdp + j]) { // j -> bId + j fix bugs.
1686  nptmp++;
1687  nj[bIdp + j] = flashCal[pvtnum]->nj[j];
1688  rho[bIdp + j] = flashCal[pvtnum]->rho[j];
1689  xi[bIdp + j] = flashCal[pvtnum]->xi[j];
1690  mu[bIdp + j] = flashCal[pvtnum]->mu[j];
1691  vj[bIdp + j] = flashCal[pvtnum]->v[j];
1692 
1693  // Derivatives
1694  muP[bIdp + j] = flashCal[pvtnum]->muP[j];
1695  xiP[bIdp + j] = flashCal[pvtnum]->xiP[j];
1696  rhoP[bIdp + j] = flashCal[pvtnum]->rhoP[j];
1697  for (USI i = 0; i < numCom; i++) {
1698  xij[bIdp * numCom + j * numCom + i] =
1699  flashCal[pvtnum]->xij[j * numCom + i];
1700  mux[bIdp * numCom + j * numCom + i] =
1701  flashCal[pvtnum]->mux[j * numCom + i];
1702  xix[bIdp * numCom + j * numCom + i] =
1703  flashCal[pvtnum]->xix[j * numCom + i];
1704  rhox[bIdp * numCom + j * numCom + i] =
1705  flashCal[pvtnum]->rhox[j * numCom + i];
1706  }
1707  }
1708  }
1709  Nt[n] = flashCal[pvtnum]->Nt;
1710  vf[n] = flashCal[pvtnum]->vf;
1711  vfp[n] = flashCal[pvtnum]->vfp;
1712 
1713  OCP_USI bIdc = n * numCom;
1714  for (USI i = 0; i < numCom; i++) {
1715  vfi[bIdc + i] = flashCal[pvtnum]->vfi[i];
1716  }
1717 
1718  len += nptmp;
1719 
1720 #ifdef OCP_NEW_FIM
1721  len *= (numCom + 1);
1722  dSdPindex[n + 1] = dSdPindex[n] + len;
1723  Dcopy(len, &dSec_dPri[0] + dSdPindex[n], &flashCal[pvtnum]->dXsdXp[0]);
1724 #else
1725  Dcopy(lendSdP, &dSec_dPri[0] + n * lendSdP, &flashCal[pvtnum]->dXsdXp[0]);
1726 #endif // OCP_NEW_FIM
1727 
1728  // test
1729  phaseNum[n] = nptmp - 1; // So water must exist!!!
1730 
1731  if (comps) {
1732  ePEC[n] = flashCal[pvtnum]->GetErrorPEC();
1733  if (nptmp == 3) {
1734  // num of hydrocarbon phase equals 2
1735  // Calculate Ks
1736  // IMPORTANT!!!
1737  // Ks will change as long as nums of hydroncarbon phase equals 2, and it
1738  // will has an effect on phase spliting calculation as a intial value. So
1739  // you should not expect to obtain the exact same result with identical P,
1740  // T, Ni if the final mixture contains 2 hydroncarbon phase.
1741  OCP_USI bIdc1 = n * numCom_1;
1742  for (USI i = 0; i < numCom_1; i++) {
1743  Ks[bIdc1 + i] =
1744  flashCal[pvtnum]->xij[i] / flashCal[pvtnum]->xij[numCom + i];
1745  }
1746  }
1747 
1748  if (flashCal[pvtnum]->GetFtype() == 0) {
1749  flagSkip[n] = flashCal[pvtnum]->GetFlagSkip();
1750  if (flagSkip[n]) {
1751  minEigenSkip[n] = flashCal[pvtnum]->GetMinEigenSkip();
1752  for (USI j = 0; j < numPhase - 1; j++) {
1753  if (phaseExist[bIdp + j]) {
1754  for (USI i = 0; i < numCom - 1; i++) {
1755  ziSkip[bIdc + i] = flashCal[pvtnum]->xij[j * numCom + i];
1756  }
1757  break;
1758  }
1759  }
1760  PSkip[n] = P[n];
1761  }
1762  }
1763  if (miscible) {
1764  surTen[n] = flashCal[pvtnum]->GetSurTen();
1765  }
1766  }
1767 }
1768 
1769 void Bulk::PassFlashValueDeriv_n(const OCP_USI& n)
1770 {
1771  OCP_FUNCNAME;
1772 
1773  OCP_USI bIdp = n * numPhase;
1774  USI pvtnum = PVTNUM[n];
1775  USI nptmp = 0;
1776  USI len = 0;
1777 
1778  for (USI j = 0; j < numPhase; j++) {
1779  // Important! Saturation must be passed no matter if the phase exists. This is
1780  // because it will be used to calculate relative permeability and capillary
1781  // pressure at each time step. Make sure that all saturations are updated at
1782  // each step!
1783  S[bIdp + j] = flashCal[pvtnum]->S[j];
1784  dSNR[bIdp + j] = S[bIdp + j] - dSNR[bIdp + j];
1785  if (phaseExist[bIdp + j]) {
1786  NRdSSP +=
1787  (dSNR[bIdp + j] - dSNRP[bIdp + j]) * (dSNR[bIdp + j] - dSNRP[bIdp + j]);
1788  if (fabs(maxNRdSSP) < fabs(dSNR[bIdp + j] - dSNRP[bIdp + j])) {
1789  maxNRdSSP = dSNR[bIdp + j] - dSNRP[bIdp + j];
1790  index_maxNRdSSP = n;
1791  }
1792  }
1793 
1794  phaseExist[bIdp + j] = flashCal[pvtnum]->phaseExist[j];
1795  pEnumCom[bIdp + j] = flashCal[pvtnum]->pEnumCom[j];
1796  len += pEnumCom[bIdp + j];
1797  if (phaseExist[bIdp + j]) { // j -> bId + j fix bugs.
1798  nptmp++;
1799  nj[bIdp + j] = flashCal[pvtnum]->nj[j];
1800  rho[bIdp + j] = flashCal[pvtnum]->rho[j];
1801  xi[bIdp + j] = flashCal[pvtnum]->xi[j];
1802  mu[bIdp + j] = flashCal[pvtnum]->mu[j];
1803  vj[bIdp + j] = flashCal[pvtnum]->v[j];
1804 
1805  // Derivatives
1806  muP[bIdp + j] = flashCal[pvtnum]->muP[j];
1807  xiP[bIdp + j] = flashCal[pvtnum]->xiP[j];
1808  rhoP[bIdp + j] = flashCal[pvtnum]->rhoP[j];
1809  for (USI i = 0; i < numCom; i++) {
1810  xij[bIdp * numCom + j * numCom + i] =
1811  flashCal[pvtnum]->xij[j * numCom + i];
1812  mux[bIdp * numCom + j * numCom + i] =
1813  flashCal[pvtnum]->mux[j * numCom + i];
1814  xix[bIdp * numCom + j * numCom + i] =
1815  flashCal[pvtnum]->xix[j * numCom + i];
1816  rhox[bIdp * numCom + j * numCom + i] =
1817  flashCal[pvtnum]->rhox[j * numCom + i];
1818  }
1819  }
1820  }
1821  Nt[n] = flashCal[pvtnum]->Nt;
1822  vf[n] = flashCal[pvtnum]->vf;
1823  vfp[n] = flashCal[pvtnum]->vfp;
1824 
1825  OCP_USI bIdc = n * numCom;
1826  for (USI i = 0; i < numCom; i++) {
1827  vfi[bIdc + i] = flashCal[pvtnum]->vfi[i];
1828  }
1829 
1830  len += nptmp;
1831  resIndex[n + 1] = resIndex[n] + len;
1832  Dcopy(len, &res_n[0] + resIndex[n], &flashCal[pvtnum]->res[0]);
1833  len *= (numCom + 1);
1834  dSdPindex[n + 1] = dSdPindex[n] + len;
1835  Dcopy(len, &dSec_dPri[0] + dSdPindex[n], &flashCal[pvtnum]->dXsdXp[0]);
1836 
1837  resPc[n] = flashCal[pvtnum]->resPc;
1838 
1839  // test
1840  phaseNum[n] = nptmp - 1; // So water must exist!!!
1841 
1842  if (comps) {
1843  ePEC[n] = flashCal[pvtnum]->GetErrorPEC();
1844  if (nptmp == 3) {
1845  // num of hydrocarbon phase equals 2
1846  // Calculate Ks
1847  // IMPORTANT!!!
1848  // Ks will change as long as nums of hydroncarbon phase equals 2, and it
1849  // will has an effect on phase spliting calculation as a intial value. So
1850  // you should not expect to obtain the exact same result with identical P,
1851  // T, Ni if the final mixture contains 2 hydroncarbon phase.
1852  OCP_USI bIdc1 = n * numCom_1;
1853  for (USI i = 0; i < numCom_1; i++) {
1854  Ks[bIdc1 + i] =
1855  flashCal[pvtnum]->xij[i] / flashCal[pvtnum]->xij[numCom + i];
1856  }
1857  }
1858 
1859  if (flashCal[pvtnum]->GetFtype() == 0) {
1860  flagSkip[n] = flashCal[pvtnum]->GetFlagSkip();
1861  if (flagSkip[n]) {
1862  minEigenSkip[n] = flashCal[pvtnum]->GetMinEigenSkip();
1863  for (USI j = 0; j < numPhase - 1; j++) {
1864  if (phaseExist[bIdp + j]) {
1865  for (USI i = 0; i < numCom - 1; i++) {
1866  ziSkip[bIdc + i] = flashCal[pvtnum]->xij[j * numCom + i];
1867  }
1868  break;
1869  }
1870  }
1871  PSkip[n] = P[n];
1872  }
1873  }
1874  if (miscible) {
1875  surTen[n] = flashCal[pvtnum]->GetSurTen();
1876  }
1877  }
1878 }
1879 
1880 void Bulk::ResetFlash()
1881 {
1882  OCP_FUNCNAME;
1883 
1884  phaseExist = lphaseExist;
1885  S = lS;
1886  rho = lrho;
1887  xi = lxi;
1888  xij = lxij;
1889  mu = lmu;
1890  vj = lvj;
1891  vf = lvf;
1892  vfp = lvfp;
1893  vfi = lvfi;
1894 }
1895 
1897 void Bulk::CalKrPc()
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 }
1928 
1929 void Bulk::CalKrPcDeriv()
1930 {
1931  OCP_FUNCNAME;
1932 
1933  if (!miscible) {
1934  OCP_DBL tmp = 0;
1935  for (OCP_USI n = 0; n < numBulk; n++) {
1936  OCP_USI bId = n * numPhase;
1937  flow[SATNUM[n]]->CalKrPcDeriv(&S[bId], &kr[bId], &Pc[bId],
1938  &dKr_dS[bId * numPhase],
1939  &dPcj_dS[bId * numPhase], 0, tmp, tmp);
1940  for (USI j = 0; j < numPhase; j++) Pj[bId + j] = P[n] + Pc[bId + j];
1941  }
1942  } else {
1943  for (OCP_USI n = 0; n < numBulk; n++) {
1944  OCP_USI bId = n * numPhase;
1945  flow[SATNUM[n]]->CalKrPcDeriv(
1946  &S[bId], &kr[bId], &Pc[bId], &dKr_dS[bId * numPhase],
1947  &dPcj_dS[bId * numPhase], surTen[n], Fk[n], Fp[n]);
1948  for (USI j = 0; j < numPhase; j++)
1949  Pj[n * numPhase + j] = P[n] + Pc[n * numPhase + j];
1950  }
1951  }
1952 
1953  if (ScalePcow) {
1954  // correct
1955  USI Wid = phase2Index[WATER];
1956  for (USI n = 0; n < numBulk; n++) {
1957  Pc[n * numPhase + Wid] *= ScaleValuePcow[n];
1958  Pj[n * numPhase + Wid] = P[n] + Pc[n * numPhase + Wid];
1959  }
1960  }
1961 }
1962 
1963 void Bulk::CalVpore()
1964 {
1965  OCP_FUNCNAME;
1966 
1967  for (OCP_USI n = 0; n < numBulk; n++) {
1968  OCP_DBL dP = rockC1 * (P[n] - rockPref);
1969  rockVp[n] = rockVpInit[n] * (1 + dP);
1970  // rockVp[n] = rockVpInit[n] * (1 + dP + dP * dP / 2);
1971  }
1972 }
1973 
1974 OCP_DBL Bulk::CalFPR() const
1975 {
1976  OCP_FUNCNAME;
1977 
1978  OCP_DBL ptmp = 0;
1979  OCP_DBL vtmp = 0;
1980  OCP_DBL tmp = 0;
1981 
1982  if (numPhase == 3) {
1983  for (OCP_USI n = 0; n < numBulk; n++) {
1984  tmp = rockVp[n] * (1 - S[n * numPhase + 2]);
1985  ptmp += P[n] * tmp;
1986  vtmp += tmp;
1987  }
1988  } else if (numPhase < 3) {
1989  for (OCP_USI n = 0; n < numBulk; n++) {
1990  tmp = rockVp[n] * (S[n * numPhase]);
1991  ptmp += P[n] * tmp;
1992  vtmp += tmp;
1993  }
1994  } else {
1995  OCP_ABORT("Number of phases is out of range!");
1996  }
1997  return ptmp / vtmp;
1998 }
1999 
2000 void Bulk::CalMaxChange()
2001 {
2002  OCP_FUNCNAME;
2003 
2004  dPmax = 0;
2005  dNmax = 0;
2006  dSmax = 0;
2007  dVmax = 0;
2008  OCP_DBL tmp = 0;
2009  OCP_USI id;
2010  OCP_USI ndPmax, ndNmax, ndSmax, ndVmax;
2011 
2012  for (OCP_USI n = 0; n < numBulk; n++) {
2013 
2014  // dP
2015  tmp = fabs(P[n] - lP[n]);
2016  if (dPmax < tmp) {
2017  dPmax = tmp;
2018  ndPmax = n;
2019  }
2020 
2021  // dS
2022  for (USI j = 0; j < numPhase; j++) {
2023  id = n * numPhase + j;
2024  tmp = fabs(S[id] - lS[id]);
2025  if (dSmax < tmp) {
2026  dSmax = tmp;
2027  ndSmax = n;
2028  }
2029  }
2030 
2031  // dN
2032  for (USI i = 0; i < numCom; i++) {
2033  id = n * numCom + i;
2034 
2035  tmp = fabs(max(Ni[id], lNi[id]));
2036  if (tmp > TINY) {
2037  tmp = fabs(Ni[id] - lNi[id]) / tmp;
2038  if (dNmax < tmp) {
2039  dNmax = tmp;
2040  ndNmax = n;
2041  }
2042  }
2043  }
2044 
2045  tmp = fabs(vf[n] - rockVp[n]) / rockVp[n];
2046  if (dVmax < tmp) {
2047  dVmax = tmp;
2048  ndVmax = n;
2049  }
2050  }
2051 
2052  // cout << scientific << setprecision(6);
2053  // cout << setw(15) << dPmax << setw(15) << dNmax << setw(15) << dSmax << setw(15)
2054  // << dVmax << endl; cout << setw(15) << ndPmax << setw(15) << ndNmax << setw(15) <<
2055  // ndSmax << setw(15) << ndVmax << endl;
2056 
2057  // cout << ndSmax << endl;
2058  // cout << "old " << lP[ndSmax] << endl;
2059  // tmp = Dnorm1(numCom_1, &lNi[ndSmax * numCom]);
2060  // for (USI i = 0; i < numCom_1; i++) {
2061  // cout << lNi[ndSmax * numCom + i] << " ";
2062  // }
2063  // cout << endl;
2064  // for (USI i = 0; i < numCom_1; i++) {
2065  // cout << lNi[ndSmax * numCom + i] / tmp << " ";
2066  // }
2067  // cout << endl;
2068  // for (USI j = 0; j < numPhase - 1; j++) {
2069  // cout << lS[ndSmax * numPhase + j] << " ";
2070  // for (USI i = 0; i < numCom_1; i++) {
2071  // cout << lxij[(ndSmax * numPhase + j) * numCom + i] << " ";
2072  // }
2073  // cout << endl;
2074  // }
2075  // cout << "new " << P[ndSmax] << endl;
2076  // tmp = Dnorm1(numCom_1, &Ni[ndSmax * numCom]);
2077  // for (USI i = 0; i < numCom_1; i++) {
2078  // cout << Ni[ndSmax * numCom + i] << " ";
2079  // }
2080  // cout << endl;
2081  // for (USI i = 0; i < numCom_1; i++) {
2082  // cout << Ni[ndSmax * numCom + i] / tmp << " ";
2083  // }
2084  // cout << endl;
2085  // for (USI j = 0; j < numPhase - 1; j++) {
2086  // cout << S[ndSmax * numPhase + j] << " ";
2087  // for (USI i = 0; i < numCom_1; i++) {
2088  // cout << xij[(ndSmax * numPhase + j) * numCom + i] << " ";
2089  // }
2090  // cout << endl;
2091  // }
2092 }
2093 
2095 bool Bulk::CheckP() const
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 }
2112 
2114 bool Bulk::CheckNi()
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 }
2158 
2160 bool Bulk::CheckVe(const OCP_DBL& Vlim) const
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 }
2181 
2182 void Bulk::CheckDiff()
2183 {
2184  OCP_FUNCNAME;
2185 
2186  OCP_DBL tmp;
2187  OCP_USI id;
2188  for (OCP_USI n = 0; n < numBulk; n++) {
2189  for (USI j = 0; j < numPhase; j++) {
2190  id = n * numPhase + j;
2191  tmp = fabs(phaseExist[id] - lphaseExist[id]);
2192  if (tmp >= 1E-10) {
2193  cout << "Difference in phaseExist\t" << tmp << "\n";
2194  }
2195  if (lphaseExist[id] || phaseExist[id]) {
2196  tmp = fabs(S[id] - lS[id]);
2197  if (tmp >= 1E-10) {
2198  cout << scientific << setprecision(16);
2199  cout << "Bulk[" << n << "]" << endl;
2200  cout << "Saturation" << endl;
2201  for (USI j = 0; j < numPhase; j++) {
2202  cout << S[n * numPhase + j] << " " << lS[n * numPhase + j]
2203  << endl;
2204  }
2205  cout << "Pressure" << endl;
2206  cout << P[n] << " " << lP[n] << endl;
2207  cout << "Ni" << endl;
2208  for (USI i = 0; i < numCom; i++) {
2209  cout << Ni[n * numCom + i] << " " << lNi[n * numCom + i]
2210  << endl;
2211  }
2212  cout << "PhaseNum" << endl;
2213  cout << phaseNum[n] << " " << lphaseNum[n] << endl;
2214  cout << "minEigenSkip" << endl;
2215  cout << minEigenSkip[n] << " " << lminEigenSkip[n] << endl;
2216  cout << "flagSkip" << endl;
2217  cout << flagSkip[n] << " " << lflagSkip[n] << endl;
2218  cout << "PSkip" << endl;
2219  cout << PSkip[n] << " " << lPSkip[n] << endl;
2220  cout << "ziSkip" << endl;
2221  for (USI i = 0; i < numCom; i++) {
2222  cout << ziSkip[n * numCom + i] << " "
2223  << lziSkip[n * numCom + i] << endl;
2224  }
2225  cout << "Ks" << endl;
2226  for (USI i = 0; i < numCom_1; i++) {
2227  cout << Ks[n * numCom_1 + i] << " " << lKs[n * numCom_1 + i]
2228  << endl;
2229  }
2230 
2231  cout << "Difference in S\t" << tmp << " " << phaseExist[id]
2232  << "\n";
2233  }
2234  tmp = fabs(xi[id] - lxi[id]);
2235  if (tmp >= 1E-10) {
2236  cout << "Difference in Xi\t" << tmp << " " << phaseExist[id]
2237  << "\n";
2238  }
2239  tmp = fabs(rho[id] - lrho[id]);
2240  if (tmp >= 1E-10) {
2241  cout << "Difference in rho\t" << tmp << " " << phaseExist[id]
2242  << "\n";
2243  }
2244  }
2245  }
2246  }
2247 }
2248 
2249 void Bulk::CheckSat() const
2250 {
2251  OCP_FUNCNAME;
2252 
2253  OCP_DBL tmp;
2254  for (OCP_USI n = 0; n < numBulk; n++) {
2255  tmp = 0.0;
2256  for (USI j = 0; j < numPhase; j++) {
2257  if (phaseExist[n * numPhase + j]) {
2258  if (S[n * numPhase + j] < 0) {
2259  OCP_ABORT("Negative volume!");
2260  }
2261  tmp += S[n * numPhase + j];
2262  }
2263  }
2264  if (fabs(tmp - 1) > TINY) {
2265  OCP_ABORT("Saturation is greater than 1!");
2266  }
2267  }
2268 }
2269 
2270 USI Bulk::GetMixMode() const
2271 {
2272  OCP_FUNCNAME;
2273 
2274  if (blackOil)
2275  return BLKOIL;
2276  else if (comps)
2277  return EOS_PVTW;
2278  else
2279  OCP_ABORT("Mixture model is not supported!");
2280 }
2281 
2282 void Bulk::CalSomeInfo(const Grid& myGrid) const
2283 {
2284  // test
2285  OCP_DBL depthMax = 0;
2286  OCP_USI ndepa = 0;
2287  OCP_DBL depthMin = 1E8;
2288  OCP_USI ndepi = 0;
2289  OCP_DBL dxMax = 0;
2290  OCP_USI nxa = 0;
2291  OCP_DBL dxMin = 1E8;
2292  OCP_USI nxi = 0;
2293  OCP_DBL dyMax = 0;
2294  OCP_USI nya = 0;
2295  OCP_DBL dyMin = 1E8;
2296  OCP_USI nyi = 0;
2297  OCP_DBL dzMax = 0;
2298  OCP_USI nza = 0;
2299  OCP_DBL dzMin = 1E8;
2300  OCP_USI nzi = 0;
2301  OCP_DBL RVMax = 0;
2302  OCP_USI nRVa = 0;
2303  OCP_DBL RVMin = 1E8;
2304  OCP_USI nRVi = 0;
2305  OCP_DBL RVPMax = 0;
2306  OCP_USI nRVPa = 0;
2307  OCP_DBL RVPMin = 1E8;
2308  OCP_USI nRVPi = 0;
2309  OCP_DBL PerxMax = 0;
2310  OCP_USI nPerxa = 0;
2311  OCP_DBL PerxMin = 1E8;
2312  OCP_USI nPerxi = 0;
2313  USI I, J, K;
2314  const USI sp = myGrid.GetNumDigitIJK();
2315  for (OCP_USI n = 0; n < numBulk; n++) {
2316  // if (!activeMap_G2B[nn].IsAct())
2317  // continue;
2318  // OCP_USI n = activeMap_G2B[nn].GetId();
2319  if (depthMax < depth[n]) {
2320  depthMax = depth[n];
2321  ndepa = n;
2322  }
2323  if (depthMin > depth[n]) {
2324  depthMin = depth[n];
2325  ndepi = n;
2326  }
2327  if (dxMax < dx[n]) {
2328  dxMax = dx[n];
2329  nxa = n;
2330  }
2331  if (dxMin > dx[n]) {
2332  dxMin = dx[n];
2333  nxi = n;
2334  }
2335  if (dyMax < dy[n]) {
2336  dyMax = dy[n];
2337  nya = n;
2338  }
2339  if (dyMin > dy[n]) {
2340  dyMin = dy[n];
2341  nyi = n;
2342  }
2343  if (dzMax < dz[n]) {
2344  dzMax = dz[n];
2345  nza = n;
2346  }
2347  if (dzMin > dz[n]) {
2348  dzMin = dz[n];
2349  nzi = n;
2350  }
2351  OCP_DBL tmp = myGrid.v[myGrid.activeMap_B2G[n]];
2352  if (RVMax < tmp) {
2353  RVMax = tmp;
2354  nRVa = n;
2355  }
2356  if (RVMin > tmp) {
2357  RVMin = tmp;
2358  nRVi = n;
2359  }
2360  tmp = rockVp[n];
2361  if (RVPMax < tmp) {
2362  RVPMax = tmp;
2363  nRVPa = n;
2364  }
2365  if (RVPMin > tmp) {
2366  RVPMin = tmp;
2367  nRVPi = n;
2368  }
2369  tmp = rockKx[n];
2370  if (PerxMax < tmp) {
2371  PerxMax = tmp;
2372  nPerxa = n;
2373  }
2374  if (PerxMin > tmp && fabs(tmp) > 1E-8) {
2375  PerxMin = tmp;
2376  nPerxi = n;
2377  }
2378  }
2379 
2380  cout << "---------------------" << endl
2381  << "BULK" << endl
2382  << "---------------------" << endl;
2383  myGrid.GetIJKGrid(I, J, K, ndepa);
2384  cout << " Depthmax = " << depthMax << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2385  myGrid.GetIJKGrid(I, J, K, ndepi);
2386  cout << " Depthmin = " << depthMin << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2387  myGrid.GetIJKGrid(I, J, K, nxa);
2388  cout << " DXmax = " << dxMax << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2389  myGrid.GetIJKGrid(I, J, K, nxi);
2390  cout << " DXmin = " << dxMin << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2391  myGrid.GetIJKGrid(I, J, K, nya);
2392  cout << " DYmax = " << dyMax << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2393  myGrid.GetIJKGrid(I, J, K, nyi);
2394  cout << " DYmin = " << dyMin << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2395  myGrid.GetIJKGrid(I, J, K, nza);
2396  cout << " DZmax = " << dzMax << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2397  myGrid.GetIJKGrid(I, J, K, nzi);
2398  cout << " DZmin = " << dzMin << " feet " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2399  myGrid.GetIJKGrid(I, J, K, nRVa);
2400  cout << " RVmax = " << RVMax / CONV1 << " rb " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2401  myGrid.GetIJKGrid(I, J, K, nRVi);
2402  cout << " RVmin = " << RVMin / CONV1 << " rb " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2403  myGrid.GetIJKGrid(I, J, K, nRVPa);
2404  cout << " RVmax = " << RVPMax / CONV1 << " rb " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2405  myGrid.GetIJKGrid(I, J, K, nRVPi);
2406  cout << " RVmin = " << RVPMin / CONV1 << " rb " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2407  myGrid.GetIJKGrid(I, J, K, nPerxa);
2408  cout << " Perxmax = " << PerxMax << " " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2409  myGrid.GetIJKGrid(I, J, K, nPerxi);
2410  cout << " Perxmin = " << scientific << PerxMin << " " << GetIJKformat(to_string(I), to_string(J), to_string(K), sp) << endl;
2411 }
2412 
2414 // IMPEC
2416 
2417 void Bulk::AllocateAuxIMPEC()
2418 {
2419  OCP_FUNCNAME;
2420 
2421  // IMPEC var
2422  vfi.resize(numBulk * numCom);
2423  vfp.resize(numBulk);
2424  cfl.resize(numBulk * numPhase);
2425 
2426  // Last Step
2427  lP.resize(numBulk);
2428  lPj.resize(numBulk * numPhase);
2429  lPc.resize(numBulk * numPhase);
2430  lphaseExist.resize(numBulk * numPhase);
2431  lS.resize(numBulk * numPhase);
2432  lrho.resize(numBulk * numPhase);
2433  lxi.resize(numBulk * numPhase);
2434  lxij.resize(numBulk * numPhase * numCom);
2435  lNi.resize(numBulk * numCom);
2436  lmu.resize(numBulk * numPhase);
2437  lkr.resize(numBulk * numPhase);
2438  lvj.resize(numBulk * numPhase);
2439  lvf.resize(numBulk);
2440  lNt.resize(numBulk);
2441  lvfi.resize(numBulk * numCom);
2442  lvfp.resize(numBulk);
2443  lrockVp.resize(numBulk);
2444 }
2445 
2446 void Bulk::GetSolIMPEC(const vector<OCP_DBL>& u)
2447 {
2448  OCP_FUNCNAME;
2449 
2450  for (OCP_USI n = 0; n < numBulk; n++) {
2451  P[n] = u[n];
2452  for (USI j = 0; j < numPhase; j++) {
2453  OCP_USI id = n * numPhase + j;
2454  if (phaseExist[id]) Pj[id] = P[n] + Pc[id];
2455  }
2456  }
2457 }
2458 
2459 OCP_DBL Bulk::CalCFL() const
2460 {
2461  OCP_FUNCNAME;
2462 
2463  OCP_DBL tmp = 0;
2464  OCP_USI id;
2465  for (OCP_USI n = 0; n < numBulk; n++) {
2466  for (USI j = 0; j < numPhase; j++) {
2467  id = n * numPhase + j;
2468  if (phaseExist[id]) {
2469  if (vj[id] <= 0) continue; // temp
2470  cfl[id] /= vj[id];
2471 
2472 #ifdef DEBUG
2473  if (!isfinite(cfl[id])) {
2474  OCP_ABORT("cfl is nan!");
2475  }
2476 #endif // DEBUG
2477 
2478  if (tmp < cfl[id]) tmp = cfl[id];
2479  }
2480  }
2481  }
2482  return tmp;
2483 }
2484 
2485 void Bulk::UpdateLastStepIMPEC()
2486 {
2487  OCP_FUNCNAME;
2488  lphaseNum = phaseNum;
2489  lminEigenSkip = minEigenSkip;
2490  lflagSkip = flagSkip;
2491  lziSkip = ziSkip;
2492  lPSkip = PSkip;
2493  lKs = Ks;
2494 
2495  lP = P;
2496  lPj = Pj;
2497  lPc = Pc;
2498  lphaseExist = phaseExist;
2499  lS = S;
2500  lrho = rho;
2501  lxi = xi;
2502  lxij = xij;
2503  lNi = Ni;
2504  lmu = mu;
2505  lkr = kr;
2506  lvj = vj;
2507  lvf = vf;
2508  lvfi = vfi;
2509  lvfp = vfp;
2510  lrockVp = rockVp;
2511  lNt = Nt;
2512 
2513  if (miscible) {
2514  lsurTen = surTen;
2515  }
2516 }
2517 
2519 // FIM
2521 
2522 void Bulk::AllocateAuxFIM()
2523 {
2524  OCP_FUNCNAME;
2525 
2526  // FIM var
2527  nj.resize(numBulk * numPhase);
2528  vfi.resize(numBulk * numCom);
2529  vfp.resize(numBulk);
2530  muP.resize(numBulk * numPhase);
2531  xiP.resize(numBulk * numPhase);
2532  rhoP.resize(numBulk * numPhase);
2533  mux.resize(numBulk * numCom * numPhase);
2534  xix.resize(numBulk * numCom * numPhase);
2535  rhox.resize(numBulk * numCom * numPhase);
2536  lendSdP = (numCom + 1) * (numCom + 1) * numPhase;
2537  dSec_dPri.resize(numBulk * lendSdP);
2538  res_n.resize((numPhase + numPhase * numCom) * numBulk);
2539  resPc.resize(numBulk);
2540  dSdPindex.resize(numBulk + 1, 0);
2541  resIndex.resize(numBulk + 1, 0);
2542  pEnumCom.resize(numBulk * numPhase);
2543  dKr_dS.resize(numBulk * numPhase * numPhase);
2544  dPcj_dS.resize(numBulk * numPhase * numPhase);
2545 
2546  // Last Step
2547  lP.resize(numBulk);
2548  lPj.resize(numBulk * numPhase);
2549  lPc.resize(numBulk * numPhase);
2550  lphaseExist.resize(numBulk * numPhase);
2551  lS.resize(numBulk * numPhase);
2552  lnj.resize(numBulk * numPhase);
2553  lrho.resize(numBulk * numPhase);
2554  lxi.resize(numBulk * numPhase);
2555  lxij.resize(numBulk * numPhase * numCom);
2556  lNi.resize(numBulk * numCom);
2557  lmu.resize(numBulk * numPhase);
2558  lkr.resize(numBulk * numPhase);
2559  lvj.resize(numBulk * numPhase);
2560  lvf.resize(numBulk);
2561  lNt.resize(numBulk);
2562  lvfi.resize(numBulk * numCom);
2563  lvfp.resize(numBulk);
2564  lrockVp.resize(numBulk);
2565  lmuP.resize(numBulk * numPhase);
2566  lxiP.resize(numBulk * numPhase);
2567  lrhoP.resize(numBulk * numPhase);
2568  lmux.resize(numBulk * numCom * numPhase);
2569  lxix.resize(numBulk * numCom * numPhase);
2570  lrhox.resize(numBulk * numCom * numPhase);
2571  ldSec_dPri.resize(numBulk * lendSdP);
2572  lres_n.resize((numPhase + numPhase * numCom) * numBulk);
2573  ldSdPindex.resize(numBulk + 1, 0);
2574  lresIndex.resize(numBulk + 1, 0);
2575  lpEnumCom.resize(numBulk * numPhase);
2576  ldKr_dS.resize(numBulk * numPhase * numPhase);
2577  ldPcj_dS.resize(numBulk * numPhase * numPhase);
2578 
2579  dSNR.resize(numBulk * numPhase);
2580  dSNRP.resize(numBulk * numPhase);
2581  dNNR.resize(numBulk * numCom);
2582  dPNR.resize(numBulk);
2583 
2584  NRstep.resize(numBulk);
2585 }
2586 
2587 void Bulk::GetSolFIM(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
2588  const OCP_DBL& dSmaxlim)
2589 {
2590  OCP_FUNCNAME;
2591 
2592  dSNR = S;
2593  NRphaseNum = phaseNum;
2594 
2595  NRdSmaxP = 0;
2596  NRdPmax = 0;
2597  NRdNmax = 0;
2598  OCP_DBL dP;
2599  USI row0 = numPhase;
2600  USI row = numPhase * (numCom + 1);
2601  USI col = numCom + 1;
2602  USI bsize = row * col;
2603  OCP_USI n_np_j;
2604  vector<OCP_DBL> dtmp(row, 0);
2605  OCP_DBL chopmin = 1;
2606  OCP_DBL choptmp = 0;
2607 
2608  for (OCP_USI n = 0; n < numBulk; n++) {
2609  //const vector<OCP_DBL>& scm = satcm[SATNUM[n]];
2610 
2611  chopmin = 1;
2612  // compute the chop
2613  fill(dtmp.begin(), dtmp.end(), 0.0);
2614 #ifdef OCP_NEW_FIM
2615  DaAxpby((dSdPindex[n + 1] - dSdPindex[n]) / col, col, 1, dSec_dPri.data() + dSdPindex[n],
2616  u.data() + n * col, 1, dtmp.data());
2617  /*DaAxpby(phaseNum[n] + 1, col, 1, dSec_dPri.data() + dSdPindex[n],
2618  u.data() + n * col, 1, dtmp.data());*/
2619  const bool newFIM = true;
2620 #else
2621  DaAxpby(row0, col, 1, dSec_dPri.data() + n * bsize, u.data() + n * col, 1,
2622  dtmp.data());
2623  const bool newFIM = false;
2624 #endif // OCP_NEW_FIM
2625 
2626  USI js = 0;
2627  for (USI j = 0; j < numPhase; j++) {
2628  if (!phaseExist[n * numPhase + j] && newFIM) {
2629  continue;
2630  }
2631  n_np_j = n * numPhase + j;
2632 
2633  choptmp = 1;
2634  if (fabs(dtmp[js]) > dSmaxlim) {
2635  choptmp = dSmaxlim / fabs(dtmp[js]);
2636  } else if (S[n_np_j] + dtmp[js] < 0.0) {
2637  choptmp = 0.9 * S[n_np_j] / fabs(dtmp[js]);
2638  }
2639 
2640  //if (fabs(S[n_np_j] - scm[j]) > TINY &&
2641  // (S[n_np_j] - scm[j]) / (choptmp * dtmp[js]) < 0)
2642  // choptmp *= min(1.0, -((S[n_np_j] - scm[j]) / (choptmp * dtmp[js])));
2643 
2644  chopmin = min(chopmin, choptmp);
2645  js++;
2646  }
2647 
2648  // dS
2649  js = 0;
2650  for (USI j = 0; j < numPhase; j++) {
2651  if (!phaseExist[n * numPhase + j] && newFIM) {
2652  dSNRP[n * numPhase + j] = 0;
2653  continue;
2654  }
2655  dSNRP[n * numPhase + j] = chopmin * dtmp[js];
2656  if (fabs(NRdSmaxP) < fabs(dSNRP[n * numPhase + j]))
2657  NRdSmaxP = dSNRP[n * numPhase + j];
2658  js++;
2659  }
2660 
2661  // dxij ---- Compositional model only
2662  if (phaseNum[n] == 2 && comps) {
2663  bool tmpflag = true;
2664  OCP_USI bId = 0;
2665  for (USI j = 0; j < 2; j++) {
2666  bId = n * numPhase * numCom + j * numCom;
2667  for (USI i = 0; i < numCom_1; i++) {
2668  xij[bId + i] += chopmin * dtmp[js];
2669  js++;
2670  if (xij[bId + i] < 0) tmpflag = false;
2671  }
2672  }
2673  if (tmpflag || true) {
2674  bId = n * numPhase * numCom;
2675  for (USI i = 0; i < numCom_1; i++) {
2676  Ks[n * numCom_1 + i] = xij[bId + i] / xij[bId + numCom + i];
2677  }
2678  }
2679  }
2680 
2681 
2682  // dP
2683  dP = u[n * col];
2684  // choptmp = dPmaxlim / fabs(dP);
2685  // chopmin = min(chopmin, choptmp);
2686  if (fabs(NRdPmax) < fabs(dP)) NRdPmax = dP;
2687  P[n] += dP; // seems better
2688  dPNR[n] = dP;
2689 
2690 
2691  // dNi
2692  NRstep[n] = chopmin;
2693  for (USI i = 0; i < numCom; i++) {
2694  dNNR[n * numCom + i] = u[n * col + 1 + i] * chopmin;
2695  if (fabs(NRdNmax) < fabs(dNNR[n * numCom + i]) / Nt[n])
2696  NRdNmax = dNNR[n * numCom + i] / Nt[n];
2697 
2698  Ni[n * numCom + i] += dNNR[n * numCom + i];
2699  }
2700  // cout << scientific << setprecision(6) << dP << " " << n << endl;
2701  }
2702 }
2703 
2704 void Bulk::GetSolFIM_n(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
2705  const OCP_DBL& dSmaxlim)
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 }
2857 
2858 void Bulk::GetSol01FIM(const vector<OCP_DBL>& u) { GetSolAIMc(u, 1, 1); }
2859 
2860 void Bulk::CalRelResFIM(ResFIM& resFIM) const
2861 {
2862  OCP_FUNCNAME;
2863 
2864  // OCP_USI tmpid01 = -1;
2865  // OCP_USI tmpid02 = -1;
2866  OCP_DBL tmp;
2867 
2868  const USI len = numCom + 1;
2869  for (OCP_USI n = 0; n < numBulk; n++) {
2870 
2871  eV[n] = 0;
2872  for (USI i = 0; i < len; i++) {
2873  tmp = fabs(resFIM.res[n * len + i] / rockVp[n]);
2874  if (resFIM.maxRelRes_v < tmp) {
2875  resFIM.maxRelRes_v = tmp;
2876  resFIM.maxId_v = n;
2877  }
2878  eV[n] += tmp * tmp;
2879  }
2880  eV[n] = sqrt(eV[n]);
2881 
2882  eN[n] = 0;
2883  for (USI i = 1; i < len; i++) {
2884  tmp = fabs(resFIM.res[n * len + i] / Nt[n]);
2885  if (resFIM.maxRelRes_mol < tmp) {
2886  resFIM.maxRelRes_mol = tmp;
2887  resFIM.maxId_mol = n;
2888  }
2889  eN[n] += tmp * tmp;
2890  }
2891  eN[n] = sqrt(eN[n]);
2892  }
2893  // cout << scientific;
2894  // if (tmpid01 < numBulk) {
2895  // cout << "maxRelRes_v: " << tmpid01 << " " << S[tmpid01 * numPhase] << " "
2896  // << S[tmpid01 * numPhase + 1] << " " << S[tmpid01 * numPhase + 2] <<
2897  // " " << resFIM.maxRelRes_v << endl;
2898  // }
2899  // if (tmpid02 < numBulk) {
2900  // cout << "maxRelRes_mol: " << tmpid02 << " " << S[tmpid02 * numPhase] << " "
2901  // << S[tmpid02 * numPhase + 1] << " " << S[tmpid02 * numPhase + 2] <<
2902  // " " << resFIM.maxRelRes_mol << endl;
2903  // }
2904 }
2905 
2906 void Bulk::ShowRes(const vector<OCP_DBL>& res) const
2907 {
2908  OCP_USI bId;
2909 
2910  const USI len = numCom + 1;
2911  for (OCP_USI n = 0; n < numBulk; n++) {
2912  bId = n * len;
2913  cout << endl;
2914  cout << "Bulk[" << setw(3) << n << "] ";
2915  cout << "Vp " << rockVp[n] << " ";
2916  cout << "Nt " << Nt[n] << " ";
2917  cout << "NRstep " << NRstep[n] << " ";
2918  cout << "P " << P[n] << " ";
2919  cout << dPNR[n] << " " << dPNR[n] / P[n] << " ";
2920  cout << (P[n] - lP[n]) / P[n] << " ";
2921  cout << "PhaseNum " << phaseNum[n] << " ";
2922  cout << endl;
2923  for (USI i = 0; i < len; i++) {
2924 
2925  cout << "[" << setw(2) << i << "]"
2926  << " ";
2927  cout << -res[bId + i] << " ";
2928  cout << fabs(res[bId + i] / rockVp[n]) << " ";
2929  if (i > 0) {
2930  cout << fabs(res[bId + i] / Nt[n]) << " ";
2931  cout << Ni[n * numCom + i - 1] << " " << lNi[n * numCom + i - 1]
2932  << " ";
2933  cout << dNNR[n * numCom + i - 1] << " "
2934  << dNNR[n * numCom + i - 1] /
2935  (Ni[n * numCom + i - 1] - dNNR[n * numCom + i - 1])
2936  << " ";
2937  cout << (Ni[n * numCom + i - 1] - lNi[n * numCom + i - 1]) /
2938  lNi[n * numCom + i - 1]
2939  << " ";
2940  // cout << vfi[n * numCom + i] << " ";
2941  } else {
2942  cout << rockVp[n] - vf[n] << " ";
2943  }
2944 
2945  cout << endl;
2946  }
2947  }
2948 }
2949 
2950 void Bulk::ResetFIM()
2951 {
2952  OCP_FUNCNAME;
2953 
2954  phaseNum = lphaseNum;
2955  minEigenSkip = lminEigenSkip;
2956  flagSkip = lflagSkip;
2957  ziSkip = lziSkip;
2958  PSkip = lPSkip;
2959  Ks = lKs;
2960 
2961  P = lP;
2962  Pj = lPj;
2963  Pc = lPc;
2964  phaseExist = lphaseExist;
2965  S = lS;
2966  nj = lnj;
2967  rho = lrho;
2968  xi = lxi;
2969  xij = lxij;
2970  Ni = lNi;
2971  mu = lmu;
2972  kr = lkr;
2973  vj = lvj;
2974  vf = lvf;
2975  Nt = lNt;
2976  vfi = lvfi;
2977  vfp = lvfp;
2978  rockVp = lrockVp;
2979  muP = lmuP;
2980  xiP = lxiP;
2981  rhoP = lrhoP;
2982  mux = lmux;
2983  xix = lxix;
2984  rhox = lrhox;
2985  dSec_dPri = ldSec_dPri;
2986  res_n = lres_n;
2987  resPc = lresPc;
2988  dSdPindex = ldSdPindex;
2989  resIndex = lresIndex;
2990  pEnumCom = lpEnumCom;
2991  dKr_dS = ldKr_dS;
2992  dPcj_dS = ldPcj_dS;
2993 
2994  if (miscible) {
2995  surTen = lsurTen;
2996  }
2997 }
2998 
2999 void Bulk::UpdateLastStepFIM()
3000 {
3001  OCP_FUNCNAME;
3002  lphaseNum = phaseNum;
3003  lminEigenSkip = minEigenSkip;
3004  lflagSkip = flagSkip;
3005  lziSkip = ziSkip;
3006  lPSkip = PSkip;
3007  lKs = Ks;
3008 
3009  lP = P;
3010  lPj = Pj;
3011  lPc = Pc;
3012  lphaseExist = phaseExist;
3013  lS = S;
3014  lnj = nj;
3015  lrho = rho;
3016  lxi = xi;
3017  lxij = xij;
3018  lNi = Ni;
3019  lmu = mu;
3020  lkr = kr;
3021  lvj = vj;
3022  lvf = vf;
3023  lNt = Nt;
3024  lvfi = vfi;
3025  lvfp = vfp;
3026  lrockVp = rockVp;
3027  lmuP = muP;
3028  lxiP = xiP;
3029  lrhoP = rhoP;
3030  lmux = mux;
3031  lxix = xix;
3032  lrhox = rhox;
3033  ldSec_dPri = dSec_dPri;
3034  lres_n = res_n;
3035  lresPc = resPc;
3036  ldSdPindex = dSdPindex;
3037  lresIndex = resIndex;
3038  lpEnumCom = pEnumCom;
3039  ldKr_dS = dKr_dS;
3040  ldPcj_dS = dPcj_dS;
3041 
3042  if (miscible) {
3043  lsurTen = surTen;
3044  }
3045 }
3046 
3047 OCP_DBL Bulk::CalNRdSmax(OCP_USI& index)
3048 {
3049  NRdSmax = 0;
3050  OCP_USI len = numBulk * numPhase;
3051  for (USI n = 0; n < len; n++) {
3052  if (fabs(NRdSmax) < fabs(dSNR[n])) {
3053  NRdSmax = dSNR[n];
3054  index = n;
3055  }
3056  }
3057  index /= numPhase;
3058  return NRdSmax;
3059 }
3060 
3061 OCP_DBL Bulk::GetNRdPmax() { return NRdPmax; }
3062 
3063 OCP_DBL Bulk::GetNRdSmaxP() { return NRdSmaxP; }
3064 
3065 OCP_DBL Bulk::GetNRdNmax() { return NRdNmax; }
3066 
3067 void Bulk::CorrectNi(const vector<OCP_DBL>& res)
3068 {
3069  for (OCP_USI n = 0; n < numBulk; n++) {
3070  for (USI i = 0; i < numCom; i++) {
3071  Ni[n * numCom + i] += res[n * (numCom + 1) + i];
3072  if (Ni[n * numCom + i] < 0) {
3073  Ni[n * numCom + i] = 1E-8 * Nt[n];
3074  }
3075  }
3076  }
3077 }
3078 
3079 void Bulk::OutputInfo(const OCP_USI& n) const
3080 {
3081  OCP_USI bIdC = n * numCom;
3082  OCP_USI bIdP = n * numPhase;
3083  OCP_USI bIdPC = bIdP * numCom;
3084 
3085  cout << "------------------------------" << endl;
3086  cout << "Bulk[" << n << "]" << endl;
3087  cout << fixed << setprecision(3);
3088  for (USI i = 0; i < numCom; i++) {
3089  cout << Ni[bIdC + i] << " ";
3090  }
3091  cout << endl << P[n] << " " << T;
3092  cout << endl;
3093 
3094  cout << fixed << setprecision(8);
3095  if (phaseExist[bIdP + 0]) {
3096  for (USI i = 0; i < numCom; i++) {
3097  cout << xij[bIdPC + i] << " ";
3098  }
3099  } else {
3100  for (USI i = 0; i < numCom; i++) {
3101  cout << 0.000000 << " ";
3102  }
3103  }
3104  cout << phaseExist[bIdP + 0] << " ";
3105  cout << xi[bIdP + 0] << " ";
3106  cout << S[bIdP + 0] << " ";
3107  cout << endl;
3108 
3109  if (phaseExist[bIdP + 1]) {
3110  for (USI i = 0; i < numCom; i++) {
3111  cout << xij[bIdPC + numCom + i] << " ";
3112  }
3113  } else {
3114  for (USI i = 0; i < numCom; i++) {
3115  cout << 0.000000 << " ";
3116  }
3117  }
3118 
3119  cout << phaseExist[bIdP + 1] << " ";
3120  cout << xi[bIdP + 1] << " ";
3121  cout << S[bIdP + 1] << " ";
3122  cout << endl;
3123 
3124  cout << fixed << setprecision(3);
3125 
3126  cout << vf[n] << " " << rockVp[n] << " ";
3127  cout << fixed << setprecision(12);
3128  cout << fabs(vf[n] - rockVp[n]) / rockVp[n] << endl;
3129  cout << fixed << setprecision(3);
3130  cout << vj[bIdP] << " " << vj[bIdP + 1] << " " << vj[bIdP + 2] << endl;
3131  cout << "------------------------------" << endl;
3132 }
3133 
3135 // For AIMt
3137 
3138 void Bulk::AllocateAuxAIM(const OCP_DBL& ratio)
3139 {
3140  OCP_FUNCNAME;
3141 
3142  maxNumFIMBulk = numBulk * ratio;
3143  if (maxNumFIMBulk < wellBulkId.capacity()) {
3144  maxNumFIMBulk = wellBulkId.capacity();
3145  }
3146  FIMBulk.reserve(maxNumFIMBulk);
3147  map_Bulk2FIM.resize(numBulk, -1);
3148 
3149  muP.resize(maxNumFIMBulk * numPhase);
3150  xiP.resize(maxNumFIMBulk * numPhase);
3151  rhoP.resize(maxNumFIMBulk * numPhase);
3152  mux.resize(maxNumFIMBulk * numCom * numPhase);
3153  xix.resize(maxNumFIMBulk * numCom * numPhase);
3154  rhox.resize(maxNumFIMBulk * numCom * numPhase);
3155  lendSdP = (numCom + 1) * (numCom + 1) * numPhase;
3156  dSec_dPri.resize(maxNumFIMBulk * lendSdP);
3157  dKr_dS.resize(maxNumFIMBulk * numPhase * numPhase);
3158  dPcj_dS.resize(maxNumFIMBulk * numPhase * numPhase);
3159 
3160  lmuP.resize(maxNumFIMBulk * numPhase);
3161  lxiP.resize(maxNumFIMBulk * numPhase);
3162  lrhoP.resize(maxNumFIMBulk * numPhase);
3163  lmux.resize(maxNumFIMBulk * numCom * numPhase);
3164  lxix.resize(maxNumFIMBulk * numCom * numPhase);
3165  lrhox.resize(maxNumFIMBulk * numCom * numPhase);
3166  ldSec_dPri.resize(maxNumFIMBulk * lendSdP);
3167  ldKr_dS.resize(maxNumFIMBulk * numPhase * numPhase);
3168  ldPcj_dS.resize(maxNumFIMBulk * numPhase * numPhase);
3169 
3170  FIMNi.resize(maxNumFIMBulk * numCom);
3171 }
3172 
3173 void Bulk::FlashDerivAIM(const bool& IfAIMs)
3174 {
3175  OCP_FUNCNAME;
3176 
3177  USI ftype;
3178  OCP_USI n, bIdc;
3179  OCP_DBL Ntw;
3180  OCP_DBL minEig;
3181  dSec_dPri.clear();
3182 
3183  USI len;
3184  if (IfAIMs) {
3185  len = FIMBulk.size();
3186  } else {
3187  len = numFIMBulk;
3188  }
3189 
3190  for (OCP_USI k = 0; k < len; k++) {
3191  n = FIMBulk[k];
3192  bIdc = n * numCom;
3193  ftype = 1;
3194  if (flagSkip[n]) {
3195  minEig = minEigenSkip[n];
3196  if (fabs(1 - PSkip[n] / P[n]) >= minEig / 10) {
3197  ftype = 0;
3198  }
3199  if (ftype == 1) {
3200 
3201  Ntw = Nt[n] - Ni[bIdc + numCom_1];
3202  for (USI i = 0; i < numCom_1; i++) {
3203  // cout << fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) << " ";
3204  if (fabs(Ni[bIdc + i] / Ntw - ziSkip[bIdc + i]) >= minEig / 10) {
3205  ftype = 0;
3206  break;
3207  }
3208  }
3209  }
3210  // cout << n << endl;
3211  } else {
3212  ftype = 0;
3213  }
3214 
3215  flashCal[PVTNUM[n]]->FlashDeriv(P[n], T, &Ni[n * numCom], ftype, phaseNum[n],
3216  &Ks[n * numCom_1]);
3217  PassFlashValueDerivAIM(k);
3218  }
3219 
3220 #ifdef DEBUG
3221  CheckSat();
3222 #endif // DEBUG
3223 }
3224 
3225 void Bulk::PassFlashValueDerivAIM(const OCP_USI& fn)
3226 {
3227  OCP_FUNCNAME;
3228 
3229  OCP_USI n = FIMBulk[fn];
3230  OCP_USI bIdp = n * numPhase;
3231  OCP_USI fnp = fn * numPhase;
3232  USI pvtnum = PVTNUM[n];
3233  USI nptmp = 0;
3234  for (USI j = 0; j < numPhase; j++) {
3235  phaseExist[bIdp + j] = flashCal[pvtnum]->phaseExist[j];
3236  // Important! Saturation must be passed no matter if the phase exists. This is
3237  // because it will be used to calculate relative permeability and capillary
3238  // pressure at each time step. Make sure that all saturations are updated at
3239  // each step!
3240  S[bIdp + j] = flashCal[pvtnum]->S[j];
3241  if (phaseExist[bIdp + j]) { // j -> bId + j fix bugs.
3242  nptmp++;
3243  rho[bIdp + j] = flashCal[pvtnum]->rho[j];
3244  xi[bIdp + j] = flashCal[pvtnum]->xi[j];
3245  mu[bIdp + j] = flashCal[pvtnum]->mu[j];
3246  vj[bIdp + j] = flashCal[pvtnum]->v[j];
3247 
3248  // Derivatives
3249  muP[fnp + j] = flashCal[pvtnum]->muP[j];
3250  xiP[fnp + j] = flashCal[pvtnum]->xiP[j];
3251  rhoP[fnp + j] = flashCal[pvtnum]->rhoP[j];
3252  for (USI i = 0; i < numCom; i++) {
3253  xij[bIdp * numCom + j * numCom + i] =
3254  flashCal[pvtnum]->xij[j * numCom + i];
3255  mux[fnp * numCom + j * numCom + i] =
3256  flashCal[pvtnum]->mux[j * numCom + i];
3257  xix[fnp * numCom + j * numCom + i] =
3258  flashCal[pvtnum]->xix[j * numCom + i];
3259  rhox[fnp * numCom + j * numCom + i] =
3260  flashCal[pvtnum]->rhox[j * numCom + i];
3261  }
3262  }
3263  }
3264  Nt[n] = flashCal[pvtnum]->Nt;
3265  vf[n] = flashCal[pvtnum]->vf;
3266  vfp[n] = flashCal[pvtnum]->vfp;
3267 
3268  OCP_USI bIdc = n * numCom;
3269  for (USI i = 0; i < numCom; i++) {
3270  vfi[bIdc + i] = flashCal[pvtnum]->vfi[i];
3271  }
3272  dSec_dPri.insert(dSec_dPri.end(), flashCal[pvtnum]->dXsdXp.begin(),
3273  flashCal[pvtnum]->dXsdXp.end());
3274 
3275  phaseNum[n] = nptmp - 1; // water is excluded
3276 
3277  if (comps) {
3278  if (nptmp == 3) {
3279  // num of hydrocarbon phase equals 2
3280  // Calculate Ks
3281  OCP_USI bIdc1 = n * numCom_1;
3282  for (USI i = 0; i < numCom_1; i++) {
3283  Ks[bIdc1 + i] =
3284  flashCal[pvtnum]->xij[i] / flashCal[pvtnum]->xij[numCom + i];
3285  }
3286  }
3287 
3288  if (flashCal[pvtnum]->GetFtype() == 0) {
3289  flagSkip[n] = flashCal[pvtnum]->GetFlagSkip();
3290  if (flagSkip[n]) {
3291  minEigenSkip[n] = flashCal[pvtnum]->GetMinEigenSkip();
3292  for (USI j = 0; j < numPhase - 1; j++) {
3293  if (phaseExist[bIdp + j]) {
3294  for (USI i = 0; i < numCom - 1; i++) {
3295  ziSkip[bIdc + i] = flashCal[pvtnum]->xij[j * numCom + i];
3296  }
3297  break;
3298  }
3299  }
3300  PSkip[n] = P[n];
3301  }
3302  }
3303  if (miscible) {
3304  surTen[n] = flashCal[pvtnum]->GetSurTen();
3305  }
3306  }
3307 }
3308 
3309 void Bulk::CalKrPcDerivAIM(const bool& IfAIMs)
3310 {
3311  OCP_FUNCNAME;
3312 
3313  USI len;
3314  if (IfAIMs) {
3315  len = FIMBulk.size();
3316  } else {
3317  len = numFIMBulk;
3318  }
3319  OCP_DBL tmp;
3320  for (OCP_USI fn = 0; fn < len; fn++) {
3321  OCP_USI n = FIMBulk[fn];
3322  OCP_USI fnp = fn * numPhase;
3323  OCP_USI bId = n * numPhase;
3324  flow[SATNUM[n]]->CalKrPcDeriv(&S[bId], &kr[bId], &Pc[bId],
3325  &dKr_dS[fnp * numPhase], &dPcj_dS[fnp * numPhase],
3326  0, tmp, tmp);
3327  for (USI j = 0; j < numPhase; j++) Pj[bId + j] = P[n] + Pc[bId + j];
3328  }
3329 }
3330 
3331 void Bulk::CalRelResAIMt(ResFIM& resFIM) const
3332 {
3333  OCP_FUNCNAME;
3334 
3335  // OCP_USI tmpid01 = -1;
3336  // OCP_USI tmpid02 = -1;
3337  OCP_DBL tmp;
3338 
3339  const USI len = numCom + 1;
3340  for (OCP_USI fn = 0; fn < numFIMBulk; fn++) {
3341  OCP_USI n = FIMBulk[fn];
3342 
3343  for (USI i = 0; i < len; i++) {
3344  tmp = fabs(resFIM.res[fn * len + i] / rockVp[n]);
3345  if (resFIM.maxRelRes_v < tmp) {
3346  resFIM.maxRelRes_v = tmp;
3347  // tmpid01 = n;
3348  }
3349  }
3350  for (USI i = 1; i < len; i++) {
3351  tmp = fabs(resFIM.res[fn * len + i] / Nt[n]);
3352  if (resFIM.maxRelRes_mol < tmp) {
3353  resFIM.maxRelRes_mol = tmp;
3354  // tmpid02 = n;
3355  }
3356  }
3357  }
3358 }
3359 
3360 void Bulk::GetSolAIMt(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
3361  const OCP_DBL& dSmaxlim)
3362 {
3363  OCP_FUNCNAME;
3364 
3365  NRdSmaxP = 0;
3366  NRdPmax = 0;
3367  OCP_DBL dP;
3368  USI row = numPhase * (numCom + 1);
3369  USI col = numCom + 1;
3370  USI bsize = row * col;
3371  vector<OCP_DBL> dtmp(row, 0);
3372  OCP_DBL chopmin = 1;
3373  OCP_DBL choptmp = 0;
3374 
3375  for (OCP_USI fn = 0; fn < numFIMBulk; fn++) {
3376  OCP_USI n = FIMBulk[fn];
3377 
3378  chopmin = 1;
3379  // compute the chop
3380  fill(dtmp.begin(), dtmp.end(), 0.0);
3381  DaAxpby(row, col, 1, dSec_dPri.data() + fn * bsize, u.data() + fn * col, 1,
3382  dtmp.data());
3383 
3384  for (USI j = 0; j < numPhase; j++) {
3385 
3386  choptmp = 1;
3387  if (fabs(dtmp[j]) > dSmaxlim) {
3388  choptmp = dSmaxlim / fabs(dtmp[j]);
3389  } else if (S[n * numPhase + j] + dtmp[j] < 0.0) {
3390  choptmp = 0.9 * S[n * numPhase + j] / fabs(dtmp[j]);
3391  }
3392 
3393  chopmin = min(chopmin, choptmp);
3394  NRdSmaxP = max(NRdSmaxP, choptmp * fabs(dtmp[j]));
3395  }
3396  dP = u[fn * col];
3397  choptmp = dPmaxlim / fabs(dP);
3398  chopmin = min(chopmin, choptmp);
3399  dP *= chopmin;
3400  NRdPmax = max(NRdPmax, fabs(dP));
3401  P[n] += dP; // seems better
3402 
3404  // for (USI i = 0; i < numCom; i++) {
3405  // if (Ni[n * numCom + i] + u[n * ncol + 1 + i] < 0) {
3406  // chopmin = 0.9 * min(chopmin, fabs(Ni[n * numCom + i] / u[n * ncol + 1
3407  // + i]));
3408 
3409  // //if (chopmin < 0 || !isfinite(chopmin)) {
3410  // // OCP_ABORT("Wrong Chop!");
3411  // //}
3412  // }
3413  //}
3414 
3415  for (USI i = 0; i < numCom; i++) {
3416  Ni[n * numCom + i] += u[fn * col + 1 + i] * chopmin;
3417 
3418  // if (Ni[n * numCom + i] < 0) {
3419  // cout << Ni[n * numCom + i] << " " << u[n * ncol + 1 + i] * chopmin <<
3420  // " " << chopmin << endl;
3421  //}
3422  }
3423  }
3424 }
3425 
3426 void Bulk::CalRelResAIMs(ResFIM& resFIM) const
3427 {
3428  OCP_FUNCNAME;
3429 
3430  // OCP_USI tmpid01 = -1;
3431  // OCP_USI tmpid02 = -1;
3432  OCP_DBL tmp;
3433 
3434  const USI len = numCom + 1;
3435  for (OCP_USI fn = 0; fn < numFIMBulk; fn++) {
3436  OCP_USI n = FIMBulk[fn];
3437 
3438  for (USI i = 0; i < len; i++) {
3439  tmp = fabs(resFIM.res[n * len + i] / rockVp[n]);
3440  if (resFIM.maxRelRes_v < tmp) {
3441  resFIM.maxRelRes_v = tmp;
3442  // tmpid01 = n;
3443  }
3444  }
3445  for (USI i = 1; i < len; i++) {
3446  tmp = fabs(resFIM.res[n * len + i] / Nt[n]);
3447  if (resFIM.maxRelRes_mol < tmp) {
3448  resFIM.maxRelRes_mol = tmp;
3449  // tmpid02 = n;
3450  }
3451  }
3452  }
3453 }
3454 
3455 void Bulk::GetSolAIMs(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
3456  const OCP_DBL& dSmaxlim)
3457 {
3458  OCP_FUNCNAME;
3459 
3460  NRdSmaxP = 0;
3461  NRdPmax = 0;
3462  OCP_DBL dP;
3463  USI row = numPhase * (numCom + 1);
3464  USI col = numCom + 1;
3465  USI bsize = row * col;
3466  vector<OCP_DBL> dtmp(row, 0);
3467  OCP_DBL chopmin = 1;
3468  OCP_DBL choptmp = 0;
3469  OCP_INT bIde;
3470 
3471  for (OCP_USI n = 0; n < numBulk; n++) {
3472  bIde = map_Bulk2FIM[n];
3473 
3474  if (bIde < 0 || bIde >= numFIMBulk) {
3475  // IMPEC Bulk
3476  // Method 1
3477  P[n] += u[n * col];
3478  // Method 2
3479  // P[n] = u[n * ncol];
3480  if (P[n] < 0) {
3481  cout << "";
3482  }
3483  } else {
3484  // FIM Bulk
3485  chopmin = 1;
3486  // compute the chop
3487  fill(dtmp.begin(), dtmp.end(), 0.0);
3488  DaAxpby(row, col, 1, dSec_dPri.data() + bIde * bsize, u.data() + bIde * col,
3489  1, dtmp.data());
3490 
3491  for (USI j = 0; j < numPhase; j++) {
3492  choptmp = 1;
3493  if (fabs(dtmp[j]) > dSmaxlim) {
3494  choptmp = dSmaxlim / fabs(dtmp[j]);
3495  } else if (S[n * numPhase + j] + dtmp[j] < 0.0) {
3496  choptmp = 0.9 * S[n * numPhase + j] / fabs(dtmp[j]);
3497  }
3498 
3499  chopmin = min(chopmin, choptmp);
3500  NRdSmaxP = max(NRdSmaxP, choptmp * fabs(dtmp[j]));
3501  }
3502  dP = u[n * col];
3503  choptmp = dPmaxlim / fabs(dP);
3504  chopmin = min(chopmin, choptmp);
3505  dP *= chopmin;
3506  NRdPmax = max(NRdPmax, fabs(dP));
3507  P[n] += dP; // seems better
3508 
3510  // for (USI i = 0; i < numCom; i++) {
3511  // if (Ni[n * numCom + i] + u[n * ncol + 1 + i] < 0) {
3512  // chopmin = 0.9 * min(chopmin, fabs(Ni[n * numCom + i] / u[n * ncol
3513  // + 1 + i]));
3514 
3515  // //if (chopmin < 0 || !isfinite(chopmin)) {
3516  // // OCP_ABORT("Wrong Chop!");
3517  // //}
3518  // }
3519  //}
3520 
3521  for (USI i = 0; i < numCom; i++) {
3522  Ni[n * numCom + i] += u[n * col + 1 + i] * chopmin;
3523 
3524  // if (Ni[n * numCom + i] < 0) {
3525  // cout << Ni[n * numCom + i] << " " << u[n * ncol + 1 + i] *
3526  // chopmin << " " << chopmin << endl;
3527  //}
3528  }
3529  }
3530  }
3531 }
3532 
3533 void Bulk::UpdateLastStepAIM()
3534 {
3535  lmuP = muP;
3536  lxiP = xiP;
3537  lrhoP = rhoP;
3538  lmux = mux;
3539  lxix = xix;
3540  lrhox = rhox;
3541  ldPcj_dS = dPcj_dS;
3542  ldKr_dS = dKr_dS;
3543  ldSec_dPri = dSec_dPri;
3544 }
3545 
3546 void Bulk::ResetFIMBulk()
3547 {
3548  muP = lmuP;
3549  xiP = lxiP;
3550  rhoP = lrhoP;
3551  mux = lmux;
3552  xix = lxix;
3553  rhox = lrhox;
3554  dPcj_dS = ldPcj_dS;
3555  dKr_dS = ldKr_dS;
3556  dSec_dPri = ldSec_dPri;
3557 }
3558 
3559 void Bulk::ShowFIMBulk(const bool& flag) const
3560 {
3561  cout << numFIMBulk << " " << fixed << setprecision(3)
3562  << numFIMBulk * 100.0 / numBulk << "%" << endl;
3563  if (flag) {
3564  for (USI n = 0; n < numFIMBulk; n++) {
3565  cout << setw(6) << FIMBulk[n] << " ";
3566  if ((n + 1) % 10 == 0) {
3567  cout << endl;
3568  }
3569  }
3570  cout << endl;
3571  }
3572  if (false) {
3573  for (USI n = 0; n < numBulk; n++) {
3574  cout << setw(6) << map_Bulk2FIM[n] << " ";
3575  if ((n + 1) % 10 == 0) {
3576  cout << endl;
3577  }
3578  }
3579  cout << endl;
3580  }
3581 }
3582 
3583 bool Bulk::CheckNiFIMBulk() const
3584 {
3585  OCP_FUNCNAME;
3586 
3587  OCP_USI len = numBulk * numCom;
3588  for (OCP_USI n = 0; n < len; n++) {
3589  if (Ni[n] < 0.0) {
3590  OCP_USI bId = n / numCom;
3591  USI cId = n - bId * numCom;
3592  std::ostringstream NiStringSci;
3593  NiStringSci << std::scientific << Ni[n];
3594  OCP_WARNING("Negative Ni: Ni[" + std::to_string(cId) + "] in Bulk[" +
3595  std::to_string(bId) + "] = " + NiStringSci.str());
3596  return false;
3597  }
3598  }
3599  return true;
3600 }
3601 
3602 void Bulk::InFIMNi()
3603 {
3604  OCP_USI bIdb, bIdf;
3605  for (OCP_USI n = 0; n < numFIMBulk; n++) {
3606  bIdf = n * numCom;
3607  bIdb = FIMBulk[n] * numCom;
3608  for (USI i = 0; i < numCom; i++) {
3609  FIMNi[bIdf + i] = Ni[bIdb + i];
3610  }
3611  }
3612 }
3613 
3614 void Bulk::OutFIMNi()
3615 {
3616  OCP_USI bIdb, bIdf;
3617  for (OCP_USI n = 0; n < numFIMBulk; n++) {
3618  bIdf = n * numCom;
3619  bIdb = FIMBulk[n] * numCom;
3620  for (USI i = 0; i < numCom; i++) {
3621  Ni[bIdb + i] = FIMNi[bIdf + i];
3622  }
3623  }
3624 }
3625 
3626 void Bulk::AllocateAuxAIMc()
3627 {
3628  cfl.resize(numBulk * numPhase);
3629  map_Bulk2FIM.resize(numBulk, -1);
3630 }
3631 
3632 void Bulk::FlashAIMc()
3633 {
3634  if (comps) {
3635  FlashCOMPAIMc();
3636  } else {
3637  FlashBLKOILAIMc();
3638  }
3639 }
3640 
3641 void Bulk::FlashBLKOILAIMc()
3642 {
3643  for (OCP_USI n = 0; n < numBulk; n++) {
3644  if (map_Bulk2FIM[n] > -1) {
3645  // FIM bulk
3646  continue;
3647  }
3648 
3649  flashCal[PVTNUM[n]]->Flash(P[n], T, &Ni[n * numCom], 0, 0, 0);
3650  PassFlashValueAIMc(n);
3651  }
3652 }
3653 
3654 void Bulk::FlashCOMPAIMc()
3655 {
3656  USI ftype;
3657  OCP_USI bId;
3658  OCP_DBL Ntw;
3659  OCP_DBL minEig;
3660  // cout << endl << "==================================" << endl;
3661  for (OCP_USI n = 0; n < numBulk; n++) {
3662  if (map_Bulk2FIM[n] > -1) {
3663  // FIM bulk
3664  continue;
3665  }
3666 
3667  ftype = 1;
3668  if (flagSkip[n]) {
3669  minEig = minEigenSkip[n];
3670  if (fabs(1 - PSkip[n] / P[n]) >= minEig / 10) {
3671  ftype = 0;
3672  }
3673  // cout << setprecision(2) << scientific << minEig / 10 << " " << fabs(1 -
3674  // lP[n] / P[n]) << " ";
3675  if (ftype == 1) {
3676  bId = n * numCom;
3677  Ntw = Nt[n] - Ni[bId + numCom - 1];
3678  for (USI i = 0; i < numCom - 1; i++) {
3679  // cout << fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) << " ";
3680  if (fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) >= minEig / 10) {
3681  ftype = 0;
3682  break;
3683  }
3684  }
3685  }
3686  // cout << n << endl;
3687  } else {
3688  ftype = 0;
3689  }
3690  flashCal[PVTNUM[n]]->Flash(P[n], T, &Ni[n * numCom], ftype, phaseNum[n],
3691  &Ks[n * numCom_1]);
3692  PassFlashValueAIMc(n);
3693  }
3694 }
3695 
3696 void Bulk::FlashAIMc01()
3697 {
3698  if (comps) {
3699  FlashCOMPAIMc01();
3700  } else {
3701  FlashBLKOILAIMc01();
3702  }
3703 }
3704 
3705 void Bulk::FlashBLKOILAIMc01()
3706 {
3707  for (OCP_USI n = 0; n < numBulk; n++) {
3708  if (map_Bulk2FIM[n] > -1) {
3709  // FIM bulk
3710  continue;
3711  }
3712 
3713  flashCal[PVTNUM[n]]->Flash(P[n], T, &Ni[n * numCom], 0, 0, 0);
3714  PassFlashValue(n);
3715  }
3716 }
3717 
3718 void Bulk::FlashCOMPAIMc01()
3719 {
3720  USI ftype;
3721  OCP_USI bId;
3722  OCP_DBL Ntw;
3723  OCP_DBL minEig;
3724  // cout << endl << "==================================" << endl;
3725  for (OCP_USI n = 0; n < numBulk; n++) {
3726  if (map_Bulk2FIM[n] > -1) {
3727  // FIM bulk
3728  continue;
3729  }
3730 
3731  ftype = 1;
3732  if (flagSkip[n]) {
3733  minEig = minEigenSkip[n];
3734  if (fabs(1 - PSkip[n] / P[n]) >= minEig / 10) {
3735  ftype = 0;
3736  }
3737  // cout << setprecision(2) << scientific << minEig / 10 << " " << fabs(1 -
3738  // lP[n] / P[n]) << " ";
3739  if (ftype == 1) {
3740  bId = n * numCom;
3741  Ntw = Nt[n] - Ni[bId + numCom - 1];
3742  for (USI i = 0; i < numCom - 1; i++) {
3743  // cout << fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) << " ";
3744  if (fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) >= minEig / 10) {
3745  ftype = 0;
3746  break;
3747  }
3748  }
3749  }
3750  // cout << n << endl;
3751  } else {
3752  ftype = 0;
3753  }
3754  flashCal[PVTNUM[n]]->Flash(P[n], T, &Ni[n * numCom], ftype, phaseNum[n],
3755  &Ks[n * numCom_1]);
3756  PassFlashValue(n);
3757  }
3758 }
3759 
3760 void Bulk::FlashDerivAIMc()
3761 {
3762  if (comps) {
3763  FlashDerivCOMPAIMc();
3764  } else {
3765  FlashDerivBLKOILAIMc();
3766  }
3767 }
3768 
3769 void Bulk::FlashDerivBLKOILAIMc()
3770 {
3771  // dSec_dPri.clear();
3772  for (auto& n : FIMBulk) {
3773  flashCal[PVTNUM[n]]->FlashDeriv(P[n], T, &Ni[n * numCom], 0, 0, 0);
3774  PassFlashValueDeriv(n);
3775  }
3776 }
3777 
3778 void Bulk::FlashDerivCOMPAIMc()
3779 {
3780  USI ftype;
3781  OCP_USI bId;
3782  OCP_DBL Ntw;
3783  OCP_DBL minEig;
3784  // dSec_dPri.clear();
3785  for (auto& n : FIMBulk) {
3786  ftype = 1;
3787  if (flagSkip[n]) {
3788  minEig = minEigenSkip[n];
3789  if (fabs(1 - PSkip[n] / P[n]) >= minEig / 10) {
3790  ftype = 0;
3791  }
3792  // cout << setprecision(2) << scientific << minEig / 10 << " " << fabs(1 -
3793  // lP[n] / P[n]) << " ";
3794  if (ftype == 1) {
3795  bId = n * numCom;
3796  Ntw = Nt[n] - Ni[bId + numCom - 1];
3797  for (USI i = 0; i < numCom - 1; i++) {
3798  // cout << fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) << " ";
3799  if (fabs(Ni[bId + i] / Ntw - ziSkip[bId + i]) >= minEig / 10) {
3800  ftype = 0;
3801  break;
3802  }
3803  }
3804  }
3805  // cout << n << endl;
3806  } else {
3807  ftype = 0;
3808  }
3809 
3810  flashCal[PVTNUM[n]]->FlashDeriv(P[n], T, &Ni[n * numCom], ftype, phaseNum[n],
3811  &Ks[n * numCom_1]);
3812  PassFlashValueDeriv(n);
3813  // lphaseNum[n];
3814  }
3815 }
3816 
3817 void Bulk::CalKrPcAIMc()
3818 {
3819  OCP_DBL tmp = 0;
3820  for (OCP_USI n = 0; n < numBulk; n++) {
3821  if (map_Bulk2FIM[n] > -1) {
3822  // FIM bulk
3823  continue;
3824  }
3825  OCP_USI bId = n * numPhase;
3826  flow[SATNUM[n]]->CalKrPc(&S[bId], &kr[bId], &Pc[bId], 0, tmp, tmp);
3827  for (USI j = 0; j < numPhase; j++)
3828  Pj[n * numPhase + j] = P[n] + Pc[n * numPhase + j];
3829  }
3830 }
3831 
3832 void Bulk::CalKrPcDerivAIMc()
3833 {
3834  OCP_DBL tmp;
3835  for (auto& n : FIMBulk) {
3836  OCP_USI bId = n * numPhase;
3837  flow[SATNUM[n]]->CalKrPcDeriv(&S[bId], &kr[bId], &Pc[bId],
3838  &dKr_dS[bId * numPhase], &dPcj_dS[bId * numPhase],
3839  0, tmp, tmp);
3840  for (USI j = 0; j < numPhase; j++) Pj[bId + j] = P[n] + Pc[bId + j];
3841  }
3842 }
3843 
3844 void Bulk::GetSolAIMc(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
3845  const OCP_DBL& dSmaxlim)
3846 {
3847  NRdPmax = 0;
3848  NRdNmax = 0;
3849  OCP_DBL tmp;
3850  USI row = numPhase * (numCom + 1);
3851  USI col = numCom + 1;
3852  OCP_DBL chopmin = 1;
3853  OCP_DBL choptmp = 0;
3854 
3855  for (OCP_USI n = 0; n < numBulk; n++) {
3856  // cout << "[" << n << "]" << endl;
3857  // cout << scientific << P[n] << " " << u[n * ncol] << " " << u[n * ncol] /
3858  // P[n] << endl;
3859 
3860  dPNR[n] = u[n * col];
3861  tmp = fabs(dPNR[n] / P[n]);
3862  P[n] += dPNR[n];
3863 
3864  if (tmp > NRdPmax) {
3865  NRdPmax = tmp;
3866  }
3867 
3868  for (USI i = 0; i < numCom; i++) {
3869  /*cout << Ni[n * numCom + i] << " " << u[n * ncol + 1 + i] << " "
3870  << u[n * ncol + 1 + i] / Ni[n * numCom + i] << endl;*/
3871 
3872  if (Ni[n * numCom + i] + u[n * col + 1 + i] < 0) {
3873  dNNR[n * numCom + i] = -Ni[n * numCom + i];
3874  tmp = 1;
3875  Ni[n * numCom + i] = 1E-8 * Nt[n];
3876  // cout << Ni[n * numCom + i] << " " << u[n * ncol + 1 + i] << endl;
3877  } else if (u[n * col + 1 + i] > Ni[n * numCom + i]) {
3878  dNNR[n * numCom + i] = Ni[n * numCom + i];
3879  tmp = 1;
3880  Ni[n * numCom + i] += dNNR[n * numCom + i];
3881  } else {
3882  dNNR[n * numCom + i] = u[n * col + 1 + i];
3883  tmp = fabs(dNNR[n * numCom + i]) / Ni[n * numCom + i];
3884  Ni[n * numCom + i] += dNNR[n * numCom + i];
3885  }
3886  if (tmp > NRdNmax) {
3887  NRdNmax = tmp;
3888  }
3889  }
3890  }
3891 }
3892 
3893 void Bulk::GetSolAIMc01(const vector<OCP_DBL>& u, const OCP_DBL& dPmaxlim,
3894  const OCP_DBL& dSmaxlim)
3895 {
3896  NRdSmaxP = 0;
3897  NRdPmax = 0;
3898  OCP_DBL dP;
3899  USI row0 = numPhase;
3900  USI row = numPhase * (numCom + 1);
3901  USI col = numCom + 1;
3902  USI bsize = row * col;
3903  vector<OCP_DBL> dtmp(row0, 0);
3904  OCP_DBL chopmin = 1;
3905  OCP_DBL choptmp = 0;
3906 
3907  for (OCP_USI n = 0; n < numBulk; n++) {
3908  if (map_Bulk2FIM[n] < 0) {
3909  // IMPEC Bulk
3910  // Pressure
3911  dP = u[n * col];
3912  NRdPmax = max(NRdPmax, fabs(dP));
3913  P[n] += dP; // seems better
3914  dPNR[n] = dP;
3915  NRstep[n] = 1;
3916  // Ni
3917  for (USI i = 0; i < numCom; i++) {
3918  dNNR[n * numCom + i] = u[n * col + 1 + i];
3919  Ni[n * numCom + i] += dNNR[n * numCom + i];
3920  }
3921  continue;
3922  }
3923 
3924  chopmin = 1;
3925  // compute the chop
3926  fill(dtmp.begin(), dtmp.end(), 0.0);
3927  DaAxpby(row0, col, 1, dSec_dPri.data() + n * bsize, u.data() + n * col, 1,
3928  dtmp.data());
3929 
3930  for (USI j = 0; j < numPhase; j++) {
3931  choptmp = 1;
3932  if (fabs(dtmp[j]) > dSmaxlim) {
3933  choptmp = dSmaxlim / fabs(dtmp[j]);
3934  } else if (S[n * numPhase + j] + dtmp[j] < 0.0) {
3935  choptmp = 0.9 * S[n * numPhase + j] / fabs(dtmp[j]);
3936  }
3937 
3938  chopmin = min(chopmin, choptmp);
3939  NRdSmaxP = max(NRdSmaxP, choptmp * fabs(dtmp[j]));
3940  }
3941  dP = u[n * col];
3942  choptmp = dPmaxlim / fabs(dP);
3943  chopmin = min(chopmin, choptmp);
3944  NRdPmax = max(NRdPmax, fabs(dP));
3945  P[n] += dP; // seems better
3946  dPNR[n] = dP;
3947  NRstep[n] = chopmin;
3948 
3949  for (USI i = 0; i < numCom; i++) {
3950  dNNR[n * numCom + i] = u[n * col + 1 + i] * chopmin;
3951  Ni[n * numCom + i] += dNNR[n * numCom + i];
3952  }
3953  }
3954 }
3955 
3956 void Bulk::UpdatePj()
3957 {
3958  for (OCP_USI n = 0; n < numBulk; n++) {
3959  for (USI j = 0; j < numPhase; j++) {
3960  Pj[n * numPhase + j] = P[n] + Pc[n * numPhase + j];
3961  }
3962  }
3963 }
3964 
3965 /*----------------------------------------------------------------------------*/
3966 /* Brief Change History of This File */
3967 /*----------------------------------------------------------------------------*/
3968 /* Author Date Actions */
3969 /*----------------------------------------------------------------------------*/
3970 /* Shizhe Li Oct/01/2021 Create file */
3971 /* Chensong Zhang Oct/15/2021 Format file */
3972 /* Chensong Zhang Jan/09/2022 Update Doxygen */
3973 /*----------------------------------------------------------------------------*/
Bulk.hpp
Bulk class declaration.
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
BLKOIL
const USI BLKOIL
Mixture model = black-oil.
Definition: OCPConst.hpp:88
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:36
PHASE_ODGW01
const USI PHASE_ODGW01
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:102
EOS_PVTW
const USI EOS_PVTW
Mixture model = equation-of-state.
Definition: OCPConst.hpp:89
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
GAS
const USI GAS
Fluid type = gas.
Definition: OCPConst.hpp:83
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
CONV1
const OCP_DBL CONV1
1 bbl = CONV1 ft3
Definition: OCPConst.hpp:59
WATER
const USI WATER
Fluid type = water.
Definition: OCPConst.hpp:84
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
OIL
const USI OIL
Fluid type = oil.
Definition: OCPConst.hpp:82
PHASE_W
const USI PHASE_W
Phase type = water only.
Definition: OCPConst.hpp:96
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:24
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
OCP_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:25
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
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:191
BOMixture_OW
Definition: MixtureBO.hpp:140
Bulk::AllocateAuxAIM
void AllocateAuxAIM(const OCP_DBL &ratio)
Allocate memory for auxiliary variables used for AIMt.
Definition: Bulk.cpp:3138
Bulk::InitSjPcComp
void InitSjPcComp(const USI &tabrow, const Grid &myGrid)
Calculate initial equilibrium for compositional model according to EQUIL.
Definition: Bulk.cpp:812
Bulk::PassFlashValue
void PassFlashValue(const OCP_USI &n)
Pass values from Flash to Bulk after Flash calculation.
Definition: Bulk.cpp:1529
Bulk::FlashBLKOIL
void FlashBLKOIL()
Perform flash calculation with Ni in Black Oil Model.
Definition: Bulk.cpp:1355
Bulk::CalFlashType
USI CalFlashType(const OCP_USI &n) const
determine which flash type will be used
Definition: Bulk.cpp:1480
Bulk::CalRelResAIMs
void CalRelResAIMs(ResFIM &resFIM) const
Calculate relative resiual for AIMs, parts related to FIM are considered.
Definition: Bulk.cpp:3426
Bulk::GetSolIMPEC
void GetSolIMPEC(const vector< OCP_DBL > &u)
Update P and Pj after linear system is solved.
Definition: Bulk.cpp:2446
Bulk::GetNRdPmax
OCP_DBL GetNRdPmax()
Return NRdPmax.
Definition: Bulk.cpp:3061
Bulk::CheckDiff
void CheckDiff()
Check difference from last time step, for Debug and Test.
Definition: Bulk.cpp:2182
Bulk::FlashCOMPAIMc
void FlashCOMPAIMc()
Perform flash calculation with Ni in Compositional Model.
Definition: Bulk.cpp:3654
Bulk::InitFlash
void InitFlash(const bool &flag=false)
Perform flash calculation with saturations.
Definition: Bulk.cpp:1280
Bulk::FlashAIMc
void FlashAIMc()
Perform flash calculation with Ni.
Definition: Bulk.cpp:3632
Bulk::CheckVe
bool CheckVe(const OCP_DBL &Vlim) const
Check if relative volume error is out of range, return false if so.
Definition: Bulk.cpp:2160
Bulk::AllocateAuxIMPEC
void AllocateAuxIMPEC()
Allocate memory for auxiliary variables used for IMPEC.
Definition: Bulk.cpp:2417
Bulk::CalMaxChange
void CalMaxChange()
Calculate max change of some variables.
Definition: Bulk.cpp:2000
Bulk::GetMixMode
USI GetMixMode() const
Return the mixture mode.
Definition: Bulk.cpp:2270
Bulk::OutFIMNi
void OutFIMNi()
FIMNi -> Ni in FIM Bulk.
Definition: Bulk.cpp:3614
Bulk::UpdateLastStepFIM
void UpdateLastStepFIM()
Update values of last step for FIM.
Definition: Bulk.cpp:2999
Bulk::PassFlashValueDeriv
void PassFlashValueDeriv(const OCP_USI &n)
Pass derivative values from Flash to Bulk after Flash calculation.
Definition: Bulk.cpp:1655
Bulk::FlashDeriv
void FlashDeriv()
Perform flash calculation with Ni and calculate derivatives.
Definition: Bulk.cpp:1401
Bulk::GetNRdSmaxP
OCP_DBL GetNRdSmaxP()
Return NRdSmaxP.
Definition: Bulk.cpp:3063
Bulk::CalVpore
void CalVpore()
Calculate volume of pore with pressure.
Definition: Bulk.cpp:1963
Bulk::CalRelResFIM
void CalRelResFIM(ResFIM &resFIM) const
Calculate relative resiual for FIM.
Definition: Bulk.cpp:2860
Bulk::CheckSetup
void CheckSetup() const
Check if error occurs in Setup.
Definition: Bulk.cpp:342
Bulk::ResetFIM
void ResetFIM()
Reset FIM.
Definition: Bulk.cpp:2950
Bulk::CalKrPcDerivAIMc
void CalKrPcDerivAIMc()
Calculate relative permeability and capillary pressure and their derivatives.
Definition: Bulk.cpp:3832
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
Bulk::GetSolAIMt
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.
Definition: Bulk.cpp:3360
Bulk::Flash
void Flash()
Perform flash calculation with Ni.
Definition: Bulk.cpp:1340
Bulk::CalKrPcDeriv
void CalKrPcDeriv()
Calculate relative permeability and capillary pressure and their derivatives.
Definition: Bulk.cpp:1929
Bulk::CheckInitVpore
void CheckInitVpore() const
Check initial pore volume.
Definition: Bulk.cpp:349
Bulk::CalNRdSmax
OCP_DBL CalNRdSmax(OCP_USI &index)
Calculate some auxiliary variable, for example, dSmax.
Definition: Bulk.cpp:3047
Bulk::GetSol01FIM
void GetSol01FIM(const vector< OCP_DBL > &u)
Get the solution for FIM after a Newton iteration???
Definition: Bulk.cpp:2858
Bulk::CheckNiFIMBulk
bool CheckNiFIMBulk() const
Check if negative Ni occurs, return false if so.
Definition: Bulk.cpp:3583
Bulk::FlashDerivBLKOILAIMc
void FlashDerivBLKOILAIMc()
Perform flash calculation with Ni in Black Oil Model.
Definition: Bulk.cpp:3769
Bulk::CalCFL
OCP_DBL CalCFL() const
Calculate the CFL number.
Definition: Bulk.cpp:2459
Bulk::FlashDerivBLKOIL
void FlashDerivBLKOIL()
Perform flash calculation with Ni in Black Oil Model.
Definition: Bulk.cpp:1428
Bulk::Setup
void Setup(const Grid &myGrid)
Allocate memory for bulk data of grid.
Definition: Bulk.cpp:183
Bulk::FlashDerivAIM
void FlashDerivAIM(const bool &IfAIMs)
Perform flash calculation with Ni and calculate derivatives.
Definition: Bulk.cpp:3173
Bulk::CalKrPcDerivAIM
void CalKrPcDerivAIM(const bool &IfAIMs)
Calculate relative permeability and capillary pressure and their derivatives.
Definition: Bulk.cpp:3309
Bulk::FlashDerivCOMPAIMc
void FlashDerivCOMPAIMc()
Perform flash calculation with Ni in Compositional Model.
Definition: Bulk.cpp:3778
Bulk::ResetFlash
void ResetFlash()
Reset variables in flash calculations.
Definition: Bulk.cpp:1880
Bulk::FlashDerivAIMc
void FlashDerivAIMc()
Perform flash calculation with Ni and calculate derivatives.
Definition: Bulk.cpp:3760
Bulk::InputParam
void InputParam(ParamReservoir &rs_param)
Input param from internal data structure ParamReservoir.
Definition: Bulk.cpp:24
Bulk::CalKrPcAIMc
void CalKrPcAIMc()
Calculate relative permeability and capillary pressure with saturation.
Definition: Bulk.cpp:3817
Bulk::InFIMNi
void InFIMNi()
Ni in FIM Bulk -> FIMNi.
Definition: Bulk.cpp:3602
Bulk::CalKrPc
void CalKrPc()
Calculate relative permeability and capillary pressure with saturation.
Definition: Bulk.cpp:1897
Bulk::CalFPR
OCP_DBL CalFPR() const
Calculate average pressure in reservoir.
Definition: Bulk.cpp:1974
Bulk::InitSjPcBo
void InitSjPcBo(const USI &tabrow)
Calculate initial equilibrium for blkoil model according to EQUIL.
Definition: Bulk.cpp:371
Bulk::UpdateLastStepIMPEC
void UpdateLastStepIMPEC()
Update value of last step for IMPEC.
Definition: Bulk.cpp:2485
Bulk::CheckNi
bool CheckNi()
Check if negative Ni occurs, return false if so.
Definition: Bulk.cpp:2114
Bulk::InitFlashDer
void InitFlashDer()
Perform flash calculation with saturations and calculate derivatives.
Definition: Bulk.cpp:1302
Bulk::AllocateAuxFIM
void AllocateAuxFIM()
Allocate memory for auxiliary variables used for FIM.
Definition: Bulk.cpp:2522
Bulk::GetSolFIM_n
void GetSolFIM_n(const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
Definition: Bulk.cpp:2704
Bulk::FlashBLKOILAIMc
void FlashBLKOILAIMc()
Perform flash calculation with Ni in Black Oil Model.
Definition: Bulk.cpp:3641
Bulk::CheckVpore
void CheckVpore() const
Check pore volume.
Definition: Bulk.cpp:360
Bulk::CheckP
bool CheckP() const
Check if negative P occurs, return false if so.
Definition: Bulk.cpp:2095
Bulk::CalRelResAIMt
void CalRelResAIMt(ResFIM &resFIM) const
Calculate relative resiual for local FIM.
Definition: Bulk.cpp:3331
Bulk::FlashDerivCOMP
void FlashDerivCOMP()
Perform flash calculation with Ni in Compositional Model.
Definition: Bulk.cpp:1440
Bulk::GetSolFIM
void GetSolFIM(const vector< OCP_DBL > &u, const OCP_DBL &dPmaxlim, const OCP_DBL &dSmaxlim)
Get the solution for FIM after a Newton iteration.
Definition: Bulk.cpp:2587
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
Grid
Basic information of computational grid, including the rock properties.
Definition: Grid.hpp:76
Grid::GetIJKGrid
void GetIJKGrid(USI &i, USI &j, USI &k, const OCP_USI &n) const
Return the 3D coordinate for object grid with Grid index.
Definition: Grid.cpp:346
OCPTable::Eval_All
USI Eval_All(const USI &j, const OCP_DBL &val, vector< OCP_DBL > &outdata, vector< OCP_DBL > &slope)
Definition: OCPTable.cpp:53
OCPTable::Display
void Display() const
Display the data of table on screen.
Definition: OCPTable.cpp:219
OCPTable::Setup
void Setup(const vector< vector< OCP_DBL >> &src)
Setup tables from existing data of table.
Definition: OCPTable.cpp:32
OCPTable::GetCol
vector< OCP_DBL > & GetCol(const USI &j)
return the jth column in table to modify or use.
Definition: OCPTable.hpp:54
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
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
Generated on Sat Oct 15 2022 00:08:52 for OpenCAEPoro by doxygen 1.9.1