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