OpenCAEPoro  v0.5.0
A simulator for multicomponent porous media flow
Well.cpp
Go to the documentation of this file.
1 
12 #include "Well.hpp"
13 #include <cmath>
14 
15 void Well::InputPerfo(const WellParam& well)
16 {
17  OCP_FUNCNAME;
18 
19  numPerf = well.I_perf.size();
20  perf.resize(numPerf);
21  for (USI p = 0; p < numPerf; p++) {
22  perf[p].I = well.I_perf[p] - 1;
23  perf[p].J = well.J_perf[p] - 1;
24  perf[p].K = well.K_perf[p] - 1;
25  perf[p].WI = well.WI[p];
26  perf[p].radius = well.diameter[p] / 2.0;
27  perf[p].kh = well.kh[p];
28  perf[p].skinFactor = well.skinFactor[p];
29  if (well.direction[p] == "X" || well.direction[p] == "x") {
30  perf[p].direction = X_DIRECTION;
31  } else if (well.direction[p] == "Y" || well.direction[p] == "y") {
32  perf[p].direction = Y_DIRECTION;
33  } else if (well.direction[p] == "Z" || well.direction[p] == "z") {
34  perf[p].direction = Z_DIRECTION;
35  } else {
36  OCP_ABORT("Wrong direction of perforations!");
37  }
38  }
39 }
40 
41 void Well::Setup(const Grid& gd, const Bulk& bk, const vector<SolventINJ>& sols)
42 {
43  OCP_FUNCNAME;
44 
45  numCom = bk.numCom;
46  numPhase = bk.numPhase;
47  flashCal = bk.flashCal;
48 
49  qi_lbmol.resize(numCom);
50  prodWeight.resize(numCom);
51  prodRate.resize(numPhase);
52 
53  for (auto& opt : optSet) {
54  if (!opt.state) continue;
55  if (!bk.ifThermal) {
56  opt.injTemp = bk.rsTemp;
57  }
58  flashCal[0]->SetupWellOpt(opt, sols, Psurf, Tsurf);
59  }
60 
61  // Perf
62  USI pp = 0;
63  for (USI p = 0; p < numPerf; p++) {
64  OCP_USI pId = perf[p].K * gd.nx * gd.ny + perf[p].J * gd.nx + perf[p].I;
65  if (gd.map_All2Flu[pId].IsAct()) {
66 
67  perf[pp] = perf[p];
68  perf[pp].state = OPEN;
69  perf[pp].location = gd.map_All2Act[pId].GetId();
70  perf[pp].depth = bk.depth[perf[pp].location];
71  perf[pp].multiplier = 1;
72  perf[pp].qi_lbmol.resize(numCom);
73  perf[pp].transj.resize(numPhase);
74  perf[pp].qj_ft3.resize(numPhase);
75  pp++;
76  } else {
77  OCP_WARNING("Perforation is in non-fluid Bulk!");
78  }
79  }
80  numPerf = pp;
81  perf.resize(numPerf);
82  // dG
83  dG.resize(numPerf, 0);
84  ldG = dG;
85 
86  if (depth < 0) depth = perf[0].depth;
87 
88  CalWI_Peaceman(bk);
89  // test
90  // ShowPerfStatus(bk);
91 }
92 
93 void Well::CalWI_Peaceman(const Bulk& myBulk)
94 {
95  OCP_FUNCNAME;
96 
97  // this fomular needs to be carefully checked !
98  // especially the dz
99 
100  for (USI p = 0; p < numPerf; p++) {
101  if (perf[p].WI > 0) {
102  break;
103  } else {
104  const OCP_USI Idb = perf[p].location;
105  const OCP_DBL dx = myBulk.dx[Idb];
106  const OCP_DBL dy = myBulk.dy[Idb];
107  const OCP_DBL dz = myBulk.dz[Idb] * myBulk.ntg[Idb];
108  OCP_DBL ro = 0;
109  switch (perf[p].direction) {
110  case X_DIRECTION:
111  {
112  const OCP_DBL kykz = myBulk.rockKy[Idb] * myBulk.rockKz[Idb];
113  const OCP_DBL ky_kz = myBulk.rockKy[Idb] / myBulk.rockKz[Idb];
114  assert(kykz > 0);
115  ro = 0.28 * pow((dy * dy * pow(1 / ky_kz, 0.5) +
116  dz * dz * pow(ky_kz, 0.5)),
117  0.5);
118  ro /= (pow(ky_kz, 0.25) + pow(1 / ky_kz, 0.25));
119 
120  if (perf[p].kh < 0) {
121  perf[p].kh = (dx * pow(kykz, 0.5));
122  }
123  break;
124  }
125  case Y_DIRECTION:
126  {
127  const OCP_DBL kzkx = myBulk.rockKz[Idb] * myBulk.rockKx[Idb];
128  const OCP_DBL kz_kx = myBulk.rockKz[Idb] / myBulk.rockKx[Idb];
129  assert(kzkx > 0);
130  ro = 0.28 * pow((dz * dz * pow(1 / kz_kx, 0.5) +
131  dx * dx * pow(kz_kx, 0.5)),
132  0.5);
133  ro /= (pow(kz_kx, 0.25) + pow(1 / kz_kx, 0.25));
134 
135  if (perf[p].kh < 0) {
136  perf[p].kh = (dy * pow(kzkx, 0.5));
137  }
138  break;
139  }
140  case Z_DIRECTION:
141  {
142  const OCP_DBL kxky = myBulk.rockKx[Idb] * myBulk.rockKy[Idb];
143  const OCP_DBL kx_ky = myBulk.rockKx[Idb] / myBulk.rockKy[Idb];
144  assert(kxky > 0);
145  ro = 0.28 * pow((dx * dx * pow(1 / kx_ky, 0.5) +
146  dy * dy * pow(kx_ky, 0.5)),
147  0.5);
148  ro /= (pow(kx_ky, 0.25) + pow(1 / kx_ky, 0.25));
149 
150  if (perf[p].kh < 0) {
151  perf[p].kh = (dz * pow(kxky, 0.5));
152  }
153  break;
154  }
155  default:
156  OCP_ABORT("Wrong direction of perforations!");
157  }
158  perf[p].WI = CONV2 * (2 * PI) * perf[p].kh /
159  (log(ro / perf[p].radius) + perf[p].skinFactor);
160  }
161  }
162 }
163 
164 void Well::CalTrans(const Bulk& myBulk)
165 {
166  OCP_FUNCNAME;
167 
168  if (opt.type == INJ) {
169  for (USI p = 0; p < numPerf; p++) {
170  perf[p].transINJ = 0;
171  OCP_USI k = perf[p].location;
172  OCP_DBL temp = CONV1 * perf[p].WI * perf[p].multiplier;
173 
174  // single phase
175  for (USI j = 0; j < numPhase; j++) {
176  perf[p].transj[j] = 0;
177  OCP_USI id = k * numPhase + j;
178  if (myBulk.phaseExist[id]) {
179  perf[p].transj[j] = temp * myBulk.kr[id] / myBulk.mu[id];
180  perf[p].transINJ += perf[p].transj[j];
181  }
182  }
183  if (ifUseUnweight) {
184  perf[p].transINJ = perf[p].WI;
185  }
186  }
187  } else {
188  for (USI p = 0; p < numPerf; p++) {
189  OCP_USI k = perf[p].location;
190  OCP_DBL temp = CONV1 * perf[p].WI * perf[p].multiplier;
191 
192  // multi phase
193  for (USI j = 0; j < numPhase; j++) {
194  perf[p].transj[j] = 0;
195  OCP_USI id = k * numPhase + j;
196  if (myBulk.phaseExist[id]) {
197  perf[p].transj[j] = temp * myBulk.kr[id] / myBulk.mu[id];
198  }
199  }
200  }
201  }
202 }
203 
204 void Well::CalFlux(const Bulk& myBulk, const OCP_BOOL ReCalXi)
205 {
206  OCP_FUNCNAME;
207 
208  // cout << name << endl;
209  fill(qi_lbmol.begin(), qi_lbmol.end(), 0.0);
210 
211  if (opt.type == INJ) {
212 
213  for (USI p = 0; p < numPerf; p++) {
214  perf[p].P = bhp + dG[p];
215  OCP_USI k = perf[p].location;
216  OCP_DBL dP = myBulk.P[k] - perf[p].P;
217 
218  perf[p].qt_ft3 = perf[p].transINJ * dP;
219 
220  if (ReCalXi) {
221  USI pvtnum = myBulk.PVTNUM[k];
222  perf[p].xi = myBulk.flashCal[pvtnum]->XiPhase(
223  perf[p].P, opt.injTemp, &opt.injZi[0], opt.injProdPhase);
224  }
225  for (USI i = 0; i < numCom; i++) {
226  perf[p].qi_lbmol[i] = perf[p].qt_ft3 * perf[p].xi * opt.injZi[i];
227  qi_lbmol[i] += perf[p].qi_lbmol[i];
228  }
229  }
230  } else {
231 
232  for (USI p = 0; p < numPerf; p++) {
233  perf[p].P = bhp + dG[p];
234  OCP_USI k = perf[p].location;
235  perf[p].qt_ft3 = 0;
236  fill(perf[p].qi_lbmol.begin(), perf[p].qi_lbmol.end(), 0.0);
237  fill(perf[p].qj_ft3.begin(), perf[p].qj_ft3.end(), 0.0);
238 
239  for (USI j = 0; j < numPhase; j++) {
240  OCP_USI id = k * numPhase + j;
241  if (myBulk.phaseExist[id]) {
242  OCP_DBL dP = myBulk.Pj[id] - perf[p].P;
243 
244  perf[p].qj_ft3[j] = perf[p].transj[j] * dP;
245  perf[p].qt_ft3 += perf[p].qj_ft3[j];
246 
247  OCP_DBL xi = myBulk.xi[id];
248  OCP_DBL xij;
249  for (USI i = 0; i < numCom; i++) {
250  xij = myBulk.xij[id * numCom + i];
251  perf[p].qi_lbmol[i] += perf[p].qj_ft3[j] * xi * xij;
252  }
253  }
254  }
255  for (USI i = 0; i < numCom; i++) qi_lbmol[i] += perf[p].qi_lbmol[i];
256  }
257  }
258 }
259 
263 OCP_DBL Well::CalInjRateMaxBHP(const Bulk& myBulk)
264 {
265  OCP_FUNCNAME;
266 
267  OCP_DBL qj = 0;
268  OCP_DBL Pwell = opt.maxBHP;
269 
270  for (USI p = 0; p < numPerf; p++) {
271 
272  OCP_DBL Pperf = Pwell + dG[p];
273  OCP_USI k = perf[p].location;
274 
275  USI pvtnum = myBulk.PVTNUM[k];
276  OCP_DBL xi = myBulk.flashCal[pvtnum]->XiPhase(Pperf, opt.injTemp, &opt.injZi[0],
277  opt.injProdPhase);
278 
279  OCP_DBL dP = Pperf - myBulk.P[k];
280  qj += perf[p].transINJ * xi * dP;
281  }
282  return qj;
283 }
284 
288 OCP_DBL Well::CalProdRateMinBHP(const Bulk& myBulk)
289 {
290  OCP_FUNCNAME;
291 
292  OCP_DBL qj = 0;
293  OCP_DBL Pwell = opt.minBHP;
294 
295  vector<OCP_DBL> tmpQi_lbmol(numCom, 0);
296  vector<OCP_DBL> tmpQj(numPhase, 0);
297 
298  for (USI p = 0; p < numPerf; p++) {
299 
300  OCP_DBL Pperf = Pwell + dG[p];
301  OCP_USI k = perf[p].location;
302 
303  for (USI j = 0; j < numPhase; j++) {
304  OCP_USI id = k * numPhase + j;
305  if (myBulk.phaseExist[id]) {
306  OCP_DBL dP = myBulk.Pj[id] - Pperf;
307  OCP_DBL temp = perf[p].transj[j] * myBulk.xi[id] * dP;
308  for (USI i = 0; i < numCom; i++) {
309  tmpQi_lbmol[i] += myBulk.xij[id * numCom + i] * temp;
310  }
311  }
312  }
313  }
314  flashCal[0]->CalProdRate(Psurf, Tsurf, &tmpQi_lbmol[0], tmpQj);
315  for (USI j = 0; j < numPhase; j++) {
316  qj += tmpQj[j] * opt.prodPhaseWeight[j];
317  }
318 
319  return qj;
320 }
321 
322 void Well::CalInjQj(const Bulk& myBulk, const OCP_DBL& dt)
323 {
324  OCP_FUNCNAME;
325 
326  OCP_DBL qj = 0;
327 
328  for (USI i = 0; i < numCom; i++) {
329  qj += qi_lbmol[i];
330  }
331  if (opt.fluidType == "WAT") {
332  WWIR = -qj / opt.factorINJ;
333  WWIT += WWIR * dt;
334  } else {
335  WGIR = -qj / opt.factorINJ; // Mscf or lbmol -> Mscf
336  WGIT += WGIR * dt;
337  }
338 }
339 
340 void Well::CalProdQj(const Bulk& myBulk, const OCP_DBL& dt)
341 {
342 
343  flashCal[0]->CalProdRate(Psurf, Tsurf, &qi_lbmol[0], prodRate);
344  USI iter = 0;
345  if (myBulk.oil) WOPR = prodRate[iter++];
346  if (myBulk.gas) WGPR = prodRate[iter++];
347  if (myBulk.water) WWPR = prodRate[iter++];
348 
349  WOPT += WOPR * dt;
350  WGPT += WGPR * dt;
351  WWPT += WWPR * dt;
352 }
353 
357 void Well::CaldG(const Bulk& myBulk)
358 {
359  OCP_FUNCNAME;
360 
361  if (opt.type == INJ)
362  CalInjdG(myBulk);
363  else
364  CalProddG01(myBulk);
365 }
366 
367 void Well::CalInjdG(const Bulk& myBulk)
368 {
369  OCP_FUNCNAME;
370 
371  const OCP_DBL maxlen = 10;
372  USI seg_num = 0;
373  OCP_DBL seg_len = 0;
374  vector<OCP_DBL> dGperf(numPerf, 0);
375 
376  if (depth <= perf.front().depth) {
377  // Well is higher
378  for (OCP_INT p = numPerf - 1; p >= 0; p--) {
379  if (p == 0) {
380  seg_num = ceil(fabs((perf[0].depth - depth) / maxlen));
381  if (seg_num == 0) continue;
382  seg_len = (perf[0].depth - depth) / seg_num;
383  } else {
384  seg_num = ceil(fabs((perf[p].depth - perf[p - 1].depth) / maxlen));
385  if (seg_num == 0) continue;
386  seg_len = (perf[p].depth - perf[p - 1].depth) / seg_num;
387  }
388  OCP_USI n = perf[p].location;
389  perf[p].P = bhp + dG[p];
390  OCP_DBL Pperf = perf[p].P;
391  OCP_DBL Ptmp = Pperf;
392 
393  USI pvtnum = myBulk.PVTNUM[n];
394  for (USI i = 0; i < seg_num; i++) {
395  Ptmp -= myBulk.flashCal[pvtnum]->RhoPhase(
396  Ptmp, 0, opt.injTemp, opt.injZi.data(), opt.injProdPhase) *
397  GRAVITY_FACTOR * seg_len;
398  }
399  dGperf[p] = Pperf - Ptmp;
400  }
401  dG[0] = dGperf[0];
402  for (USI p = 1; p < numPerf; p++) {
403  dG[p] = dG[p - 1] + dGperf[p];
404  }
405  } else if (depth >= perf[numPerf - 1].depth) {
406  // Well is lower
407  for (USI p = 0; p < numPerf; p++) {
408  if (p == numPerf - 1) {
409  seg_num = ceil(fabs((depth - perf[numPerf - 1].depth) / maxlen));
410  if (seg_num == 0) continue;
411  seg_len = (depth - perf[numPerf - 1].depth) / seg_num;
412  } else {
413  seg_num = ceil(fabs((perf[p + 1].depth - perf[p].depth) / maxlen));
414  if (seg_num == 0) continue;
415  seg_len = (perf[p + 1].depth - perf[p].depth) / seg_num;
416  }
417  OCP_USI n = perf[p].location;
418  perf[p].P = bhp + dG[p];
419  OCP_DBL Pperf = perf[p].P;
420  OCP_DBL Ptmp = Pperf;
421 
422  USI pvtnum = myBulk.PVTNUM[n];
423  for (USI i = 0; i < seg_num; i++) {
424  Ptmp += myBulk.flashCal[pvtnum]->RhoPhase(
425  Ptmp, 0, opt.injTemp, opt.injZi.data(), opt.injProdPhase) *
426  GRAVITY_FACTOR * seg_len;
427  }
428  dGperf[p] = Ptmp - Pperf;
429  }
430  dG[numPerf - 1] = dGperf[numPerf - 1];
431  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
432  dG[p] = dG[p + 1] + dGperf[p];
433  }
434  }
435 }
436 
437 // Use transj
438 void Well::CalProddG01(const Bulk& myBulk)
439 {
440  OCP_FUNCNAME;
441 
442  const OCP_DBL maxlen = 10;
443  USI seg_num = 0;
444  OCP_DBL seg_len = 0;
445  vector<OCP_DBL> dGperf(numPerf, 0);
446  vector<OCP_DBL> tmpNi(numCom, 0);
447  OCP_DBL rhotmp, qtacc, rhoacc;
448 
449  if (depth <= perf.front().depth) {
450  // Well is higher
451  for (OCP_INT p = numPerf - 1; p >= 0; p--) {
452  if (p == 0) {
453  seg_num = ceil(fabs((perf[0].depth - depth) / maxlen));
454  if (seg_num == 0) continue;
455  seg_len = (perf[0].depth - depth) / seg_num;
456  } else {
457  seg_num = ceil(fabs((perf[p].depth - perf[p - 1].depth) / maxlen));
458  if (seg_num == 0) continue;
459  seg_len = (perf[p].depth - perf[p - 1].depth) / seg_num;
460  }
461  OCP_USI n = perf[p].location;
462  perf[p].P = bhp + dG[p];
463  OCP_DBL Pperf = perf[p].P;
464  OCP_DBL Ptmp = Pperf;
465 
466  // fill(tmpNi.begin(), tmpNi.end(), 0.0);
467  for (USI j = 0; j < numPhase; j++) {
468  const OCP_USI n_np_j = n * numPhase + j;
469  if (!myBulk.phaseExist[n_np_j]) continue;
470  for (USI k = 0; k < numCom; k++) {
471  tmpNi[k] += (myBulk.P[n] - perf[p].P) * perf[p].transj[j] *
472  myBulk.xi[n_np_j] * myBulk.xij[n_np_j * numCom + k];
473  }
474  }
475  OCP_DBL tmpSum = Dnorm1(numCom, &tmpNi[0]);
476  if (tmpSum < TINY) {
477  for (USI i = 0; i < numCom; i++) {
478  tmpNi[i] = myBulk.Ni[n * numCom + i];
479  }
480  }
481 
482  USI pvtnum = myBulk.PVTNUM[n];
483  for (USI i = 0; i < seg_num; i++) {
484  qtacc = rhoacc = 0;
485  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T[n], tmpNi.data());
486  for (USI j = 0; j < myBulk.numPhase; j++) {
487  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
488  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
489  qtacc += myBulk.flashCal[pvtnum]->vj[j];
490  rhoacc += myBulk.flashCal[pvtnum]->vj[j] * rhotmp;
491  }
492  }
493  Ptmp -= rhoacc / qtacc * GRAVITY_FACTOR * seg_len;
494  }
495  dGperf[p] = Pperf - Ptmp;
496  }
497  dG[0] = dGperf[0];
498  for (USI p = 1; p < numPerf; p++) {
499  dG[p] = dG[p - 1] + dGperf[p];
500  }
501  } else if (depth >= perf[numPerf - 1].depth) {
502  // Well is lower
503  for (USI p = 0; p < numPerf; p++) {
504  if (p == numPerf - 1) {
505  seg_num = ceil(fabs((depth - perf[numPerf - 1].depth) / maxlen));
506  if (seg_num == 0) continue;
507  seg_len = (depth - perf[numPerf - 1].depth) / seg_num;
508  } else {
509  seg_num = ceil(fabs((perf[p + 1].depth - perf[p].depth) / maxlen));
510  if (seg_num == 0) continue;
511  seg_len = (perf[p + 1].depth - perf[p].depth) / seg_num;
512  }
513  OCP_USI n = perf[p].location;
514  perf[p].P = bhp + dG[p];
515  OCP_DBL Pperf = perf[p].P;
516  OCP_DBL Ptmp = Pperf;
517 
518  // fill(tmpNi.begin(), tmpNi.end(), 0.0);
519  for (USI j = 0; j < numPhase; j++) {
520  const OCP_USI n_np_j = n * numPhase + j;
521  if (!myBulk.phaseExist[n_np_j]) continue;
522  for (USI k = 0; k < numCom; k++) {
523  tmpNi[k] += (myBulk.P[n] - perf[p].P) * perf[p].transj[j] *
524  myBulk.xi[n_np_j] * myBulk.xij[n_np_j * numCom + k];
525  }
526  }
527  OCP_DBL tmpSum = Dnorm1(numCom, &tmpNi[0]);
528  if (tmpSum < TINY) {
529  for (USI i = 0; i < numCom; i++) {
530  tmpNi[i] = myBulk.Ni[n * numCom + i];
531  }
532  }
533 
534  USI pvtnum = myBulk.PVTNUM[n];
535  for (USI i = 0; i < seg_num; i++) {
536  qtacc = rhoacc = 0;
537  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T[n], tmpNi.data());
538  for (USI j = 0; j < myBulk.numPhase; j++) {
539  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
540  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
541  qtacc += myBulk.flashCal[pvtnum]->vj[j];
542  rhoacc += myBulk.flashCal[pvtnum]->vj[j] * rhotmp;
543  }
544  }
545  Ptmp += rhoacc / qtacc * GRAVITY_FACTOR * seg_len;
546  }
547  dGperf[p] = Ptmp - Pperf;
548  }
549  dG[numPerf - 1] = dGperf[numPerf - 1];
550  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
551  dG[p] = dG[p + 1] + dGperf[p];
552  }
553  }
554 }
555 
556 // Use bulk
557 void Well::CalProddG02(const Bulk& myBulk)
558 {
559  OCP_FUNCNAME;
560 
561  const OCP_DBL maxlen = 10;
562  USI seg_num = 0;
563  OCP_DBL seg_len = 0;
564  vector<OCP_DBL> dGperf(numPerf, 0);
565  vector<OCP_DBL> tmpNi(numCom, 0);
566  OCP_DBL rhotmp, qtacc, rhoacc;
567 
568  if (depth <= perf.front().depth) {
569  // Well is higher
570  for (OCP_INT p = numPerf - 1; p >= 0; p--) {
571  if (p == 0) {
572  seg_num = ceil(fabs((perf[0].depth - depth) / maxlen));
573  if (seg_num == 0) continue;
574  seg_len = (perf[0].depth - depth) / seg_num;
575  } else {
576  seg_num = ceil(fabs((perf[p].depth - perf[p - 1].depth) / maxlen));
577  if (seg_num == 0) continue;
578  seg_len = (perf[p].depth - perf[p - 1].depth) / seg_num;
579  }
580  OCP_USI n = perf[p].location;
581  perf[p].P = bhp + dG[p];
582  OCP_DBL Pperf = perf[p].P;
583  OCP_DBL Ptmp = Pperf;
584 
585  fill(tmpNi.begin(), tmpNi.end(), 0.0);
586  for (USI j = 0; j < numPhase; j++) {
587  OCP_USI id = n * numPhase + j;
588  if (!myBulk.phaseExist[id]) continue;
589  for (USI k = 0; k < numCom; k++) {
590  tmpNi[k] += (perf[p].transj[j] > 0) * myBulk.xi[id] *
591  myBulk.xij[id * numCom + k];
592  }
593  }
594  OCP_DBL tmpSum = Dnorm1(numCom, &tmpNi[0]);
595  if (tmpSum < TINY) {
596  for (USI i = 0; i < numCom; i++) {
597  tmpNi[i] = myBulk.Ni[n * numCom + i];
598  }
599  }
600 
601  USI pvtnum = myBulk.PVTNUM[n];
602  for (USI i = 0; i < seg_num; i++) {
603  qtacc = rhoacc = 0;
604  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T[n], tmpNi.data());
605  for (USI j = 0; j < myBulk.numPhase; j++) {
606  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
607  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
608  qtacc += myBulk.flashCal[pvtnum]->vj[j];
609  rhoacc += myBulk.flashCal[pvtnum]->vj[j] * rhotmp;
610  }
611  }
612  Ptmp -= rhoacc / qtacc * GRAVITY_FACTOR * seg_len;
613  }
614  dGperf[p] = Pperf - Ptmp;
615  }
616  dG[0] = dGperf[0];
617  for (USI p = 1; p < numPerf; p++) {
618  dG[p] = dG[p - 1] + dGperf[p];
619  }
620  } else if (depth >= perf[numPerf - 1].depth) {
621  // Well is lower
622  for (USI p = 0; p < numPerf; p++) {
623  if (p == numPerf - 1) {
624  seg_num = ceil(fabs((depth - perf[numPerf - 1].depth) / maxlen));
625  if (seg_num == 0) continue;
626  seg_len = (depth - perf[numPerf - 1].depth) / seg_num;
627  } else {
628  seg_num = ceil(fabs((perf[p + 1].depth - perf[p].depth) / maxlen));
629  if (seg_num == 0) continue;
630  seg_len = (perf[p + 1].depth - perf[p].depth) / seg_num;
631  }
632  OCP_USI n = perf[p].location;
633  perf[p].P = bhp + dG[p];
634  OCP_DBL Pperf = perf[p].P;
635  OCP_DBL Ptmp = Pperf;
636 
637  fill(tmpNi.begin(), tmpNi.end(), 0.0);
638  for (USI j = 0; j < numPhase; j++) {
639  OCP_USI id = n * numPhase + j;
640  if (!myBulk.phaseExist[id]) continue;
641  for (USI k = 0; k < numCom; k++) {
642  tmpNi[k] += (perf[p].transj[j] > 0) * myBulk.xi[id] *
643  myBulk.xij[id * numCom + k];
644  }
645  }
646  OCP_DBL tmpSum = Dnorm1(numCom, &tmpNi[0]);
647  if (tmpSum < TINY) {
648  for (USI i = 0; i < numCom; i++) {
649  tmpNi[i] = myBulk.Ni[n * numCom + i];
650  }
651  }
652 
653  USI pvtnum = myBulk.PVTNUM[n];
654  for (USI i = 0; i < seg_num; i++) {
655  qtacc = rhoacc = 0;
656  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T[n], tmpNi.data());
657  for (USI j = 0; j < myBulk.numPhase; j++) {
658  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
659  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
660  qtacc += myBulk.flashCal[pvtnum]->vj[j];
661  rhoacc += myBulk.flashCal[pvtnum]->vj[j] * rhotmp;
662  }
663  }
664  Ptmp += rhoacc / qtacc * GRAVITY_FACTOR * seg_len;
665  }
666  dGperf[p] = Ptmp - Pperf;
667  }
668  dG[numPerf - 1] = dGperf[numPerf - 1];
669  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
670  dG[p] = dG[p + 1] + dGperf[p];
671  }
672  }
673 }
674 
675 // Use qi_lbmol
676 void Well::CalProddG(const Bulk& myBulk)
677 {
678  OCP_FUNCNAME;
679 
680  const OCP_DBL maxlen = 5;
681  USI seg_num = 0;
682  OCP_DBL seg_len = 0;
683  vector<OCP_DBL> tmpNi(numCom, 0);
684  vector<OCP_DBL> dGperf(numPerf, 0);
685  OCP_DBL qtacc = 0;
686  OCP_DBL rhoacc = 0;
687  OCP_DBL rhotmp = 0;
688 
689  if (depth <= perf.front().depth) {
690  // Well is higher
691 
692  // check qi_lbmol ---- test
693  if (perf[numPerf - 1].state == CLOSE) {
694  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
695  if (perf[p].state == OPEN) {
696  for (USI i = 0; i < numCom; i++) {
697  perf[numPerf - 1].qi_lbmol[i] = perf[p].qi_lbmol[i];
698  }
699  break;
700  }
701  }
702  }
703 
704  for (OCP_INT p = numPerf - 1; p >= 0; p--) {
705 
706  if (p == 0) {
707  seg_num = ceil(fabs((perf[0].depth - depth) / maxlen));
708  if (seg_num == 0) continue;
709  seg_len = (perf[0].depth - depth) / seg_num;
710  } else {
711  seg_num = ceil(fabs((perf[p].depth - perf[p - 1].depth) / maxlen));
712  if (seg_num == 0) continue;
713  seg_len = (perf[p].depth - perf[p - 1].depth) / seg_num;
714  }
715 
716  OCP_USI n = perf[p].location;
717  perf[p].P = bhp + dG[p];
718  OCP_DBL Pperf = perf[p].P;
719  OCP_DBL Ptmp = Pperf;
720 
721  USI pvtnum = myBulk.PVTNUM[n];
722  // tmpNi.assign(numCom, 0);
723  // for (OCP_INT p1 = numPerf - 1; p1 >= p; p1--) {
724  // for (USI i = 0; i < numCom; i++) {
725  // tmpNi[i] += perf[p1].qi_lbmol[i];
726  // }
727  //}
728 
729  for (USI i = 0; i < numCom; i++) {
730  tmpNi[i] += perf[p].qi_lbmol[i];
731  }
732 
733  // check tmpNi
734  for (auto& v : tmpNi) {
735  v = fabs(v);
736  }
737 
738  for (USI k = 0; k < seg_num; k++) {
739  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T[n], tmpNi.data());
740  for (USI j = 0; j < myBulk.numPhase; j++) {
741  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
742  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
743  qtacc += myBulk.flashCal[pvtnum]->vj[j] / seg_num;
744  rhoacc += myBulk.flashCal[pvtnum]->vj[j] * rhotmp *
745  GRAVITY_FACTOR / seg_num;
746 #ifdef DEBUG
747  if (rhotmp <= 0 || !isfinite(rhotmp)) {
748  OCP_ABORT("Wrong rho " + to_string(rhotmp));
749  }
750 #endif // DEBUG
751  }
752  }
753  Ptmp -= rhoacc / qtacc * seg_len;
754  }
755  dGperf[p] = Pperf - Ptmp;
756  }
757  dG[0] = dGperf[0];
758  for (USI p = 1; p < numPerf; p++) {
759  dG[p] = dG[p - 1] + dGperf[p];
760  }
761  } else if (depth >= perf.back().depth) {
762  // Well is lower
763 
764  // check qi_lbmol ---- test
765  if (perf[0].state == CLOSE) {
766  for (USI p = 1; p <= numPerf; p++) {
767  if (perf[p].state == OPEN) {
768  for (USI i = 0; i < numCom; i++) {
769  perf[numPerf - 1].qi_lbmol[i] = perf[p].qi_lbmol[i];
770  }
771  break;
772  }
773  }
774  }
775 
776  for (USI p = 0; p < numPerf; p++) {
777  if (p == numPerf - 1) {
778  seg_num = ceil(fabs((depth - perf[numPerf - 1].depth) / maxlen));
779  if (seg_num == 0) continue;
780  seg_len = (depth - perf[numPerf - 1].depth) / seg_num;
781  } else {
782  seg_num = ceil(fabs((perf[p + 1].depth - perf[p].depth) / maxlen));
783  if (seg_num == 0) continue;
784  seg_len = (perf[p + 1].depth - perf[p].depth) / seg_num;
785  }
786 
787  OCP_USI n = perf[p].location;
788  perf[p].P = bhp + dG[p];
789  OCP_DBL Pperf = perf[p].P;
790  OCP_DBL Ptmp = Pperf;
791 
792  USI pvtnum = myBulk.PVTNUM[n];
793  fill(tmpNi.begin(), tmpNi.end(), 0.0);
794  for (OCP_INT p1 = numPerf - 1; p1 - p >= 0; p1--) {
795  for (USI i = 0; i < numCom; i++) {
796  tmpNi[i] += perf[p1].qi_lbmol[i];
797  }
798  }
799 
800  // check tmpNi
801  for (auto& v : tmpNi) {
802  v = fabs(v);
803  }
804 
805  for (USI k = 0; k < seg_num; k++) {
806  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T[n], tmpNi.data());
807  for (USI j = 0; j < numPhase; j++) {
808  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
809  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
810  qtacc += myBulk.flashCal[pvtnum]->vj[j] / seg_num;
811  rhoacc += myBulk.flashCal[pvtnum]->vj[j] * rhotmp *
812  GRAVITY_FACTOR / seg_num;
813  }
814  }
815  Ptmp += rhoacc / qtacc * seg_len;
816  }
817  dGperf[p] = Ptmp - Pperf;
818  }
819  dG[numPerf - 1] = dGperf[numPerf - 1];
820  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
821  dG[p] = dG[p + 1] + dGperf[p];
822  }
823 
824  } else {
825  OCP_ABORT("Wrong well position!");
826  }
827 }
828 
829 void Well::CalProdWeight(const Bulk& myBulk) const
830 {
831 
832  if (opt.type == PROD) {
833  // in some cases, qi_lbmol may be zero, so use other methods
834  OCP_DBL qt = 0;
835  OCP_BOOL flag = OCP_TRUE;
836  for (USI i = 0; i < myBulk.numCom; i++) {
837  qt += qi_lbmol[i];
838  if (qi_lbmol[i] < 0) flag = OCP_FALSE;
839  }
840  if (qt > TINY && flag) {
841  flashCal[0]->CalProdWeight(Psurf, Tsurf, &qi_lbmol[0], opt.prodPhaseWeight,
842  prodWeight);
843  } else {
844  vector<OCP_DBL> tmpNi(numCom, 0);
845  for (USI p = 0; p < numPerf; p++) {
846  OCP_USI n = perf[p].location;
847 
848  for (USI j = 0; j < numPhase; j++) {
849  OCP_USI id = n * numPhase + j;
850  if (!myBulk.phaseExist[id]) continue;
851  for (USI k = 0; k < numCom; k++) {
852  tmpNi[k] += perf[p].transj[j] * myBulk.xi[id] *
853  myBulk.xij[id * numCom + k];
854  }
855  }
856  }
857  qt = Dnorm1(numCom, &tmpNi[0]);
858  flashCal[0]->CalProdWeight(Psurf, Tsurf, &tmpNi[0], opt.prodPhaseWeight,
859  prodWeight);
860  }
861  }
862 }
863 
864 void Well::CalReInjFluid(const Bulk& myBulk, vector<OCP_DBL>& myZi)
865 {
866  CalTrans(myBulk);
867 
868  for (USI p = 0; p < numPerf; p++) {
869  OCP_USI n = perf[p].location;
870 
871  for (USI j = 0; j < numPhase; j++) {
872  OCP_USI id = n * numPhase + j;
873  if (!myBulk.phaseExist[id]) continue;
874  for (USI k = 0; k < numCom; k++) {
875  myZi[k] +=
876  perf[p].transj[j] * myBulk.xi[id] * myBulk.xij[id * numCom + k];
877  }
878  }
879  }
880 }
881 
882 void Well::CorrectBHP()
883 {
884  if (opt.type == PROD && opt.optMode == BHP_MODE) {
885  bhp = opt.minBHP;
886  } else if (opt.type == INJ && opt.optMode == BHP_MODE) {
887  bhp = opt.maxBHP;
888  }
889 }
890 
893 void Well::CheckOptMode(const Bulk& myBulk)
894 {
895  OCP_FUNCNAME;
896  if (opt.initOptMode == BHP_MODE) {
897  if (opt.type == INJ) {
898  OCP_DBL q = CalInjRateMaxBHP(myBulk);
899  // for INJ well, maxRate has been switch to lbmols
900  OCP_DBL tarRate = opt.maxRate;
901  if (opt.reInj) {
902  if (opt.reInjPhase == GAS)
903  tarRate = WGIR;
904  else if (opt.reInjPhase == WATER)
905  tarRate = WWIR;
906  }
907  if (q > tarRate) {
908  opt.optMode = RATE_MODE;
909  } else {
910  opt.optMode = BHP_MODE;
911  bhp = opt.maxBHP;
912  }
913  } else {
914  opt.optMode = BHP_MODE;
915  bhp = opt.minBHP;
916  }
917  } else {
918  if (opt.type == INJ) {
919  OCP_DBL q = CalInjRateMaxBHP(myBulk);
920  // for INJ well, maxRate has been switch to lbmols
921  OCP_DBL tarRate = opt.maxRate;
922  if (opt.reInj) {
923  if (opt.reInjPhase == GAS)
924  tarRate = WGIR;
925  else if (opt.reInjPhase == WATER)
926  tarRate = WWIR;
927  }
928 
929  if (q > tarRate) {
930  opt.optMode = opt.initOptMode;
931  } else {
932  opt.optMode = BHP_MODE;
933  bhp = opt.maxBHP;
934  }
935  } else {
936  OCP_DBL q = CalProdRateMinBHP(myBulk);
937  // cout << q << endl;
938  if (q > opt.maxRate) {
939  opt.optMode = opt.initOptMode;
940  } else {
941  opt.optMode = BHP_MODE;
942  bhp = opt.minBHP;
943  }
944  }
945  }
946 }
947 
948 OCP_INT Well::CheckP(const Bulk& myBulk)
949 {
950  OCP_FUNCNAME;
951  // 0 : all correct
952  // 1 : negative P
953  // 2 : outlimited P
954  // 3 : crossflow happens
955 
956  if (bhp < 0) {
957  cout << "### WARNING: Well " << name << " BHP = " << bhp << endl;
958  return WELL_NEGATIVE_PRESSURE;
959  }
960  for (USI p = 0; p < numPerf; p++) {
961  if (perf[p].state == OPEN && perf[p].P < 0) {
962 #ifdef DEBUG
963  cout << "### WARNING: Well " << name << " Perf[" << p
964  << "].P = " << perf[p].P << endl;
965 #endif // DEBUG
966  return WELL_NEGATIVE_PRESSURE;
967  }
968  }
969 
970  if (opt.type == INJ) {
971  if (opt.optMode != BHP_MODE && bhp > opt.maxBHP) {
972 #if _DEBUG
973  cout << "### WARNING: Well " << name << " switch to BHPMode" << endl;
974 #endif
975  opt.optMode = BHP_MODE;
976  bhp = opt.maxBHP;
977  return WELL_SWITCH_TO_BHPMODE;
978  }
979  } else {
980  if (opt.optMode != BHP_MODE && bhp < opt.minBHP) {
981 #if _DEBUG
982  cout << "### WARNING: Well " << name << " switch to BHPMode" << endl;
983 #endif
984  opt.optMode = BHP_MODE;
985  bhp = opt.minBHP;
986  return WELL_SWITCH_TO_BHPMODE;
987  }
988  }
989 
990  return CheckCrossFlow(myBulk);
991 }
992 
993 OCP_INT Well::CheckCrossFlow(const Bulk& myBulk)
994 {
995  OCP_FUNCNAME;
996 
997  OCP_USI k;
998  OCP_BOOL flagC = OCP_TRUE;
999 
1000  if (opt.type == PROD) {
1001  for (USI p = 0; p < numPerf; p++) {
1002  k = perf[p].location;
1003  OCP_DBL minP = myBulk.P[k];
1004  if (perf[p].state == OPEN && minP < perf[p].P) {
1005  cout << std::left << std::setw(12) << name << " "
1006  << "Well P = " << perf[p].P << ", "
1007  << "Bulk P = " << minP << endl;
1008  perf[p].state = CLOSE;
1009  perf[p].multiplier = 0;
1010  flagC = OCP_FALSE;
1011  break;
1012  } else if (perf[p].state == CLOSE && minP > perf[p].P) {
1013  perf[p].state = OPEN;
1014  perf[p].multiplier = 1;
1015  }
1016  }
1017  } else {
1018  for (USI p = 0; p < numPerf; p++) {
1019  k = perf[p].location;
1020  if (perf[p].state == OPEN && myBulk.P[k] > perf[p].P) {
1021  cout << std::left << std::setw(12) << name << " "
1022  << "Well P = " << perf[p].P << ", "
1023  << "Bulk P = " << myBulk.P[k] << endl;
1024  perf[p].state = CLOSE;
1025  perf[p].multiplier = 0;
1026  flagC = OCP_FALSE;
1027  break;
1028  } else if (perf[p].state == CLOSE && myBulk.P[k] < perf[p].P) {
1029  perf[p].state = OPEN;
1030  perf[p].multiplier = 1;
1031  }
1032  }
1033  }
1034 
1035  OCP_BOOL flag = OCP_FALSE;
1036  // check well -- if all perf are closed, open the depthest perf
1037  for (USI p = 0; p < numPerf; p++) {
1038  if (perf[p].state == OPEN) {
1039  flag = OCP_TRUE;
1040  break;
1041  }
1042  }
1043 
1044  if (!flag) {
1045  // open the deepest perf
1046  perf.back().state = OPEN;
1047  perf.back().multiplier = 1;
1048  cout << "### WARNING: All perfs of " << name
1049  << " are closed! Open the last perf!\n";
1050  }
1051 
1052  if (!flagC) {
1053  // if crossflow happens, then corresponding perforation will be closed,
1054  // the multiplier of perforation will be set to zero, so trans of well
1055  // should be recalculated!
1056  //
1057  // dG = ldG;
1058  CalTrans(myBulk);
1059  // CalFlux(myBulk);
1060  // CaldG(myBulk);
1061  // CheckOptMode(myBulk);
1062  return WELL_CROSSFLOW;
1063  }
1064 
1065  return WELL_SUCCESS;
1066 }
1067 
1068 void Well::ShowPerfStatus(const Bulk& myBulk) const
1069 {
1070  OCP_FUNCNAME;
1071 
1072  cout << fixed;
1073  cout << "----------------------------" << endl;
1074  cout << name << ": " << opt.optMode << " " << setprecision(3) << bhp << endl;
1075  for (USI p = 0; p < numPerf; p++) {
1076  vector<OCP_DBL> Qitmp(perf[p].qi_lbmol);
1077  // OCP_DBL qt = Dnorm1(myBulk.numCom, &Qitmp[0]);
1078  OCP_USI n = perf[p].location;
1079  cout << setw(3) << p << " " << perf[p].state << " " << setw(6)
1080  << perf[p].location << " " << setw(2) << perf[p].I + 1 << " " << setw(2)
1081  << perf[p].J + 1 << " " << setw(2) << perf[p].K + 1 << " " << setw(10)
1082  << setprecision(6) << perf[p].WI << " " // ccf
1083  << setprecision(3) << perf[p].radius << " " // ccf
1084  << setw(8) << setprecision(4) << perf[p].kh << " " // kh
1085  << setw(8) << setprecision(2) << perf[p].depth << " " // depth
1086  << setprecision(3) << perf[p].P << " " // Pp
1087  << setw(10) << setprecision(3) << myBulk.P[n] << " " // Pb
1088  << setw(6) << setprecision(3) << dG[p] << " " // dG
1089  << setw(8) << perf[p].qi_lbmol[myBulk.numCom - 1] << " " << setw(6)
1090  << setprecision(6) << myBulk.S[n * myBulk.numPhase + 0] << " " << setw(6)
1091  << setprecision(6) << myBulk.S[n * myBulk.numPhase + 1] << " " << setw(6)
1092  << setprecision(6) << myBulk.S[n * myBulk.numPhase + 2] << endl;
1093  }
1094 }
1095 
1096 void Well::AssembleMatReinjection_IMPEC(const Bulk& myBulk,
1097  LinearSystem& myLS,
1098  const OCP_DBL& dt,
1099  const vector<Well>& allWell,
1100  const vector<USI>& injId) const
1101 {
1102  // find Open injection well under Rate control
1103  vector<OCP_USI> tarId;
1104  for (auto& w : injId) {
1105  if (allWell[w].IsOpen() && allWell[w].opt.optMode != BHP_MODE)
1106  tarId.push_back(allWell[w].wOId + myBulk.numBulk);
1107  }
1108 
1109  USI tlen = tarId.size();
1110  if (tlen > 0) {
1111  // All inj well has the same factor
1112  const OCP_DBL factor = allWell[injId[0]].opt.reInjFactor * dt;
1113  const OCP_USI prodId = wOId + myBulk.numBulk;
1114  OCP_USI n, bId;
1115  OCP_DBL tmp, valb;
1116  OCP_DBL valw = 0;
1117  OCP_DBL rhsw = 0;
1118  for (USI p = 0; p < numPerf; p++) {
1119  n = perf[p].location;
1120  valb = 0;
1121 
1122  for (USI j = 0; j < numPhase; j++) {
1123  bId = n * numPhase + j;
1124  if (myBulk.phaseExist[bId]) {
1125  tmp = perf[p].transj[j] * myBulk.xi[bId];
1126  valb += tmp;
1127  rhsw += tmp * (myBulk.Pc[bId] - dG[p]);
1128  }
1129  }
1130  valb *= factor;
1131  for (USI t = 0; t < tlen; t++) {
1132  myLS.NewOffDiag(tarId[t], n, -valb);
1133  }
1134  valw += valb;
1135  }
1136  rhsw *= factor;
1137  // rhs and prod well
1138  for (USI t = 0; t < tlen; t++) {
1139  myLS.AddRhs(tarId[t], rhsw);
1140  myLS.NewOffDiag(tarId[t], prodId, valw);
1141  }
1142  }
1143 }
1144 
1145 void Well::AssembleMatReinjection_FIM(const Bulk& myBulk,
1146  LinearSystem& myLS,
1147  const OCP_DBL& dt,
1148  const vector<Well>& allWell,
1149  const vector<USI>& injId) const
1150 {
1151  // find Open injection well under Rate control
1152  vector<OCP_USI> tarId;
1153  for (auto& w : injId) {
1154  if (allWell[w].IsOpen() && allWell[w].opt.optMode != BHP_MODE)
1155  tarId.push_back(allWell[w].wOId + myBulk.numBulk);
1156  }
1157 
1158  USI tlen = tarId.size();
1159  if (tlen > 0) {
1160  // All inj well has the same factor
1161  const OCP_DBL factor = allWell[injId[0]].opt.reInjFactor;
1162  const OCP_USI prodId = wOId + myBulk.numBulk;
1163 
1164  const USI ncol = numCom + 1;
1165  const USI ncol2 = numPhase * numCom + numPhase;
1166  const USI bsize = ncol * ncol;
1167  const USI bsize2 = ncol * ncol2;
1168 
1169  OCP_DBL xij, xi, mu, muP, xiP, dP, transIJ, tmp;
1170  OCP_USI n_np_j;
1171 
1172  vector<OCP_DBL> bmat(bsize, 0);
1173  vector<OCP_DBL> bmat2(bsize, 0);
1174  vector<OCP_DBL> tmpMat(bsize, 0);
1175  vector<OCP_DBL> dQdXpB(bsize, 0);
1176  vector<OCP_DBL> dQdXpW(bsize, 0);
1177  vector<OCP_DBL> dQdXsB(bsize2, 0);
1178 
1179  for (USI p = 0; p < numPerf; p++) {
1180  const OCP_USI n = perf[p].location;
1181  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
1182  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
1183  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
1184 
1185  for (USI j = 0; j < numPhase; j++) {
1186  n_np_j = n * numPhase + j;
1187  if (!myBulk.phaseExist[n_np_j]) continue;
1188 
1189  dP = myBulk.Pj[n_np_j] - bhp - dG[p];
1190  xi = myBulk.xi[n_np_j];
1191  mu = myBulk.mu[n_np_j];
1192  muP = myBulk.muP[n_np_j];
1193  xiP = myBulk.xiP[n_np_j];
1194 
1195  for (USI i = 0; i < numCom; i++) {
1196  xij = myBulk.xij[n_np_j * numCom + i];
1197  // dQ / dP
1198  transIJ = perf[p].transj[j] * xi * xij;
1199  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu) +
1200  dP * perf[p].transj[j] * xij * xiP;
1201  dQdXpW[(i + 1) * ncol] += -transIJ;
1202 
1203  // dQ / dS
1204  for (USI k = 0; k < numPhase; k++) {
1205  tmp = CONV1 * perf[p].WI * perf[p].multiplier * dP / mu * xi *
1206  xij * myBulk.dKr_dS[n_np_j * numPhase + k];
1207  // capillary pressure
1208  tmp += transIJ * myBulk.dPcj_dS[n_np_j * numPhase + k];
1209  dQdXsB[(i + 1) * ncol2 + k] += tmp;
1210  }
1211  // dQ / dCij
1212  for (USI k = 0; k < numCom; k++) {
1213  tmp = dP * perf[p].transj[j] * xij *
1214  (myBulk.xix[n_np_j * numCom + k] -
1215  xi / mu * myBulk.mux[n_np_j * numCom + k]);
1216  if (k == i) {
1217  tmp += perf[p].transj[j] * xi * dP;
1218  }
1219  dQdXsB[(i + 1) * ncol2 + numPhase + j * numCom + k] += tmp;
1220  }
1221  }
1222  }
1223 
1224  // for Prod Well
1225  for (USI i = 0; i < numCom; i++) {
1226  tmpMat[0] += dQdXpW[(i + 1) * ncol] * factor;
1227  }
1228 
1229  // for Perf(bulk) of Prod Well
1230  bmat = dQdXpB;
1231  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2],
1232  1, bmat.data());
1233  fill(bmat2.begin(), bmat2.end(), 0.0);
1234  for (USI i = 0; i < numCom; i++) {
1235  // becareful '-' before factor
1236  Daxpy(ncol, -factor, bmat.data() + (i + 1) * ncol, bmat2.data());
1237  }
1238 
1239  // Insert bulk val into equations in ascending order
1240  for (USI t = 0; t < tlen; t++) {
1241  myLS.NewOffDiag(tarId[t], n, bmat2);
1242  }
1243  }
1244  // prod well
1245  for (USI t = 0; t < tlen; t++) {
1246  myLS.NewOffDiag(tarId[t], prodId, tmpMat);
1247  }
1248  }
1249 }
1250 
1251 void Well::SetPolyhedronWell(const Grid& myGrid, OCPpolyhedron& mypol)
1252 {
1253  // set a virtual point
1254  mypol.numPoints = numPerf + 1;
1255  mypol.Points.resize(mypol.numPoints);
1256 
1257  OCP_USI k;
1258  Point3D tmpP;
1259 
1260  for (USI p = 0; p < numPerf; p++) {
1261  tmpP.Reset();
1262  k = myGrid.map_Act2All[perf[p].location];
1263  for (USI i = 0; i < myGrid.polyhedronGrid[k].numPoints; i++) {
1264  tmpP += myGrid.polyhedronGrid[k].Points[i];
1265  }
1266  tmpP /= myGrid.polyhedronGrid[k].numPoints;
1267  mypol.Points[p + 1] = tmpP;
1268  }
1269 
1270  // Set virtual perf
1271  mypol.Points[0] = mypol.Points[1];
1272  mypol.Points[0].z /= 2;
1273 }
1274 
1275 /*----------------------------------------------------------------------------*/
1276 /* Brief Change History of This File */
1277 /*----------------------------------------------------------------------------*/
1278 /* Author Date Actions */
1279 /*----------------------------------------------------------------------------*/
1280 /* Shizhe Li Oct/01/2021 Create file */
1281 /* Chensong Zhang Oct/15/2021 Format file */
1282 /*----------------------------------------------------------------------------*/
DaABpbC
void DaABpbC(const int &m, const int &n, const int &k, const double &alpha, const double *A, const double *B, const double &beta, double *C)
Computes C' = alpha B'A' + beta C', all matrices are column-major.
Definition: DenseMat.cpp:75
Daxpy
void Daxpy(const int &n, const double &alpha, const double *x, double *y)
Constant times a vector plus a vector.
Definition: DenseMat.cpp:68
Dnorm1
double Dnorm1(const int &N, double *x)
Computes the L1-norm of a vector.
Definition: DenseMat.cpp:49
X_DIRECTION
const USI X_DIRECTION
x-direction
Definition: OCPConst.hpp:138
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:38
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:23
GAS
const USI GAS
Fluid type = gas.
Definition: OCPConst.hpp:92
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:27
CONV1
const OCP_DBL CONV1
1 bbl = CONV1 ft3
Definition: OCPConst.hpp:63
GRAVITY_FACTOR
const OCP_DBL GRAVITY_FACTOR
0.00694444 ft2 psi / lb
Definition: OCPConst.hpp:56
WATER
const USI WATER
Fluid type = water.
Definition: OCPConst.hpp:93
OPEN
const OCP_BOOL OPEN
Well type = open.
Definition: OCPConst.hpp:126
Y_DIRECTION
const USI Y_DIRECTION
y-direction
Definition: OCPConst.hpp:139
RATE_MODE
const USI RATE_MODE
Well option = fixed total rate???
Definition: OCPConst.hpp:130
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:25
Z_DIRECTION
const USI Z_DIRECTION
z-direction
Definition: OCPConst.hpp:140
PROD
const USI PROD
Well type = producer.
Definition: OCPConst.hpp:123
PI
const OCP_DBL PI
Pi.
Definition: OCPConst.hpp:39
OCP_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:26
INJ
const USI INJ
Well type = injector.
Definition: OCPConst.hpp:122
CLOSE
const OCP_BOOL CLOSE
Well type = closed.
Definition: OCPConst.hpp:127
BHP_MODE
const USI BHP_MODE
Well option = fixed bottom-hole-pressure.
Definition: OCPConst.hpp:135
OCP_BOOL
unsigned int OCP_BOOL
OCP_BOOL in OCP.
Definition: OCPConst.hpp:29
CONV2
const OCP_DBL CONV2
Darcy constant in field unit.
Definition: OCPConst.hpp:64
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
Well.hpp
Well class declaration.
Bulk
Physical information of each active reservoir bulk.
Definition: Bulk.hpp:104
Bulk::mu
vector< OCP_DBL > mu
Viscosity of phase: numPhase*numBulk.
Definition: Bulk.hpp:320
Bulk::Pc
vector< OCP_DBL > Pc
Capillary pressure of phase: numPhase*numBulk.
Definition: Bulk.hpp:312
Bulk::numPhase
USI numPhase
Number of phase.
Definition: Bulk.hpp:154
Bulk::T
vector< OCP_DBL > T
Temperature: numBulk.
Definition: Bulk.hpp:308
Bulk::S
vector< OCP_DBL > S
Saturation of phase: numPhase*numBulk.
Definition: Bulk.hpp:314
Bulk::oil
OCP_BOOL oil
If OCP_TRUE, oil phase could exist.
Definition: Bulk.hpp:222
Bulk::dSec_dPri
vector< OCP_DBL > dSec_dPri
d Secondary variable / d Primary variable.
Definition: Bulk.hpp:486
Bulk::ntg
vector< OCP_DBL > ntg
net to gross of bulk.
Definition: Bulk.hpp:243
Bulk::dy
vector< OCP_DBL > dy
Size of cell in y-direction: activeGridNum.
Definition: Bulk.hpp:238
Bulk::PVTNUM
vector< USI > PVTNUM
Identify PVT region in black-oil model: numBulk.
Definition: Bulk.hpp:191
Bulk::gas
OCP_BOOL gas
If OCP_TRUE, gas phase could exist.
Definition: Bulk.hpp:223
Bulk::xij
vector< OCP_DBL > xij
Nij / Nj: numPhase*numCom*numBulk.
Definition: Bulk.hpp:317
Bulk::kr
vector< OCP_DBL > kr
Relative permeability of phase: numPhase*numBulk.
Definition: Bulk.hpp:321
Bulk::depth
vector< OCP_DBL > depth
Depth of center of grid cells: activeGridNum.
Definition: Bulk.hpp:241
Bulk::Ni
vector< OCP_DBL > Ni
Moles of component: numCom*numBulk.
Definition: Bulk.hpp:306
Bulk::Pj
vector< OCP_DBL > Pj
Pressure of phase: numPhase*numBulk.
Definition: Bulk.hpp:311
Bulk::mux
vector< OCP_DBL > mux
d Muj / d xij: numPhase*numCom*numBulk.
Definition: Bulk.hpp:360
Bulk::ifThermal
OCP_BOOL ifThermal
Id OCP_TRUE, ifThermal model will be used.
Definition: Bulk.hpp:220
Bulk::muP
vector< OCP_DBL > muP
d Mu / d P: numPhase*numBulk.
Definition: Bulk.hpp:358
Bulk::dz
vector< OCP_DBL > dz
Size of cell in z-direction: activeGridNum.
Definition: Bulk.hpp:239
Bulk::dKr_dS
vector< OCP_DBL > dKr_dS
d Krj / d Sk: numPhase * numPhase * bulk.
Definition: Bulk.hpp:362
Bulk::flashCal
vector< Mixture * > flashCal
Flash calculation class.
Definition: Bulk.hpp:192
Bulk::phaseExist
vector< OCP_BOOL > phaseExist
Existence of phase: numPhase*numBulk.
Definition: Bulk.hpp:313
Bulk::xiP
vector< OCP_DBL > xiP
d xi / d P: numPhase*numBulk.
Definition: Bulk.hpp:355
Bulk::xi
vector< OCP_DBL > xi
Moles density of phase: numPhase*numBulk.
Definition: Bulk.hpp:319
Bulk::rockKz
vector< OCP_DBL > rockKz
current rock permeability along the z direction.
Definition: Bulk.hpp:249
Bulk::rockKx
vector< OCP_DBL > rockKx
current rock permeability along the x direction.
Definition: Bulk.hpp:247
Bulk::rsTemp
OCP_DBL rsTemp
Reservoir temperature.
Definition: Bulk.hpp:171
Bulk::rockKy
vector< OCP_DBL > rockKy
current rock permeability along the y direction.
Definition: Bulk.hpp:248
Bulk::water
OCP_BOOL water
If OCP_TRUE, water phase could exist.
Definition: Bulk.hpp:224
Bulk::numBulk
OCP_USI numBulk
Number of bulks (active grids).
Definition: Bulk.hpp:153
Bulk::xix
vector< OCP_DBL > xix
d Xi_j / d xij: numPhase*numCom*numBulk.
Definition: Bulk.hpp:357
Bulk::dx
vector< OCP_DBL > dx
Size of cell in x-direction: activeGridNum.
Definition: Bulk.hpp:237
Bulk::dPcj_dS
vector< OCP_DBL > dPcj_dS
d Pcj / d Sk: numPhase * numPhase * bulk.
Definition: Bulk.hpp:361
Bulk::P
vector< OCP_DBL > P
Pressure: numBulk.
Definition: Bulk.hpp:309
Bulk::numCom
USI numCom
Number of component.
Definition: Bulk.hpp:155
Grid
Definition: Grid.hpp:89
Grid::map_Act2All
vector< OCP_USI > map_Act2All
Mapping from active grid to all grid: activeGridNum.
Definition: Grid.hpp:192
Grid::nx
USI nx
Number of cells in x-direction.
Definition: Grid.hpp:145
Grid::polyhedronGrid
vector< OCPpolyhedron > polyhedronGrid
Coordinates of grid points.
Definition: Grid.hpp:219
Grid::map_All2Flu
vector< GB_Pair > map_All2Flu
Mapping from all grid to fluid grid: numGrid.
Definition: Grid.hpp:196
Grid::ny
USI ny
Number of cells in y-direction.
Definition: Grid.hpp:146
Grid::map_All2Act
vector< GB_Pair > map_All2Act
Mapping from grid to active all grid: numGrid.
Definition: Grid.hpp:193
LinearSystem
Linear solvers for discrete systems.
Definition: LinearSystem.hpp:31
LinearSystem::AddRhs
void AddRhs(const OCP_USI &n, const OCP_DBL &v)
Add a value at b[n].
Definition: LinearSystem.hpp:90
LinearSystem::NewOffDiag
void NewOffDiag(const OCP_USI &bId, const OCP_USI &eId, const OCP_DBL &v)
Push back a off-diagonal value.
Definition: LinearSystem.hpp:83
OCPpolyhedron
Record the initial grid information, all of grids are contained.
Definition: Grid.hpp:75
Point3D
A point in 3D.
Definition: CornerGrid.hpp:47
WellParam
TODO: Add Doxygen.
Definition: ParamWell.hpp:59
WellParam::direction
vector< string > direction
Direction of perforations.
Definition: ParamWell.hpp:78
WellParam::diameter
vector< OCP_DBL > diameter
Diameter of perforations.
Definition: ParamWell.hpp:75
WellParam::skinFactor
vector< OCP_DBL > skinFactor
Skin factor.
Definition: ParamWell.hpp:77
WellParam::I_perf
vector< USI > I_perf
I-index of perforation in grid.
Definition: ParamWell.hpp:71
WellParam::WI
vector< OCP_DBL > WI
Transmissibility connection factor.
Definition: ParamWell.hpp:74
WellParam::J_perf
vector< USI > J_perf
J-index of perforation in grid.
Definition: ParamWell.hpp:72
WellParam::K_perf
vector< USI > K_perf
K-index of perforation in grid.
Definition: ParamWell.hpp:73
Well::CalInjRateMaxBHP
OCP_DBL CalInjRateMaxBHP(const Bulk &myBulk)
calculate flow rate of moles of phases for injection well with maxBHP.
Definition: Well.cpp:263
Well::numPerf
USI numPerf
num of perforations belonging to this well.
Definition: Well.hpp:78
Well::WOPT
OCP_DBL WOPT
well total oil production.
Definition: Well.hpp:222
Well::AssembleMatReinjection_IMPEC
void AssembleMatReinjection_IMPEC(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt, const vector< Well > &allWell, const vector< USI > &injId) const
Definition: Well.cpp:1096
Well::numPhase
OCP_DBL numPhase
num of phases
Definition: Well.hpp:83
Well::WWIT
OCP_DBL WWIT
well total water injection.
Definition: Well.hpp:230
Well::CalProddG01
void CalProddG01(const Bulk &myBulk)
Calculate pressure difference between well and perforations for Production.
Definition: Well.cpp:438
Well::prodWeight
vector< OCP_DBL > prodWeight
maybe only a num is needed
Definition: Well.hpp:219
Well::CheckOptMode
void CheckOptMode(const Bulk &myBulk)
Check if well operation mode would be changed.
Definition: Well.cpp:893
Well::CheckP
OCP_INT CheckP(const Bulk &myBulk)
Check if abnormal Pressure occurs.
Definition: Well.cpp:948
Well::prodRate
vector< OCP_DBL > prodRate
flow rate of volume of phase outflowing
Definition: Well.hpp:218
Well::CalInjdG
void CalInjdG(const Bulk &myBulk)
Calculate pressure difference between well and perforations for Injection.
Definition: Well.cpp:367
Well::wOId
USI wOId
Definition: Well.hpp:89
Well::CalProdQj
void CalProdQj(const Bulk &myBulk, const OCP_DBL &dt)
Definition: Well.cpp:340
Well::CalProddG02
void CalProddG02(const Bulk &myBulk)
Calculate pressure difference between well and perforations for Production.
Definition: Well.cpp:557
Well::depth
OCP_DBL depth
reference depth of well.
Definition: Well.hpp:74
Well::CalReInjFluid
void CalReInjFluid(const Bulk &myBulk, vector< OCP_DBL > &myZi)
Calculate the contribution of production well to reinjection defaulted.
Definition: Well.cpp:864
Well::CalFlux
void CalFlux(const Bulk &myBulk, const OCP_BOOL ReCalXi=OCP_FALSE)
Calculate the flux for each perforations.
Definition: Well.cpp:204
Well::InputPerfo
void InputPerfo(const WellParam &well)
Input the param of perforations.
Definition: Well.cpp:15
Well::opt
WellOpt opt
well control parameters, contains current control parameters.
Definition: Well.hpp:75
Well::WGPR
OCP_DBL WGPR
well gas production rate.
Definition: Well.hpp:223
Well::WGIR
OCP_DBL WGIR
well gas injection rate.
Definition: Well.hpp:227
Well::WWPT
OCP_DBL WWPT
well total water production.
Definition: Well.hpp:226
Well::numCom
OCP_DBL numCom
num of components
Definition: Well.hpp:84
Well::WWIR
OCP_DBL WWIR
well water injection rate.
Definition: Well.hpp:229
Well::name
string name
well name
Definition: Well.hpp:70
Well::CalTrans
void CalTrans(const Bulk &myBulk)
Calculate transmissibility for each phase in perforations.
Definition: Well.cpp:164
Well::bhp
OCP_DBL bhp
Well pressure in reference depth.
Definition: Well.hpp:208
Well::IsOpen
OCP_BOOL IsOpen() const
Return the state of the well, Open or Close.
Definition: Well.hpp:149
Well::CalProddG
void CalProddG(const Bulk &myBulk)
Calculate pressure difference between well and perforations for Production.
Definition: Well.cpp:676
Well::Tsurf
OCP_DBL Tsurf
Well surface Temperature, F.
Definition: Well.hpp:87
Well::CaldG
void CaldG(const Bulk &myBulk)
Calculate pressure difference between well and perforations.
Definition: Well.cpp:357
Well::WOPR
OCP_DBL WOPR
well oil production rate.
Definition: Well.hpp:221
Well::WGPT
OCP_DBL WGPT
well total gas production.
Definition: Well.hpp:224
Well::WGIT
OCP_DBL WGIT
well total gas injection.
Definition: Well.hpp:228
Well::CorrectBHP
void CorrectBHP()
Correct BHP if opt mode is BHPMode.
Definition: Well.cpp:882
Well::ldG
vector< OCP_DBL > ldG
Last dG.
Definition: Well.hpp:214
Well::CheckCrossFlow
OCP_INT CheckCrossFlow(const Bulk &myBulk)
Check if cross flow happens.
Definition: Well.cpp:993
Well::CalWI_Peaceman
void CalWI_Peaceman(const Bulk &myBulk)
Calculate Well Index with Peaceman model.
Definition: Well.cpp:93
Well::CalProdRateMinBHP
OCP_DBL CalProdRateMinBHP(const Bulk &myBulk)
calculate flow rate of moles of phases for production well with minBHP.
Definition: Well.cpp:288
Well::Setup
void Setup(const Grid &gd, const Bulk &bk, const vector< SolventINJ > &sols)
Setup the well after Grid and Bulk finish setup.
Definition: Well.cpp:41
Well::CalInjQj
void CalInjQj(const Bulk &myBulk, const OCP_DBL &dt)
Definition: Well.cpp:322
Well::AssembleMatReinjection_FIM
void AssembleMatReinjection_FIM(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt, const vector< Well > &allWell, const vector< USI > &injId) const
Definition: Well.cpp:1145
Well::ShowPerfStatus
void ShowPerfStatus(const Bulk &myBulk) const
Display operation mode of well and state of perforations.
Definition: Well.cpp:1068
Well::optSet
vector< WellOpt > optSet
Definition: Well.hpp:76
Well::flashCal
vector< Mixture * > flashCal
from bulks's flashCal
Definition: Well.hpp:81
Well::Psurf
OCP_DBL Psurf
Well surface Pressure, psia.
Definition: Well.hpp:86
Well::qi_lbmol
vector< OCP_DBL > qi_lbmol
flow rate of moles of component inflowing/outflowing
Definition: Well.hpp:217
Well::CalProdWeight
void CalProdWeight(const Bulk &myBulk) const
Calculate the production weight.
Definition: Well.cpp:829
Well::WWPR
OCP_DBL WWPR
well water production rate.
Definition: Well.hpp:225
Well::dG
vector< OCP_DBL > dG
difference of pressure between well and perforation: numPerf.
Definition: Well.hpp:210
Well::perf
vector< Perforation > perf
information of perforation belonging to this well.
Definition: Well.hpp:79
Generated on Tue Feb 7 2023 01:49:51 for OpenCAEPoro by doxygen 1.9.1