OpenCAEPoro  0.2.0 Sep/22/2022
A simulator for multicomponent porous media flow
FlowUnit.cpp
Go to the documentation of this file.
1 
12 #include "FlowUnit.hpp"
13 #include "UtilError.hpp"
14 
15 
17 // FlowUnit_W
19 
20 void FlowUnit_W::CalKrPc(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
21  const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
22 {
23  kr_out[0] = 1;
24  pc_out[0] = 0;
25 }
26 
27 void FlowUnit_W::CalKrPcDeriv(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
28  OCP_DBL* dkrdS, OCP_DBL* dPcjdS, const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
29 {
30  kr_out[0] = 1;
31  pc_out[0] = 0;
32  dkrdS[0] = 0;
33  dPcjdS[0] = 0;
34 }
35 
36 
38 // FlowUnit_OW
40 
41 FlowUnit_OW::FlowUnit_OW(const ParamReservoir& rs_param, const USI& i)
42 {
43  SWOF.Setup(rs_param.SWOF_T.data[i]);
44  Swco = SWOF.GetCol(0)[0];
45 
46  data.resize(4, 0);
47  cdata.resize(4, 0);
48 }
49 
50 
51 void FlowUnit_OW::CalKrPc(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
52  const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
53 {
54  OCP_DBL Sw = S_in[1];
55 
56  // three phase black oil model using stone 2
57  SWOF.Eval_All(0, Sw, data, cdata);
58  OCP_DBL krw = data[1];
59  OCP_DBL kro = data[2];
60  OCP_DBL Pcwo = -data[3];
61 
62  kr_out[0] = kro;
63  kr_out[1] = krw;
64  pc_out[0] = 0;
65  pc_out[1] = Pcwo;
66 }
67 
68 void FlowUnit_OW::CalKrPcDeriv(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
69  OCP_DBL* dkrdS, OCP_DBL* dPcjdS, const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
70 {
71  OCP_DBL Sw = S_in[1];
72  SWOF.Eval_All(0, Sw, data, cdata);
73  OCP_DBL krw = data[1];
74  OCP_DBL dKrwdSw = cdata[1];
75  OCP_DBL krow = data[2];
76  OCP_DBL dKrowdSw = cdata[2];
77  OCP_DBL Pcwo = -data[3];
78  OCP_DBL dPcwdSw = -cdata[3];
79 
80  kr_out[0] = krow;
81  kr_out[1] = krw;
82  pc_out[0] = 0;
83  pc_out[1] = Pcwo;
84 
85  dkrdS[0] = 0;
86  dkrdS[1] = dKrowdSw;
87  dkrdS[2] = 0;
88  dkrdS[3] = dKrwdSw;
89 
90  dPcjdS[0] = 0;
91  dPcjdS[1] = 0;
92  dPcjdS[2] = 0;
93  dPcjdS[3] = dPcwdSw;
94 }
95 
97 // FlowUnit_OG
99 
100 FlowUnit_OG::FlowUnit_OG(const ParamReservoir& rs_param, const USI& i)
101 {
102  SGOF.Setup(rs_param.SGOF_T.data[i]);
103 
104  kroMax = SGOF.GetCol(2)[0];
105 
106  data.resize(4, 0);
107  cdata.resize(4, 0);
108 }
109 
110 
111 void FlowUnit_OG::CalKrPc(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
112  const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
113 {
114  OCP_DBL Sg = S_in[1];
115 
116  // three phase black oil model using stone 2
117  SGOF.Eval_All(0, Sg, data, cdata);
118  OCP_DBL krg = data[1];
119  OCP_DBL kro = data[2];
120  OCP_DBL Pcgo = data[3];
121 
122  kr_out[0] = kro;
123  kr_out[1] = krg;
124  pc_out[0] = 0;
125  pc_out[1] = Pcgo;
126 }
127 
128 void FlowUnit_OG::CalKrPcDeriv(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
129  OCP_DBL* dkrdS, OCP_DBL* dPcjdS, const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
130 {
131  OCP_ABORT("Not Completed Now!");
132 }
133 
135 // FlowUnit_ODGW
137 
138 OCP_DBL FlowUnit_ODGW::CalKro_Stone2(const OCP_DBL& krow, const OCP_DBL& krog,
139  const OCP_DBL& krw, const OCP_DBL& krg) const
140 {
141  // krog : oil relative permeability for a system with oil, gas and connate water
142  // krow : oil relative permeability for a system with oil and water only
143 
144  OCP_DBL kro =
145  kroMax * ((krow / kroMax + krw) * (krog / kroMax + krg) - (krw + krg));
146  if (kro < 0) kro = 0;
147 
148  return kro;
149 }
150 
151 
152 OCP_DBL FlowUnit_ODGW::CalKro_Default(const OCP_DBL& Sg, const OCP_DBL& Sw,
153  const OCP_DBL& krog, const OCP_DBL& krow) const
154 {
155  OCP_DBL tmp = Sg + Sw - Swco;
156  if (tmp <= TINY) {
157  return kroMax;
158  }
159  OCP_DBL kro = (Sg * krog + (Sw - Swco) * krow) / tmp;
160  return kro;
161 }
162 
163 
165 // FlowUnit_ODGW01
167 
168 FlowUnit_ODGW01::FlowUnit_ODGW01(const ParamReservoir& rs_param, const USI& i)
169 {
170  SWOF.Setup(rs_param.SWOF_T.data[i]);
171  SGOF.Setup(rs_param.SGOF_T.data[i]);
172 
173  kroMax = SWOF.GetCol(2)[0];
174  Swco = SWOF.GetCol(0)[0];
175 
176  data.resize(4, 0);
177  cdata.resize(4, 0);
178 
179  Generate_SWPCWG();
180 
181  Scm.resize(3, 0);
182  // oil, set to zero now
183  OCP_INT tmprow;
184  tmprow = SGOF.GetRowZero(1);
185  if (tmprow >= 0) Scm[1] = SGOF.GetCol(0)[tmprow];
186  tmprow = SWOF.GetRowZero(1);
187  if (tmprow >= 0) Scm[2] = SWOF.GetCol(0)[tmprow];
188 
189 }
190 
191 void FlowUnit_ODGW01::CalKrPc(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
192  const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
193 {
194  const OCP_DBL Sg = S_in[1];
195  const OCP_DBL Sw = S_in[2];
196 
197  // three phase black oil model using stone 2
198  SWOF.Eval_All(0, Sw, data, cdata);
199  const OCP_DBL krw = data[1];
200  const OCP_DBL krow = data[2];
201  const OCP_DBL Pcwo = -data[3];
202 
203  SGOF.Eval_All(0, Sg, data, cdata);
204  const OCP_DBL krg = data[1];
205  const OCP_DBL krog = data[2];
206  const OCP_DBL Pcgo = data[3];
207 
208  const OCP_DBL kro = CalKro_Stone2(krow, krog, krw, krg);
209  // OCP_DBL kro = CalKro_Default(Sg, Sw, krog, krow);
210 
211  kr_out[0] = kro;
212  kr_out[1] = krg;
213  kr_out[2] = krw;
214  pc_out[0] = 0;
215  pc_out[1] = Pcgo;
216  pc_out[2] = Pcwo;
217 }
218 
219 void FlowUnit_ODGW01::CalKrPcDeriv(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
220  OCP_DBL* dkrdS, OCP_DBL* dPcjdS, const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
221 {
222  OCP_DBL Sg = S_in[1];
223  OCP_DBL Sw = S_in[2];
224 
225  // three phase black oil model using stone 2
226  SWOF.Eval_All(0, Sw, data, cdata);
227  OCP_DBL krw = data[1];
228  OCP_DBL dKrwdSw = cdata[1];
229  OCP_DBL krow = data[2];
230  OCP_DBL dKrowdSw = cdata[2];
231  OCP_DBL Pcwo = -data[3];
232  OCP_DBL dPcwdSw = -cdata[3];
233 
234  SGOF.Eval_All(0, Sg, data, cdata);
235  OCP_DBL krg = data[1];
236  OCP_DBL dKrgdSg = cdata[1];
237  OCP_DBL krog = data[2];
238  OCP_DBL dKrogdSg = cdata[2];
239  OCP_DBL Pcgo = data[3];
240  OCP_DBL dPcgdSg = cdata[3];
241 
242  OCP_DBL dKrodSg{ 0 }, dKrodSw{ 0 }, kro{ 0 };
243 
244  kro = CalKro_Stone2Der(krow, krog, krw, krg, dKrwdSw, dKrowdSw, dKrgdSg, dKrogdSg,
245  dKrodSw, dKrodSg);
246  // if (kro < 0) {
247  // cout << S_in[0] << " " << S_in[1] << " " << S_in[2] << endl;
248  // kro = 0;
249  //}
250  // kro = CalKro_DefaultDer(Sg, Sw, krog, krow, dKrogdSg, dKrowdSw, dKrodSg,
251  // dKrodSw);
252 
253  kr_out[0] = kro;
254  kr_out[1] = krg;
255  kr_out[2] = krw;
256  pc_out[0] = 0;
257  pc_out[1] = Pcgo;
258  pc_out[2] = Pcwo;
259 
260  dkrdS[0] = 0;
261  dkrdS[1] = dKrodSg;
262  dkrdS[2] = dKrodSw;
263  dkrdS[3] = 0;
264  dkrdS[4] = dKrgdSg;
265  dkrdS[5] = 0;
266  dkrdS[6] = 0;
267  dkrdS[7] = 0;
268  dkrdS[8] = dKrwdSw;
269 
270  dPcjdS[0] = 0;
271  dPcjdS[1] = 0;
272  dPcjdS[2] = 0;
273  dPcjdS[3] = 0;
274  dPcjdS[4] = dPcgdSg;
275  dPcjdS[5] = 0;
276  dPcjdS[6] = 0;
277  dPcjdS[7] = 0;
278  dPcjdS[8] = dPcwdSw;
279 }
280 
281 OCP_DBL FlowUnit_ODGW01::CalKro_Stone2Der(OCP_DBL krow, OCP_DBL krog, OCP_DBL krw, OCP_DBL krg,
282  OCP_DBL dkrwdSw, OCP_DBL dkrowdSw, OCP_DBL dkrgdSg,
283  OCP_DBL dkrogdSg, OCP_DBL& out_dkrodSw,
284  OCP_DBL& out_dkrodSg) const
285 {
286  OCP_DBL kro, dkrodSw, dkrodSg;
287  kro = kroMax * ((krow / kroMax + krw) * (krog / kroMax + krg) - (krw + krg));
288 
289  dkrodSw =
290  kroMax * ((dkrowdSw / kroMax + dkrwdSw) * (krog / kroMax + krg) - (dkrwdSw));
291  dkrodSg =
292  kroMax * ((krow / kroMax + krw) * (dkrogdSg / kroMax + dkrgdSg) - (dkrgdSg));
293 
294  if (kro < 0) {
295  kro = 0;
296  dkrodSg = 0;
297  dkrodSw = 0;
298  }
299  out_dkrodSw = dkrodSw;
300  out_dkrodSg = dkrodSg;
301  return kro;
302 }
303 
304 OCP_DBL FlowUnit_ODGW01::CalKro_DefaultDer(const OCP_DBL& Sg, const OCP_DBL& Sw,
305  const OCP_DBL& krog, const OCP_DBL& krow,
306  const OCP_DBL& dkrogSg, const OCP_DBL& dkrowSw,
307  OCP_DBL& dkroSg, OCP_DBL& dkroSw) const
308 {
309  OCP_DBL tmp = Sg + Sw - Swco;
310  if (tmp <= TINY) {
311  dkroSg = 0;
312  dkroSw = 0;
313  return kroMax;
314  }
315  OCP_DBL kro = (Sg * krog + (Sw - Swco) * krow) / tmp;
316  dkroSg = (krog + Sg * dkrogSg - kro) / tmp;
317  dkroSw = (krow + (Sw - Swco) * dkrowSw - kro) / tmp;
318  return kro;
319 }
320 
321 void FlowUnit_ODGW01::Generate_SWPCWG()
322 {
323  if (SGOF.IsEmpty()) OCP_ABORT("SGOF is missing!");
324  if (SWOF.IsEmpty()) OCP_ABORT("SWOF is missing!");
325 
326  std::vector<OCP_DBL> Sw(SWOF.GetCol(0));
327  std::vector<OCP_DBL> Pcow(SWOF.GetCol(3));
328  USI n = Sw.size();
329  for (USI i = 0; i < n; i++) {
330  OCP_DBL Pcgo = SGOF.Eval(0, 1 - Sw[i], 3);
331  Pcow[i] += Pcgo; // Pcgw
332  }
333 
334  SWPCGW.PushCol(Sw);
335  SWPCGW.PushCol(Pcow);
336  SWPCGW.SetRowCol();
337 }
338 
340 // FlowUnit_ODGW01_Miscible
342 
343 void FlowUnit_ODGW01_Miscible::CalKrPc(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
344  const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
345 {
346  surTen = MySurTen;
347  if (surTen >= surTenRef || surTen < TINY) {
348  MyFk = 1; MyFp = 1;
349  FlowUnit_ODGW01::CalKrPc(S_in, kr_out, pc_out, surTen, surTen, surTen);
350  }
351  else {
352  OCP_DBL So = S_in[0];
353  OCP_DBL Sg = S_in[1];
354  OCP_DBL Sw = S_in[2];
355 
356  Fk = min(1.0, pow(surTen / surTenRef, Fkexp));
357  Fp = min(surTenPc, surTen / surTenRef);
358 
359  SWOF.Eval_All(0, Sw, data, cdata);
360  const OCP_DBL krw = data[1];
361  OCP_DBL krow = data[2];
362  const OCP_DBL Pcwo = -data[3];
363 
364  SGOF.Eval_All(0, Sg, data, cdata);
365  OCP_DBL krg = data[1];
366  OCP_DBL krog = data[2];
367  const OCP_DBL Pcgo = data[3] * Fp;
368 
369  OCP_DBL kro = CalKro_Stone2(krow, krog, krw, krg);
370  // OCP_DBL kro = CalKro_Default(Sg, Sw, krog, krow);
371  OCP_DBL krgt = SGOF.Eval(0, (1 - Sw), 1);
372  OCP_DBL krh = 0.5 * (krow + krgt);
373 
374  // from CMG, see *SIGMA
375  kro = Fk * kro + (1 - Fk) * krh * So / (1 - Sw);
376  krg = Fk * krg + (1 - Fk) * krh * Sg / (1 - Sw);
377 
378  kr_out[0] = kro;
379  kr_out[1] = krg;
380  kr_out[2] = krw;
381  pc_out[0] = 0;
382  pc_out[1] = Pcgo;
383  pc_out[2] = Pcwo;
384 
385  MyFk = Fk;
386  MyFp = Fp;
387 
388  // cout << Fk << endl;
389  }
390 }
391 
392 
393 void FlowUnit_ODGW01_Miscible::CalKrPcDeriv(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
394  OCP_DBL* dkrdS, OCP_DBL* dPcjdS, const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
395 {
396  surTen = MySurTen;
397  if (surTen >= surTenRef || surTen < TINY) {
398  MyFk = 1; MyFp = 1;
399  FlowUnit_ODGW01::CalKrPcDeriv(S_in, kr_out, pc_out, dkrdS, dPcjdS, surTen, surTen, surTen);
400  }
401  else {
402  const OCP_DBL So = S_in[0];
403  const OCP_DBL Sg = S_in[1];
404  const OCP_DBL Sw = S_in[2];
405 
406  Fk = min(1.0, pow(surTen / surTenRef, Fkexp));
407  Fp = min(surTenPc, surTen / surTenRef);
408 
409  SWOF.Eval_All(0, Sw, data, cdata);
410  const OCP_DBL krw = data[1];
411  const OCP_DBL dKrwdSw = cdata[1];
412  const OCP_DBL krow = data[2];
413  const OCP_DBL dKrowdSw = cdata[2];
414  const OCP_DBL Pcwo = -data[3];
415  const OCP_DBL dPcwdSw = -cdata[3];
416 
417  SGOF.Eval_All(0, Sg, data, cdata);
418  OCP_DBL krg = data[1];
419  const OCP_DBL dKrgdSg = cdata[1];
420  const OCP_DBL krog = data[2];
421  const OCP_DBL dKrogdSg = cdata[2];
422  const OCP_DBL Pcgo = data[3] * Fp;
423  const OCP_DBL dPcgdSg = cdata[3] * Fp;
424 
425  OCP_DBL dKrodSg{ 0 }, dKrodSw{ 0 }, kro{ 0 };
426  kro = CalKro_Stone2Der(krow, krog, krw, krg, dKrwdSw, dKrowdSw, dKrgdSg, dKrogdSg,
427  dKrodSw, dKrodSg);
428 
429  OCP_DBL dkrgd1_Sw = 0;
430  const OCP_DBL krgt = SGOF.Eval(0, (1 - Sw), 1, dkrgd1_Sw);
431  const OCP_DBL krh = 0.5 * (krow + krgt);
432 
433  // from CMG, see *SIGMA
434  kro = Fk * kro + (1 - Fk) * krh * So / (1 - Sw);
435  krg = Fk * krg + (1 - Fk) * krh * Sg / (1 - Sw);
436 
437  kr_out[0] = kro;
438  kr_out[1] = krg;
439  kr_out[2] = krw;
440  pc_out[0] = 0;
441  pc_out[1] = Pcgo;
442  pc_out[2] = Pcwo;
443 
444  OCP_DBL dkrhdSw = 0.5 * (dKrowdSw - dkrgd1_Sw);
445  OCP_DBL temp = (1 - Fk) / (1 - Sw) * (krh / (1 - Sw) + dkrhdSw);
446 
447  dkrdS[0] = (1 - Fk) * krh / (1 - Sw);
448  dkrdS[1] = Fk * dKrodSg;
449  dkrdS[2] = Fk * dKrodSw + temp * So;
450  dkrdS[3] = 0;
451  dkrdS[4] = Fk * dKrgdSg + (1 - Fk) * krh / (1 - Sw);
452  dkrdS[5] = temp * Sg;
453  dkrdS[6] = 0;
454  dkrdS[7] = 0;
455  dkrdS[8] = dKrwdSw;
456 
457  dPcjdS[0] = 0;
458  dPcjdS[1] = 0;
459  dPcjdS[2] = 0;
460  dPcjdS[3] = 0;
461  dPcjdS[4] = dPcgdSg;
462  dPcjdS[5] = 0;
463  dPcjdS[6] = 0;
464  dPcjdS[7] = 0;
465  dPcjdS[8] = dPcwdSw;
466 
467  MyFk = Fk;
468  MyFp = Fp;
469 
470  // cout << Fk << endl;
471  }
472 }
473 
475 // FlowUnit_ODGW02
477 
478 FlowUnit_ODGW02::FlowUnit_ODGW02(const ParamReservoir& rs_param, const USI& i)
479 {
480  SWFN.Setup(rs_param.SWFN_T.data[i]);
481  SGFN.Setup(rs_param.SGFN_T.data[i]);
482  SOF3.Setup(rs_param.SOF3_T.data[i]);
483 
484  kroMax = SOF3.GetCol(1).back();
485  Swco = SWFN.GetCol(0)[0];
486 
487  data.resize(3, 0);
488  cdata.resize(3, 0);
489 
490  Generate_SWPCWG();
491 }
492 
493 
494 void FlowUnit_ODGW02::CalKrPc(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
495  const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
496 {
497  OCP_DBL So = S_in[0];
498  OCP_DBL Sg = S_in[1];
499  OCP_DBL Sw = S_in[2];
500 
501  SWFN.Eval_All(0, Sw, data, cdata);
502  OCP_DBL krw = data[1];
503  OCP_DBL Pcwo = -data[2];
504 
505  SGFN.Eval_All(0, Sg, data, cdata);
506  OCP_DBL krg = data[1];
507  OCP_DBL Pcgo = data[2];
508 
509  SOF3.Eval_All(0, So, data, cdata);
510  OCP_DBL krow = data[1];
511  OCP_DBL krog = data[2];
512 
513  // OCP_DBL kro = CalKro_Stone2(krow, krog, krw, krg);
514  OCP_DBL kro = CalKro_Default(Sg, Sw, krog, krow);
515 
516 
517  kr_out[0] = kro;
518  kr_out[1] = krg;
519  kr_out[2] = krw;
520  pc_out[0] = 0;
521  pc_out[1] = Pcgo;
522  pc_out[2] = Pcwo;
523 }
524 
525 void FlowUnit_ODGW02::CalKrPcDeriv(const OCP_DBL* S_in, OCP_DBL* kr_out, OCP_DBL* pc_out,
526  OCP_DBL* dkrdS, OCP_DBL* dPcjdS, const OCP_DBL& MySurTen, OCP_DBL& MyFk, OCP_DBL& MyFp)
527 {
528  OCP_DBL So = S_in[0];
529  OCP_DBL Sg = S_in[1];
530  OCP_DBL Sw = S_in[2];
531 
532  SWFN.Eval_All(0, Sw, data, cdata);
533  OCP_DBL krw = data[1];
534  OCP_DBL dKrwdSw = cdata[1];
535  OCP_DBL Pcwo = -data[2];
536  OCP_DBL dPcwodSw = -cdata[2];
537 
538  SGFN.Eval_All(0, Sg, data, cdata);
539  OCP_DBL krg = data[1];
540  OCP_DBL dKrgdSg = cdata[1];
541  OCP_DBL Pcgo = data[2];
542  OCP_DBL dPcgodSg = cdata[2];
543 
544  SOF3.Eval_All(0, So, data, cdata);
545  OCP_DBL krow = data[1];
546  OCP_DBL dKrowdSo = cdata[1];
547  OCP_DBL krog = data[2];
548  OCP_DBL dKrogdSo = cdata[2];
549 
550  OCP_DBL dKroSo = 0;
551  OCP_DBL kro = CalKro_Stone2Der(krow, krog, krw, krg, dKrwdSw, dKrowdSo, dKrgdSg, dKrogdSo, dKroSo);
552  // OCP_DBL kro = CalKro_DefaultDer(Sg, Sw, krog, krow, dKrogdSo, dKrowdSo, dKroSo);
553 
554  kr_out[0] = kro;
555  kr_out[1] = krg;
556  kr_out[2] = krw;
557  pc_out[0] = 0;
558  pc_out[1] = Pcgo;
559  pc_out[2] = Pcwo;
560 
561  dkrdS[0] = dKroSo;
562  dkrdS[1] = 0;
563  dkrdS[2] = 0;
564  dkrdS[3] = 0;
565  dkrdS[4] = dKrgdSg;
566  dkrdS[5] = 0;
567  dkrdS[6] = 0;
568  dkrdS[7] = 0;
569  dkrdS[8] = dKrwdSw;
570 
571  dPcjdS[0] = 0;
572  dPcjdS[1] = 0;
573  dPcjdS[2] = 0;
574  dPcjdS[3] = 0;
575  dPcjdS[4] = dPcgodSg;
576  dPcjdS[5] = 0;
577  dPcjdS[6] = 0;
578  dPcjdS[7] = 0;
579  dPcjdS[8] = dPcwodSw;
580 }
581 
582 OCP_DBL FlowUnit_ODGW02::CalKro_Stone2Der(OCP_DBL krow, OCP_DBL krog, OCP_DBL krw, OCP_DBL krg,
583  OCP_DBL dkrwdSw, OCP_DBL dkrowdSo, OCP_DBL dkrgdSg,
584  OCP_DBL dkrogdSo, OCP_DBL& out_dkrodSo) const
585 {
586  OCP_DBL kro, dKrodSo;
587 
588  kro = kroMax * ((krow / kroMax + krw) * (krog / kroMax + krg) - (krw + krg));
589  dKrodSo = dkrowdSo * (krog / kroMax + krg) + dkrogdSo * (krow / kroMax + krw);
590 
591  if (kro < 0) {
592  kro = 0;
593  dKrodSo = 0;
594  }
595  out_dkrodSo = dKrodSo;
596  return kro;
597 }
598 
599 OCP_DBL FlowUnit_ODGW02::CalKro_DefaultDer(const OCP_DBL& Sg, const OCP_DBL& Sw,
600  const OCP_DBL& krog, const OCP_DBL& krow,
601  const OCP_DBL& dkrogdSo, const OCP_DBL& dkrowdSo,
602  OCP_DBL& out_dkrodSo) const
603 {
604  OCP_DBL tmp = Sg + Sw - Swco;
605  OCP_DBL kro = (Sg * krog + (Sw - Swco) * krow) / tmp;
606  out_dkrodSo = (Sg * dkrogdSo + (Sw - Swco) * dkrowdSo) / tmp;
607 
608  if (tmp <= TINY) {
609  kro = kroMax;
610  out_dkrodSo = 0;
611  }
612  return kro;
613 }
614 
615 void FlowUnit_ODGW02::Generate_SWPCWG()
616 {
617  if (SWFN.IsEmpty()) OCP_ABORT("SWFN is missing!");
618  if (SGFN.IsEmpty()) OCP_ABORT("SGFN is missing!");
619 
620  std::vector<OCP_DBL> Sw(SWFN.GetCol(0));
621  std::vector<OCP_DBL> Pcow(SWFN.GetCol(2));
622  USI n = Sw.size();
623  for (USI i = 0; i < n; i++) {
624  OCP_DBL Pcgo = SGFN.Eval(0, 1 - Sw[i], 2);
625  Pcow[i] += Pcgo; // pcgw
626  }
627 
628  SWPCGW.PushCol(Sw);
629  SWPCGW.PushCol(Pcow);
630  SWPCGW.SetRowCol();
631 }
632 
633 
634 
635 /*----------------------------------------------------------------------------*/
636 /* Brief Change History of This File */
637 /*----------------------------------------------------------------------------*/
638 /* Author Date Actions */
639 /*----------------------------------------------------------------------------*/
640 /* Shizhe Li Oct/01/2021 Create file */
641 /* Chensong Zhang Oct/15/2021 Format file */
642 /*----------------------------------------------------------------------------*/
FlowUnit.hpp
FlowUnit class declaration.
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:36
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
OCP_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:25
UtilError.hpp
Logging error and warning messages.
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47
FlowUnit_ODGW01_Miscible::CalKrPcDeriv
void CalKrPcDeriv(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, OCP_DBL *dkrdS, OCP_DBL *dPcjdS, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate derivatives of relative permeability and capillary pressure.
Definition: FlowUnit.cpp:393
FlowUnit_ODGW01_Miscible::CalKrPc
void CalKrPc(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate relative permeability and capillary pressure.
Definition: FlowUnit.cpp:343
FlowUnit_ODGW01::CalKrPcDeriv
virtual void CalKrPcDeriv(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, OCP_DBL *dkrdS, OCP_DBL *dPcjdS, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate derivatives of relative permeability and capillary pressure.
Definition: FlowUnit.cpp:219
FlowUnit_ODGW01::SGOF
OCPTable SGOF
saturation table about gas and oil.
Definition: FlowUnit.hpp:219
FlowUnit_ODGW01::cdata
vector< OCP_DBL > cdata
container to store the slopes of interpolation.
Definition: FlowUnit.hpp:222
FlowUnit_ODGW01::CalKrPc
virtual void CalKrPc(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate relative permeability and capillary pressure.
Definition: FlowUnit.cpp:191
FlowUnit_ODGW01::SWPCGW
OCPTable SWPCGW
Definition: FlowUnit.hpp:223
FlowUnit_ODGW01::SWOF
OCPTable SWOF
saturation table about water and oil.
Definition: FlowUnit.hpp:220
FlowUnit_ODGW01::data
vector< OCP_DBL > data
container to store the values of interpolation.
Definition: FlowUnit.hpp:221
FlowUnit_ODGW02::CalKrPcDeriv
void CalKrPcDeriv(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, OCP_DBL *dkrdS, OCP_DBL *dPcjdS, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate derivatives of relative permeability and capillary pressure.
Definition: FlowUnit.cpp:525
FlowUnit_ODGW02::CalKrPc
void CalKrPc(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate relative permeability and capillary pressure.
Definition: FlowUnit.cpp:494
FlowUnit_ODGW::kroMax
OCP_DBL kroMax
oil relative permeability in the presence of connate water only, used in stone2
Definition: FlowUnit.hpp:172
FlowUnit_ODGW::Scm
vector< OCP_DBL > Scm
critical saturation when phase becomes mobile / immobile
Definition: FlowUnit.hpp:174
FlowUnit_ODGW::Swco
OCP_DBL Swco
Saturation of connate water.
Definition: FlowUnit.hpp:173
FlowUnit_OG::CalKrPc
void CalKrPc(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate relative permeability and capillary pressure.
Definition: FlowUnit.cpp:111
FlowUnit_OG::CalKrPcDeriv
void CalKrPcDeriv(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, OCP_DBL *dkrdS, OCP_DBL *dPcjdS, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate derivatives of relative permeability and capillary pressure.
Definition: FlowUnit.cpp:128
FlowUnit_OW::CalKrPcDeriv
void CalKrPcDeriv(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, OCP_DBL *dkrdS, OCP_DBL *dPcjdS, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate derivatives of relative permeability and capillary pressure.
Definition: FlowUnit.cpp:68
FlowUnit_OW::CalKrPc
void CalKrPc(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate relative permeability and capillary pressure.
Definition: FlowUnit.cpp:51
FlowUnit_W::CalKrPc
void CalKrPc(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate relative permeability and capillary pressure.
Definition: FlowUnit.cpp:20
FlowUnit_W::CalKrPcDeriv
void CalKrPcDeriv(const OCP_DBL *S_in, OCP_DBL *kr_out, OCP_DBL *pc_out, OCP_DBL *dkrdS, OCP_DBL *dPcjdS, const OCP_DBL &MySurTen, OCP_DBL &MyFk, OCP_DBL &MyFp) override
Calculate derivatives of relative permeability and capillary pressure.
Definition: FlowUnit.cpp:27
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::SetRowCol
void SetRowCol()
Setup row nums and col nums of tables, initialize the bId.
Definition: OCPTable.hpp:57
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::GetRowZero
OCP_INT GetRowZero(const USI &mycol) const
return the row index of the last zero of some colnum, which is sorted in increasing order.
Definition: OCPTable.cpp:41
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
OCPTable::PushCol
void PushCol(const vector< OCP_DBL > &v)
push v into the last column of table.
Definition: OCPTable.hpp:51
ParamReservoir::SOF3_T
TableSet SOF3_T
Table set of SOF3.
Definition: ParamReservoir.hpp:198
ParamReservoir::SGOF_T
TableSet SGOF_T
Table set of SGOF.
Definition: ParamReservoir.hpp:197
ParamReservoir::SGFN_T
TableSet SGFN_T
Table set of SGFN.
Definition: ParamReservoir.hpp:196
ParamReservoir::SWOF_T
TableSet SWOF_T
Table set of SWOF.
Definition: ParamReservoir.hpp:195
ParamReservoir::SWFN_T
TableSet SWFN_T
Table set of SWFN.
Definition: ParamReservoir.hpp:194
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