OpenCAEPoro  0.2.0 Sep/22/2022
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 WellOpt::WellOpt(const WellOptParam& Optparam)
16 {
17  if (Optparam.type == "INJ") {
18  type = INJ;
19  } else if (Optparam.type == "PROD") {
20  type = PROD;
21  } else {
22  OCP_ABORT("Wrong well type!");
23  }
24 
25  if (type == INJ) {
26  fluidType = Optparam.fluidType;
27  if (fluidType == "WAT" || fluidType == "WATER") {
28  fluidType = "WAT";
29  }
30  }
31 
32  if (Optparam.state == "OPEN") {
33  state = OPEN;
34  } else if (Optparam.state == "CLOSE") {
35  state = CLOSE;
36  } else {
37  OCP_ABORT("Wrong state type!");
38  }
39 
40  if (Optparam.optMode == "RATE") {
41  optMode = RATE_MODE;
42  } else if (Optparam.optMode == "ORAT") {
43  optMode = ORATE_MODE;
44  } else if (Optparam.optMode == "GRAT") {
45  optMode = GRATE_MODE;
46  } else if (Optparam.optMode == "WRAT") {
47  optMode = WRATE_MODE;
48  } else if (Optparam.optMode == "LRAT") {
49  optMode = LRATE_MODE;
50  } else if (Optparam.optMode == "BHP") {
51  optMode = BHP_MODE;
52  } else {
53  OCP_ABORT("Wrong well option mode!");
54  }
55 
56  initOptMode = optMode;
57  maxRate = Optparam.maxRate;
58  maxBHP = Optparam.maxBHP;
59  minBHP = Optparam.minBHP;
60 }
61 
62 
63 bool WellOpt::operator !=(const WellOpt& Opt) const
64 {
65  if (this->type != Opt.type) return true;
66  if (this->state != Opt.state) return true;
67  if (this->optMode != Opt.optMode) return true;
68  if (this->initOptMode != Opt.initOptMode) return true;
69  if (this->maxRate != Opt.maxRate) return true;
70  if (this->maxBHP != Opt.maxBHP) return true;
71  if (this->minBHP != Opt.minBHP) return true;
72  for (USI i = 0; i < zi.size(); i++) {
73  if (fabs(zi[i] - Opt.zi[i]) > TINY) return true;
74  }
75  return false;
76 }
77 
78 
79 void Well::InputPerfo(const WellParam& well)
80 {
81  OCP_FUNCNAME;
82 
83  numPerf = well.I_perf.size();
84  perf.resize(numPerf);
85  for (USI p = 0; p < numPerf; p++) {
86  perf[p].I = well.I_perf[p] - 1;
87  perf[p].J = well.J_perf[p] - 1;
88  perf[p].K = well.K_perf[p] - 1;
89  perf[p].WI = well.WI[p];
90  perf[p].radius = well.diameter[p] / 2.0;
91  perf[p].kh = well.kh[p];
92  perf[p].skinFactor = well.skinFactor[p];
93  if (well.direction[p] == "X" || well.direction[p] == "x") {
94  perf[p].direction = X_DIRECTION;
95  } else if (well.direction[p] == "Y" || well.direction[p] == "y") {
96  perf[p].direction = Y_DIRECTION;
97  } else if (well.direction[p] == "Z" || well.direction[p] == "z") {
98  perf[p].direction = Z_DIRECTION;
99  } else {
100  OCP_ABORT("Wrong direction of perforations!");
101  }
102  }
103 }
104 
105 void Well::Setup(const Grid& myGrid, const Bulk& myBulk, const vector<SolventINJ>& sols)
106 {
107  OCP_FUNCNAME;
108  qi_lbmol.resize(myBulk.numCom);
109  prodWeight.resize(myBulk.numCom);
110  factor.resize(3); // oil, gas, liquid
111  Mtype = myBulk.flashCal[0]->GetType();
112  // zi
113  if (myBulk.blackOil) {
114  for (auto& opt : optSet) {
115 
116  if (!opt.state) continue;
117 
118  opt.zi.resize(myBulk.numCom, 0);
119  if (opt.type == INJ) {
120  // INJ
121  switch (myBulk.PVTmode) {
122  case PHASE_W:
123  case PHASE_OW:
124  opt.zi.back() = 1;
125  break;
126  case PHASE_ODGW:
127  case PHASE_DOGW:
128  if (opt.fluidType == "GAS")
129  opt.zi[1] = 1;
130  else
131  opt.zi[2] = 1;
132  break;
133  default:
134  OCP_ABORT("Wrong blackoil type!");
135  }
136  } else {
137  // PROD
138  switch (myBulk.PVTmode) {
139  case PHASE_W:
140  opt.zi.back() = 1;
141  break;
142  case PHASE_OW:
143  if (opt.optMode == ORATE_MODE)
144  opt.zi[0] = 1;
145  else
146  opt.zi[1] = 1;
147  break;
148  case PHASE_DOGW:
149  case PHASE_ODGW:
150  if (opt.optMode == ORATE_MODE)
151  opt.zi[0] = 1;
152  else if (opt.optMode == GRATE_MODE)
153  opt.zi[1] = 1;
154  else if (opt.optMode == WRATE_MODE)
155  opt.zi[2] = 1;
156  else if (opt.optMode == LRATE_MODE)
157  opt.zi[2] = opt.zi[0] = 1;
158  break;
159  default:
160  OCP_ABORT("Wrong blackoil type!");
161  }
162  }
163  }
164  } else if (myBulk.comps) {
165 
166  USI len = sols.size();
167 
168  for (auto& opt : optSet) {
169 
170  if (!opt.state) continue;
171 
172  if (opt.type == INJ) {
173  // INJ Well
174  if (opt.fluidType == "WAT") {
175  opt.zi.resize(myBulk.numCom, 0);
176  opt.zi.back() = 1;
177  } else {
178  for (USI i = 0; i < len; i++) {
179  if (opt.fluidType == sols[i].name) {
180  opt.zi = sols[i].data;
181  opt.zi.resize(myBulk.numCom);
182  // Convert volume units Mscf/stb to molar units lbmoles for
183  // injfluid Use flash in Bulk in surface condition
184  opt.xiINJ = myBulk.flashCal[0]->XiPhase(
185  PRESSURE_STD, TEMPERATURE_STD, &opt.zi[0]);
186  opt.maxRate *=
187  (opt.xiINJ *
188  1000); // lbmol / ft3 -> lbmol / Mscf for gas
189  break;
190  }
191  if (i == len - 1) {
192  OCP_ABORT("Wrong FluidType!");
193  }
194  }
195  }
196  }
197  // else {
198  // // PROD Well use EoS
199  // }
200  }
201  } else {
202  OCP_ABORT("Wrong mixture type!");
203  }
204 
205  // perf
206  USI pp = 0;
207  for (USI p = 0; p < numPerf; p++) {
208  OCP_USI Idg =
209  perf[p].K * myGrid.nx * myGrid.ny + perf[p].J * myGrid.nx + perf[p].I;
210  if (myGrid.activeMap_G2B[Idg].IsAct()) {
211 
212  perf[pp] = perf[p];
213  perf[pp].state = OPEN;
214  perf[pp].location = myGrid.activeMap_G2B[Idg].GetId();
215  perf[pp].depth = myBulk.depth[perf[pp].location];
216  perf[pp].multiplier = 1;
217  perf[pp].qi_lbmol.resize(myBulk.numCom);
218  perf[pp].transj.resize(myBulk.numPhase);
219  perf[pp].qj_ft3.resize(myBulk.numPhase);
220  pp++;
221  } else {
222  OCP_WARNING("Perforation is in inactive bulk!");
223  }
224  }
225  numPerf = pp;
226  perf.resize(numPerf);
227  // dG
228  dG.resize(numPerf, 0);
229  ldG = dG;
230 
231  if (depth < 0) depth = perf[0].depth;
232 
233  CalWI_Peaceman_Vertical(myBulk);
234  // test
235  // ShowPerfStatus(myBulk);
236 }
237 
238 void Well::InitBHP(const Bulk& myBulk) { BHP = myBulk.P[perf[0].location]; }
239 
240 void Well::CalWI_Peaceman_Vertical(const Bulk& myBulk)
241 {
242  OCP_FUNCNAME;
243 
244  // this fomular needs to be carefully checked !
245  // especially the dz
246 
247  for (USI p = 0; p < numPerf; p++) {
248  if (perf[p].WI > 0) {
249  break;
250  } else {
251  OCP_USI Idb = perf[p].location;
252  OCP_DBL dx = myBulk.dx[Idb];
253  OCP_DBL dy = myBulk.dy[Idb];
254  OCP_DBL dz = myBulk.dz[Idb] * myBulk.ntg[Idb];
255  OCP_DBL ro = 0;
256  switch (perf[p].direction) {
257  case X_DIRECTION: {
258  OCP_DBL kykz = myBulk.rockKy[Idb] * myBulk.rockKz[Idb];
259  OCP_DBL ky_kz = myBulk.rockKy[Idb] / myBulk.rockKz[Idb];
260  assert(kykz > 0);
261  ro =
262  0.28 *
263  pow((dy * dy * pow(1 / ky_kz, 0.5) + dz * dz * pow(ky_kz, 0.5)),
264  0.5);
265  ro /= (pow(ky_kz, 0.25) + pow(1 / ky_kz, 0.25));
266 
267  if (perf[p].kh < 0) {
268  perf[p].kh = (dx * pow(kykz, 0.5));
269  }
270  break;
271  }
272  case Y_DIRECTION: {
273  OCP_DBL kzkx = myBulk.rockKz[Idb] * myBulk.rockKx[Idb];
274  OCP_DBL kz_kx = myBulk.rockKz[Idb] / myBulk.rockKx[Idb];
275  assert(kzkx > 0);
276  ro =
277  0.28 *
278  pow((dz * dz * pow(1 / kz_kx, 0.5) + dx * dx * pow(kz_kx, 0.5)),
279  0.5);
280  ro /= (pow(kz_kx, 0.25) + pow(1 / kz_kx, 0.25));
281 
282  if (perf[p].kh < 0) {
283  perf[p].kh = (dy * pow(kzkx, 0.5));
284  }
285  break;
286  }
287  case Z_DIRECTION: {
288  OCP_DBL kxky = myBulk.rockKx[Idb] * myBulk.rockKy[Idb];
289  OCP_DBL kx_ky = myBulk.rockKx[Idb] / myBulk.rockKy[Idb];
290  assert(kxky > 0);
291  ro =
292  0.28 *
293  pow((dx * dx * pow(1 / kx_ky, 0.5) + dy * dy * pow(kx_ky, 0.5)),
294  0.5);
295  ro /= (pow(kx_ky, 0.25) + pow(1 / kx_ky, 0.25));
296 
297  if (perf[p].kh < 0) {
298  perf[p].kh = (dz * pow(kxky, 0.5));
299  }
300  break;
301  }
302  default:
303  OCP_ABORT("Wrong direction of perforations!");
304  }
305  perf[p].WI = CONV2 * (2 * PI) * perf[p].kh /
306  (log(ro / perf[p].radius) + perf[p].skinFactor);
307  }
308  }
309 }
310 
311 void Well::CalTrans(const Bulk& myBulk)
312 {
313  OCP_FUNCNAME;
314 
315  const USI np = myBulk.numPhase;
316 
317  if (opt.type == INJ) {
318  for (USI p = 0; p < numPerf; p++) {
319  perf[p].transINJ = 0;
320  OCP_USI k = perf[p].location;
321  OCP_DBL temp = CONV1 * perf[p].WI * perf[p].multiplier;
322 
323  // single phase
324  for (USI j = 0; j < np; j++) {
325  perf[p].transj[j] = 0;
326  OCP_USI id = k * np + j;
327  if (myBulk.phaseExist[id]) {
328  perf[p].transj[j] = temp * myBulk.kr[id] / myBulk.mu[id];
329  perf[p].transINJ += perf[p].transj[j];
330  }
331  }
332  }
333  } else {
334  for (USI p = 0; p < numPerf; p++) {
335  OCP_USI k = perf[p].location;
336  OCP_DBL temp = CONV1 * perf[p].WI * perf[p].multiplier;
337 
338  // multi phase
339  for (USI j = 0; j < np; j++) {
340  perf[p].transj[j] = 0;
341  OCP_USI id = k * np + j;
342  if (myBulk.phaseExist[id]) {
343  perf[p].transj[j] = temp * myBulk.kr[id] / myBulk.mu[id];
344  }
345  }
346  }
347  }
348  // check Trans
349  // if (opt.type == PROD) {
350 
351  // USI count = 0;
352  // for (USI p = 0; p < numPerf; p++) {
353  // if (perf[p].transj[0] != 0) {
354  // count++;
355  // break;
356  // }
357  // }
358  // if (count == 0) {
359  // cout << name << endl;
360  // for (USI p = 0; p < numPerf; p++) {
361  // OCP_USI k = perf[p].location;
362  // cout << "perf " << p << " " << perf[p].multiplier << " " <<
363  // myBulk.S[k * np] << " " << myBulk.kr[k * np] << " " <<
364  // myBulk.S[k * np + 1] << " " << myBulk.kr[k * np + 1] << " " <<
365  // myBulk.S[k * np + 2] << " " << myBulk.kr[k * np + 2] << " " <<
366  // endl;
367  // }
368 
369  // }
370  //}
371 }
372 
373 void Well::CalFlux(const Bulk& myBulk, const bool flag)
374 {
375  OCP_FUNCNAME;
376 
377  // cout << name << endl;
378 
379  const USI np = myBulk.numPhase;
380  const USI nc = myBulk.numCom;
381  fill(qi_lbmol.begin(), qi_lbmol.end(), 0.0);
382 
383  if (opt.type == INJ) {
384 
385  for (USI p = 0; p < numPerf; p++) {
386  perf[p].P = BHP + dG[p];
387  OCP_USI k = perf[p].location;
388  OCP_DBL dP = perf[p].P - myBulk.P[k];
389  dP *= -1.0;
390 
391  perf[p].qt_ft3 = perf[p].transINJ * dP;
392 
393  if (flag) {
394  USI pvtnum = myBulk.PVTNUM[k];
395  perf[p].xi =
396  myBulk.flashCal[pvtnum]->XiPhase(myBulk.P[k], myBulk.T, &opt.zi[0]);
397  }
398  // cout << "xi: " << setprecision(6) << perf[p].xi << endl;
399  OCP_DBL xi = perf[p].xi;
400  for (USI i = 0; i < nc; i++) {
401  perf[p].qi_lbmol[i] = perf[p].qt_ft3 * xi * opt.zi[i];
402  qi_lbmol[i] += perf[p].qi_lbmol[i];
403  }
404  }
405  } else {
406 
407  for (USI p = 0; p < numPerf; p++) {
408  perf[p].P = BHP + dG[p];
409  OCP_USI k = perf[p].location;
410  perf[p].qt_ft3 = 0;
411  fill(perf[p].qi_lbmol.begin(), perf[p].qi_lbmol.end(), 0.0);
412  fill(perf[p].qj_ft3.begin(), perf[p].qj_ft3.end(), 0.0);
413 
414  for (USI j = 0; j < np; j++) {
415  OCP_USI id = k * np + j;
416  if (myBulk.phaseExist[id]) {
417  OCP_DBL dP = myBulk.Pj[id] - perf[p].P;
418 
419  perf[p].qj_ft3[j] = perf[p].transj[j] * dP;
420  perf[p].qt_ft3 += perf[p].qj_ft3[j];
421  // cout << p << " p[" << j << "] = " << myBulk.Pj[id] << endl;
422  // cout << p << " perf = " << perf[p].p << endl;
423 
424  OCP_DBL xi = myBulk.xi[id];
425  OCP_DBL xij;
426  for (USI i = 0; i < nc; i++) {
427  xij = myBulk.xij[id * nc + i];
428  perf[p].qi_lbmol[i] += perf[p].qj_ft3[j] * xi * xij;
429  }
430  }
431  }
432 
433  // check if perf[p].qi_lbmol is zero vector
434  // bool flag = false;
435  // for (USI i = 0; i < nc; i++) {
436  // if (perf[p].qi_lbmol[i] != 0) {
437  // flag = true;
438  // break;
439  // }
440  //}
441  // if (!flag) {
442  // OCP_ABORT("Qi is all zero!");
443  //}
444 
445  for (USI i = 0; i < nc; i++) qi_lbmol[i] += perf[p].qi_lbmol[i];
446  }
447  }
448  // test
449  // cout << name << "----" << endl;
450  // vector<OCP_DBL> tmpNiP(qi_lbmol);
451  // OCP_DBL qt = Dnorm1(nc, &qi_lbmol[0]);
452  // Dscalar(myBulk.numCom, 1 / qt * 100, &tmpNiP[0]);
453  // cout << qt << " ";
454  // PrintDX(myBulk.numCom, &tmpNiP[0]);
455 }
456 
460 OCP_DBL Well::CalInjRate(const Bulk& myBulk, const bool& maxBHP)
461 {
462  OCP_FUNCNAME;
463 
464  OCP_DBL qj = 0;
465  OCP_DBL Pwell = maxBHP ? opt.maxBHP : BHP;
466 
467  for (USI p = 0; p < numPerf; p++) {
468 
469  OCP_DBL Pperf = Pwell + dG[p];
470  OCP_USI k = perf[p].location;
471 
472  OCP_DBL dP = Pperf - myBulk.P[k];
473  qj += perf[p].transINJ * perf[p].xi * dP;
474  }
475  return qj;
476 }
477 
481 OCP_DBL Well::CalProdRate(const Bulk& myBulk, const bool& minBHP)
482 {
483  OCP_FUNCNAME;
484 
485  const USI np = myBulk.numPhase;
486  const USI nc = myBulk.numCom;
487  OCP_DBL qj = 0;
488  OCP_DBL Pwell = minBHP ? opt.minBHP : BHP;
489 
490  for (USI p = 0; p < numPerf; p++) {
491 
492  OCP_DBL Pperf = Pwell + dG[p];
493  OCP_USI k = perf[p].location;
494 
495  for (USI j = 0; j < np; j++) {
496  OCP_USI id = k * np + j;
497  if (myBulk.phaseExist[id]) {
498  OCP_DBL temp = 0;
499  for (USI i = 0; i < nc; i++) {
500  temp += prodWeight[i] * myBulk.xij[id * nc + i];
501  }
502  OCP_DBL xi = myBulk.xi[id];
503  OCP_DBL dP = myBulk.Pj[id] - Pperf;
504  qj += perf[p].transj[j] * xi * dP * temp;
505  }
506  }
507  }
508  return qj;
509 }
510 
511 void Well::CalInjQi(const Bulk& myBulk, const OCP_DBL& dt)
512 {
513  OCP_FUNCNAME;
514 
515  const USI nc = myBulk.numCom;
516  OCP_DBL qj = 0;
517 
518  for (USI i = 0; i < nc; i++) {
519  qj += qi_lbmol[i];
520  }
521  if (opt.fluidType == "WAT") {
522  WWIR = -qj;
523  WWIT += WWIR * dt;
524  } else {
525  if (Mtype == BLKOIL) {
526  WGIR = -qj;
527  } else {
528  // EoS PVTW
529  WGIR = -qj / (opt.xiINJ * 1000); // lbmol -> Mscf
530  }
531  WGIT += WGIR * dt;
532  }
533 }
534 
535 void Well::CalProdQj(const Bulk& myBulk, const OCP_DBL& dt)
536 {
537  if (Mtype == BLKOIL) {
538  CalProdQiBO(myBulk);
539  } else {
540  // EoS PVTW
541  CalProdQjCOMP(myBulk);
542  }
543  WOPT += WOPR * dt;
544  WGPT += WGPR * dt;
545  WWPT += WWPR * dt;
546 }
547 
548 void Well::CalProdQjCOMP(const Bulk& myBulk)
549 {
550  USI nc = myBulk.numCom;
551  OCP_DBL qt = Dnorm1(nc, &qi_lbmol[0]);
552  WOPR = qt * factor[0];
553  WGPR = qt * factor[1];
554  WWPR = qi_lbmol[nc - 1];
555 }
556 
557 void Well::CalProdQiBO(const Bulk& myBulk)
558 {
559  OCP_FUNCNAME;
560 
561  const USI nc = myBulk.numCom;
562 
563  for (USI i = 0; i < nc; i++) {
564  if (myBulk.index2Phase[i] == OIL) {
565  WOPR = qi_lbmol[i];
566  } else if (myBulk.index2Phase[i] == GAS) {
567  WGPR = qi_lbmol[i];
568  } else if (myBulk.index2Phase[i] == WATER) {
569  WWPR = qi_lbmol[i];
570  }
571  }
572 }
573 
577 void Well::CaldG(const Bulk& myBulk)
578 {
579  OCP_FUNCNAME;
580 
581  if (opt.type == INJ)
582  CalInjdG(myBulk);
583  else
584  CalProddG(myBulk);
585 
586  // if (name == "PROD17") {
587  // for (USI p = 0; p < numPerf; p++) {
588  // cout << perf[p].state << " ";
589  // cout << setprecision(2);
590  // cout << dG[p] << " ";
591  // OCP_USI n = perf[p].location * myBulk.numPhase;
592  // cout << setprecision(6);
593  // cout << myBulk.kr[n + 0] << " ";
594  // cout << myBulk.kr[n + 1] << " ";
595  // cout << myBulk.kr[n + 2] << " ";
596  // cout << myBulk.S[n + 0] << " ";
597  // cout << myBulk.S[n + 1] << " ";
598  // cout << myBulk.S[n + 2] << " ";
599  // cout << setprecision(2);
600  // for (USI i = 0; i < myBulk.numCom; i++) {
601  // cout << perf[p].qi_lbmol[i] << " ";
602  // }
603  // cout << endl;
604  // }
605  //}
606 }
607 
608 void Well::CalInjdG(const Bulk& myBulk)
609 {
610  OCP_FUNCNAME;
611 
612  const OCP_DBL maxlen = 10;
613  USI seg_num = 0;
614  OCP_DBL seg_len = 0;
615  vector<OCP_DBL> dGperf(numPerf, 0);
616 
617  if (depth <= perf.front().depth) {
618  // Well is higher
619  for (OCP_INT p = numPerf - 1; p >= 0; p--) {
620  if (p == 0) {
621  seg_num = ceil(fabs((perf[0].depth - depth) / maxlen));
622  if (seg_num == 0) continue;
623  seg_len = (perf[0].depth - depth) / seg_num;
624  } else {
625  seg_num = ceil(fabs((perf[p].depth - perf[p - 1].depth) / maxlen));
626  if (seg_num == 0) continue;
627  seg_len = (perf[p].depth - perf[p - 1].depth) / seg_num;
628  }
629  OCP_USI n = perf[p].location;
630  perf[p].P = BHP + dG[p];
631  OCP_DBL Pperf = perf[p].P;
632  OCP_DBL Ptmp = Pperf;
633 
634  USI pvtnum = myBulk.PVTNUM[n];
635  for (USI i = 0; i < seg_num; i++) {
636  Ptmp -=
637  myBulk.flashCal[pvtnum]->RhoPhase(Ptmp, myBulk.T, opt.zi.data()) *
638  GRAVITY_FACTOR * seg_len;
639  }
640  dGperf[p] = Pperf - Ptmp;
641  }
642  dG[0] = dGperf[0];
643  for (USI p = 1; p < numPerf; p++) {
644  dG[p] = dG[p - 1] + dGperf[p];
645  }
646  } else if (depth >= perf[numPerf - 1].depth) {
647  // Well is lower
648  for (USI p = 0; p < numPerf; p++) {
649  if (p == numPerf - 1) {
650  seg_num = ceil(fabs((depth - perf[numPerf - 1].depth) / maxlen));
651  if (seg_num == 0) continue;
652  seg_len = (depth - perf[numPerf - 1].depth) / seg_num;
653  } else {
654  seg_num = ceil(fabs((perf[p + 1].depth - perf[p].depth) / maxlen));
655  if (seg_num == 0) continue;
656  seg_len = (perf[p + 1].depth - perf[p].depth) / seg_num;
657  }
658  OCP_USI n = perf[p].location;
659  perf[p].P = BHP + dG[p];
660  OCP_DBL Pperf = perf[p].P;
661  OCP_DBL Ptmp = Pperf;
662 
663  USI pvtnum = myBulk.PVTNUM[n];
664  for (USI i = 0; i < seg_num; i++) {
665  Ptmp +=
666  myBulk.flashCal[pvtnum]->RhoPhase(Ptmp, myBulk.T, opt.zi.data()) *
667  GRAVITY_FACTOR * seg_len;
668  }
669  dGperf[p] = Ptmp - Pperf;
670  }
671  dG[numPerf - 1] = dGperf[numPerf - 1];
672  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
673  dG[p] = dG[p + 1] + dGperf[p];
674  }
675  }
676 }
677 
678 void Well::CalProddG(const Bulk& myBulk)
679 {
680  OCP_FUNCNAME;
681 
682  const USI np = myBulk.numPhase;
683  const USI nc = myBulk.numCom;
684  const OCP_DBL maxlen = 5;
685  USI seg_num = 0;
686  OCP_DBL seg_len = 0;
687  vector<OCP_DBL> tmpNi(nc, 0);
688  vector<OCP_DBL> dGperf(numPerf, 0);
689  OCP_DBL qtacc = 0;
690  OCP_DBL rhoacc = 0;
691  OCP_DBL rhotmp = 0;
692 
693  if (depth <= perf.front().depth) {
694  // Well is higher
695 
696  // check qi_lbmol ---- test
697  if (perf[numPerf - 1].state == CLOSE) {
698  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
699  if (perf[p].state == OPEN) {
700  for (USI i = 0; i < nc; i++) {
701  perf[numPerf - 1].qi_lbmol[i] = perf[p].qi_lbmol[i];
702  }
703  break;
704  }
705  }
706  }
707 
708  for (OCP_INT p = numPerf - 1; p >= 0; p--) {
709 
710  if (p == 0) {
711  seg_num = ceil(fabs((perf[0].depth - depth) / maxlen));
712  if (seg_num == 0) continue;
713  seg_len = (perf[0].depth - depth) / seg_num;
714  } else {
715  seg_num = ceil(fabs((perf[p].depth - perf[p - 1].depth) / maxlen));
716  if (seg_num == 0) continue;
717  seg_len = (perf[p].depth - perf[p - 1].depth) / seg_num;
718  }
719 
720  OCP_USI n = perf[p].location;
721  perf[p].P = BHP + dG[p];
722  OCP_DBL Pperf = perf[p].P;
723  OCP_DBL Ptmp = Pperf;
724 
725  USI pvtnum = myBulk.PVTNUM[n];
726  // tmpNi.assign(nc, 0);
727  // for (OCP_INT p1 = numPerf - 1; p1 >= p; p1--) {
728  // for (USI i = 0; i < nc; i++) {
729  // tmpNi[i] += perf[p1].qi_lbmol[i];
730  // }
731  //}
732 
733  for (USI i = 0; i < nc; i++) {
734  tmpNi[i] += perf[p].qi_lbmol[i];
735  }
736 
737  // check tmpNi
738  for (auto& v : tmpNi) {
739  v = fabs(v);
740  }
741 
742  for (USI k = 0; k < seg_num; k++) {
743  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T, tmpNi.data(), 0, 0, 0);
744  for (USI j = 0; j < myBulk.numPhase; j++) {
745  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
746  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
747  qtacc += myBulk.flashCal[pvtnum]->v[j] / seg_num;
748  rhoacc += myBulk.flashCal[pvtnum]->v[j] * rhotmp *
749  GRAVITY_FACTOR / seg_num;
750 #ifdef _DEBUG
751  if (rhotmp <= 0 || !isfinite(rhotmp)) {
752  OCP_ABORT("Wrong rho " + to_string(rhotmp));
753  }
754 #endif // _DEBUG
755  }
756  }
757  Ptmp -= rhoacc / qtacc * seg_len;
758  }
759  dGperf[p] = Pperf - Ptmp;
760  }
761  dG[0] = dGperf[0];
762  for (USI p = 1; p < numPerf; p++) {
763  dG[p] = dG[p - 1] + dGperf[p];
764  }
765  } else if (depth >= perf.back().depth) {
766  // Well is lower
767 
768  // check qi_lbmol ---- test
769  if (perf[0].state == CLOSE) {
770  for (USI p = 1; p <= numPerf; p++) {
771  if (perf[p].state == OPEN) {
772  for (USI i = 0; i < nc; i++) {
773  perf[numPerf - 1].qi_lbmol[i] = perf[p].qi_lbmol[i];
774  }
775  break;
776  }
777  }
778  }
779 
780  for (USI p = 0; p < numPerf; p++) {
781  if (p == numPerf - 1) {
782  seg_num = ceil(fabs((depth - perf[numPerf - 1].depth) / maxlen));
783  if (seg_num == 0) continue;
784  seg_len = (depth - perf[numPerf - 1].depth) / seg_num;
785  } else {
786  seg_num = ceil(fabs((perf[p + 1].depth - perf[p].depth) / maxlen));
787  if (seg_num == 0) continue;
788  seg_len = (perf[p + 1].depth - perf[p].depth) / seg_num;
789  }
790 
791  OCP_USI n = perf[p].location;
792  perf[p].P = BHP + dG[p];
793  OCP_DBL Pperf = perf[p].P;
794  OCP_DBL Ptmp = Pperf;
795 
796  USI pvtnum = myBulk.PVTNUM[n];
797  fill(tmpNi.begin(), tmpNi.end(), 0.0);
798  for (OCP_INT p1 = numPerf - 1; p1 - p >= 0; p1--) {
799  for (USI i = 0; i < nc; i++) {
800  tmpNi[i] += perf[p1].qi_lbmol[i];
801  }
802  }
803 
804  // check tmpNi
805  for (auto& v : tmpNi) {
806  v = fabs(v);
807  }
808 
809  for (USI k = 0; k < seg_num; k++) {
810  myBulk.flashCal[pvtnum]->Flash(Ptmp, myBulk.T, tmpNi.data(), 0, 0, 0);
811  for (USI j = 0; j < np; j++) {
812  if (myBulk.flashCal[pvtnum]->phaseExist[j]) {
813  rhotmp = myBulk.flashCal[pvtnum]->rho[j];
814  qtacc += myBulk.flashCal[pvtnum]->v[j] / seg_num;
815  rhoacc += myBulk.flashCal[pvtnum]->v[j] * rhotmp *
816  GRAVITY_FACTOR / seg_num;
817  }
818  }
819  Ptmp += rhoacc / qtacc * seg_len;
820  }
821  dGperf[p] = Ptmp - Pperf;
822  }
823  dG[numPerf - 1] = dGperf[numPerf - 1];
824  for (OCP_INT p = numPerf - 2; p >= 0; p--) {
825  dG[p] = dG[p + 1] + dGperf[p];
826  }
827 
828  } else {
829  OCP_ABORT("Wrong well position!");
830  }
831 }
832 
833 void Well::CalProdWeight(const Bulk& myBulk) const
834 {
835  if (opt.type == PROD) {
836  if (Mtype == BLKOIL || opt.optMode == WRATE_MODE) {
837  prodWeight = opt.zi;
838  } else if (Mtype == EOS_PVTW) {
839  // use qi_lbmol, so it must be calculated before (Well::CalFlux())
840  // when reset the step, qi_lbmol may be calculated
841  USI pid = 0; // phase Id for oil
842  if (opt.optMode == GRATE_MODE) {
843  pid = 1; // phase Id for gas
844  } else if (opt.optMode == LRATE_MODE) {
845  pid = 2; // phase Id for liquid
846  }
847  // in some cases, qi_lbmol may be zero, so use other methods
848  OCP_DBL qt = 0;
849  for (USI i = 0; i < myBulk.numCom; i++) {
850  qt += qi_lbmol[i];
851  }
852  if (fabs(qt) > TINY && qt > 0) {
853 
854  /*cout << name << endl;
855  vector<OCP_DBL> tmpNiP(qi_lbmol);
856  Dscalar(myBulk.numCom, 1 / qt * 100, &tmpNiP[0]);
857  PrintDX(myBulk.numCom, &tmpNiP[0]);*/
858 
859  myBulk.flashCal[0]->Flash(PRESSURE_STD, TEMPERATURE_STD, &qi_lbmol[0],
860  0, 0, 0);
861  } else {
862  USI np = myBulk.numPhase;
863  USI nc = myBulk.numCom;
864  vector<OCP_DBL> tmpNi(nc, 0);
865  for (USI p = 0; p < numPerf; p++) {
866  OCP_USI n = perf[p].location;
867 
868  for (USI j = 0; j < np; j++) {
869  OCP_USI id = n * np + j;
870  if (!myBulk.phaseExist[id]) continue;
871  for (USI k = 0; k < nc; k++) {
872  tmpNi[k] += perf[p].transj[j] * myBulk.xi[id] *
873  myBulk.xij[id * nc + k];
874  }
875  }
876  }
877  qt = Dnorm1(nc, &tmpNi[0]);
878 
879  /*cout << name << endl;
880  vector<OCP_DBL> tmpNiP(tmpNi);
881  Dscalar(nc, 1 / qt, &tmpNiP[0]);
882  PrintDX(nc, &tmpNiP[0]);*/
883 
884  myBulk.flashCal[0]->Flash(PRESSURE_STD, TEMPERATURE_STD, &tmpNi[0], 0,
885  0, 0);
886  }
887  // test
888  // USI nc = myBulk.numCom;
889  // USI np = myBulk.numPhase;
890  // cout << "OIL : ";
891  // for (USI i = 0; i < nc; i++) {
892  // cout << myBulk.flashCal[0]->xij[i] * 100 << " ";
893  //}
894  // cout << endl << "GAS : ";
895  // for (USI i = 0; i < nc; i++) {
896  // cout << myBulk.flashCal[0]->xij[nc + i] * 100 << " ";
897  //}
898  // cout << endl;
899 
900  fill(factor.begin(), factor.end(), 0);
901  if (myBulk.flashCal[0]->phaseExist[0]) { // improve
902  OCP_DBL nuo = myBulk.flashCal[0]->xi[0] * myBulk.flashCal[0]->v[0] /
903  qt; // mole fraction of oil phase
904  factor[0] = nuo / (myBulk.flashCal[0]->xi[0] *
905  CONV1); // lbmol / ft3 -> lbmol / stb for oil
906  }
907  if (myBulk.flashCal[0]->phaseExist[1]) {
908  OCP_DBL nug = myBulk.flashCal[0]->xi[1] * myBulk.flashCal[0]->v[1] /
909  qt; // mole fraction of gas phase
910  factor[1] = nug / (myBulk.flashCal[0]->xi[1] *
911  1000); // lbmol / ft3 -> lbmol / Mscf for gas
912  }
913  factor[2] = factor[0];
914  if (myBulk.flashCal[0]->phaseExist[2]) {
915  OCP_DBL nuw = myBulk.flashCal[0]->xi[2] * myBulk.flashCal[0]->v[2] /
916  qt; // mole fraction of water phase
917  factor[2] += nuw;
918  }
919 
920  if (factor[pid] < 1E-12 || !isfinite(factor[pid])) {
921  OCP_ABORT("Wrong Condition!");
922  }
923  fill(prodWeight.begin(), prodWeight.end(), factor[pid]);
924 
925  // cout << "Factor " << factor[0] << " " << factor[1] << endl;
926  }
927  }
928 }
929 
930 void Well::CalReInjFluid(const Bulk& myBulk, vector<OCP_DBL>& myZi)
931 {
932  CalTrans(myBulk);
933 
934  USI np = myBulk.numPhase;
935  USI nc = myBulk.numCom;
936  for (USI p = 0; p < numPerf; p++) {
937  OCP_USI n = perf[p].location;
938 
939  for (USI j = 0; j < np; j++) {
940  OCP_USI id = n * np + j;
941  if (!myBulk.phaseExist[id]) continue;
942  for (USI k = 0; k < nc; k++) {
943  myZi[k] += perf[p].transj[j] * myBulk.xi[id] * myBulk.xij[id * nc + k];
944  }
945  }
946  }
947 }
948 
949 void Well::SetBHP()
950 {
951  if (opt.type == PROD && opt.optMode == BHP_MODE) {
952  BHP = opt.minBHP;
953  } else if (opt.type == INJ && opt.optMode == BHP_MODE) {
954  BHP = opt.maxBHP;
955  }
956 }
957 
959 void Well::SmoothdG()
960 {
961  OCP_FUNCNAME;
962 
963  for (USI p = 0; p < numPerf; p++) {
964  dG[p] = (ldG[p] + dG[p]) / 2; // seems better
965  // dG[p] = ldG[p] + 0.618 * (dG[p] - ldG[p]);
966  }
967 }
968 
971 void Well::CheckOptMode(const Bulk& myBulk)
972 {
973  OCP_FUNCNAME;
974  if (opt.initOptMode == BHP_MODE) {
975  if (opt.type == INJ) {
976  OCP_DBL q = CalInjRate(myBulk, true);
977  // for INJ well, maxRate has been switch to lbmols
978  OCP_DBL tarRate = opt.maxRate;
979  if (opt.reInj) {
980  if (opt.injPhase == GAS)
981  tarRate = WGIR;
982  else if (opt.injPhase == WATER)
983  tarRate = WWIR;
984  }
985  if (q > tarRate) {
986  opt.optMode = RATE_MODE;
987  }
988  else {
989  opt.optMode = BHP_MODE;
990  BHP = opt.maxBHP;
991  }
992  }
993  else {
994  opt.optMode = BHP_MODE;
995  if (opt.type == INJ)
996  BHP = opt.maxBHP;
997  else
998  BHP = opt.minBHP;
999  }
1000  //else {
1001  // OCP_DBL q = CalProdRate(myBulk, false);
1002  // // cout << q << endl;
1003  // if (q > opt.maxRate) {
1004  // opt.optMode = opt.initOptMode;
1005  // }
1006  // else {
1007  // opt.optMode = BHP_MODE;
1008  // BHP = opt.minBHP;
1009  // }
1010  //}
1011  } else {
1012  if (opt.type == INJ) {
1013  OCP_DBL q = CalInjRate(myBulk, true);
1014  // for INJ well, maxRate has been switch to lbmols
1015  OCP_DBL tarRate = opt.maxRate;
1016  if (opt.reInj) {
1017  if (opt.injPhase == GAS)
1018  tarRate = WGIR;
1019  else if (opt.injPhase == WATER)
1020  tarRate = WWIR;
1021  }
1022 
1023  if (q > tarRate) {
1024  opt.optMode = opt.initOptMode;
1025  }
1026  else {
1027  opt.optMode = BHP_MODE;
1028  BHP = opt.maxBHP;
1029  }
1030  }
1031  else {
1032  OCP_DBL q = CalProdRate(myBulk, true);
1033  // cout << q << endl;
1034  if (q > opt.maxRate) {
1035  opt.optMode = opt.initOptMode;
1036  }
1037  else {
1038  opt.optMode = BHP_MODE;
1039  BHP = opt.minBHP;
1040  }
1041  }
1042  }
1043 }
1044 
1045 OCP_INT Well::CheckP(const Bulk& myBulk)
1046 {
1047  OCP_FUNCNAME;
1048  // 0 : all correct
1049  // 1 : negative P
1050  // 2 : outlimited P
1051  // 3 : crossflow happens
1052 
1053  if (BHP < 0) {
1054  cout << "### WARNING: Negative BHP " << BHP << endl;
1055  return 1;
1056  }
1057  for (USI p = 0; p < numPerf; p++) {
1058  if (perf[p].state == OPEN && perf[p].P < 0) {
1059 #ifdef _DEBUG
1060  cout << "### WARNING: Well " << name << " Perf[" << p << "].P = " << perf[p].P
1061  << endl;
1062 #endif // _DEBUG
1063  cout << "### WARNING: Negative perforation P " << perf[p].P << endl;
1064  return 1;
1065  }
1066  }
1067 
1068  if (opt.type == INJ) {
1069  if (opt.optMode != BHP_MODE && BHP > opt.maxBHP) {
1070 #if _DEBUG
1071  cout << "### WARNING: Well " << name << " switch to BHPMode" << endl;
1072 #endif
1073  opt.optMode = BHP_MODE;
1074  BHP = opt.maxBHP;
1075  return 2;
1076  }
1077  } else {
1078  if (opt.optMode != BHP_MODE && BHP < opt.minBHP) {
1079 #if _DEBUG
1080  cout << "### WARNING: Well " << name << " switch to BHPMode" << endl;
1081 #endif
1082  opt.optMode = BHP_MODE;
1083  BHP = opt.minBHP;
1084  return 2;
1085  }
1086  }
1087 
1088  return CheckCrossFlow(myBulk);
1089 }
1090 
1091 OCP_INT Well::CheckCrossFlow(const Bulk& myBulk)
1092 {
1093  OCP_FUNCNAME;
1094 
1095  OCP_USI k;
1096  bool flagC = true;
1097 
1098  if (opt.type == PROD) {
1099  for (USI p = 0; p < numPerf; p++) {
1100  k = perf[p].location;
1101  OCP_DBL minP = myBulk.P[k];
1102  if (perf[p].state == OPEN && minP < perf[p].P) {
1103  cout << std::left << std::setw(12) << name << " "
1104  << "Well P = " << perf[p].P << ", "
1105  << "Bulk P = " << minP << endl;
1106  perf[p].state = CLOSE;
1107  perf[p].multiplier = 0;
1108  flagC = false;
1109  break;
1110  } else if (perf[p].state == CLOSE && minP > perf[p].P) {
1111  perf[p].state = OPEN;
1112  perf[p].multiplier = 1;
1113  }
1114  }
1115  } else {
1116  for (USI p = 0; p < numPerf; p++) {
1117  k = perf[p].location;
1118  if (perf[p].state == OPEN && myBulk.P[k] > perf[p].P) {
1119  cout << std::left << std::setw(12) << name << " "
1120  << "Well P = " << perf[p].P << ", "
1121  << "Bulk P = " << myBulk.P[k] << endl;
1122  perf[p].state = CLOSE;
1123  perf[p].multiplier = 0;
1124  flagC = false;
1125  break;
1126  } else if (perf[p].state == CLOSE && myBulk.P[k] < perf[p].P) {
1127  perf[p].state = OPEN;
1128  perf[p].multiplier = 1;
1129  }
1130  }
1131  }
1132 
1133  bool flag = false;
1134  // check well -- if all perf are closed, open the depthest perf
1135  for (USI p = 0; p < numPerf; p++) {
1136  if (perf[p].state == OPEN) {
1137  flag = true;
1138  break;
1139  }
1140  }
1141 
1142  if (!flag) {
1143  // open the deepest perf
1144  perf.back().state = OPEN;
1145  perf.back().multiplier = 1;
1146  cout << "### WARNING: All perfs of " << name
1147  << " are closed! Open the last perf!\n";
1148  }
1149 
1150  if (!flagC) {
1151  // if crossflow happens, then corresponding perforation will be closed,
1152  // then the multiplier of perforation will be set to zero, so trans of well
1153  // should be recalculated!
1154  //
1155  // dG = ldG;
1156  CalTrans(myBulk);
1157  // CalFlux(myBulk);
1158  // CaldG(myBulk);
1159  // SmoothdG();
1160  // CheckOptMode(myBulk);
1161  return 3;
1162  }
1163 
1164  return 0;
1165 }
1166 
1167 void Well::AllocateMat(LinearSystem& myLS) const
1168 {
1169  OCP_FUNCNAME;
1170 
1171  for (USI p = 0; p < numPerf; p++) {
1172  myLS.rowCapacity[perf[p].location]++;
1173  }
1174 }
1175 
1176 void Well::SetupWellBulk(Bulk& myBulk) const
1177 {
1178  // Attention that a bulk can only be penetrated by one well now!
1179  for (USI p = 0; p < numPerf; p++) {
1180  myBulk.wellBulkId.push_back(perf[p].location);
1181  }
1182 }
1183 
1184 void Well::ShowPerfStatus(const Bulk& myBulk) const
1185 {
1186  OCP_FUNCNAME;
1187 
1188  cout << fixed;
1189  cout << "----------------------------" << endl;
1190  cout << name << ": " << opt.optMode << " " << setprecision(3) << BHP << endl;
1191  for (USI p = 0; p < numPerf; p++) {
1192  vector<OCP_DBL> Qitmp(perf[p].qi_lbmol);
1193  // OCP_DBL qt = Dnorm1(myBulk.numCom, &Qitmp[0]);
1194  OCP_USI n = perf[p].location;
1195  cout << setw(3) << p << " " << perf[p].state << " " << setw(6)
1196  << perf[p].location << " " << setw(2) << perf[p].I + 1 << " " << setw(2)
1197  << perf[p].J + 1 << " " << setw(2) << perf[p].K + 1 << " " << setw(10)
1198  << setprecision(6) << perf[p].WI << " " // ccf
1199  << setprecision(3) << perf[p].radius << " " // ccf
1200  << setw(8) << setprecision(4) << perf[p].kh << " " // kh
1201  << setw(8) << setprecision(2) << perf[p].depth << " " // depth
1202  << setprecision(3) << perf[p].P << " " // Pp
1203  << setw(10) << setprecision(3) << myBulk.P[n] << " " // Pb
1204  << setw(6) << setprecision(3) << dG[p] << " " // dG
1205  << setw(8) << perf[p].qi_lbmol[myBulk.numCom - 1] << " " << setw(6) << setprecision(6)
1206  << myBulk.S[n * myBulk.numPhase + 0] << " " << setw(6) << setprecision(6)
1207  << myBulk.S[n * myBulk.numPhase + 1] << " " << setw(6) << setprecision(6)
1208  << myBulk.S[n * myBulk.numPhase + 2] << endl;
1209  }
1210 }
1211 
1213 // IMPEC
1215 
1216 
1217 void Well::CalCFL(const Bulk& myBulk, const OCP_DBL& dt) const
1218 {
1219  OCP_FUNCNAME;
1220 
1221  if (opt.type == PROD) {
1222  const USI np = myBulk.numPhase;
1223  for (USI p = 0; p < numPerf; p++) {
1224  if (perf[p].state == OPEN) {
1225  OCP_USI k = perf[p].location;
1226 
1227  for (USI j = 0; j < np; j++) {
1228  myBulk.cfl[k * np + j] += fabs(perf[p].qj_ft3[j]) * dt;
1229  }
1230  }
1231  }
1232  }
1233 }
1234 
1235 void Well::MassConserveIMPEC(Bulk& myBulk, const OCP_DBL& dt) const
1236 {
1237  OCP_FUNCNAME;
1238 
1239  const USI nc = myBulk.numCom;
1240 
1241  for (USI p = 0; p < numPerf; p++) {
1242  OCP_USI k = perf[p].location;
1243  for (USI i = 0; i < nc; i++) {
1244  myBulk.Ni[k * nc + i] -= perf[p].qi_lbmol[i] * dt;
1245  }
1246  }
1247 }
1248 
1249 void Well::AssembleMatINJ_IMPEC(const Bulk& myBulk, LinearSystem& myLS,
1250  const OCP_DBL& dt) const
1251 {
1252  OCP_FUNCNAME;
1253 
1254  const USI nc = myBulk.numCom;
1255  const OCP_USI wId = myLS.dim;
1256  // important !
1257  myLS.dim++;
1258 
1259  for (USI p = 0; p < numPerf; p++) {
1260  OCP_USI k = perf[p].location;
1261 
1262  OCP_DBL Vfi_zi = 0;
1263  for (USI i = 0; i < nc; i++) {
1264  Vfi_zi += myBulk.vfi[k * nc + i] * opt.zi[i];
1265  }
1266 
1267  OCP_DBL valw = dt * perf[p].xi * perf[p].transINJ;
1268  OCP_DBL bw = valw * dG[p];
1269  OCP_DBL valb = valw * Vfi_zi;
1270  OCP_DBL bb = valb * dG[p];
1271 
1272  // Bulk to Well
1273 
1274  // diag
1275  USI ptr = myLS.diagPtr[k];
1276  myLS.val[k][ptr] += valb;
1277  // off diag
1278  myLS.colId[k].push_back(wId);
1279  myLS.val[k].push_back(-valb);
1280  // b
1281  myLS.b[k] += bb;
1282 
1283  // Well to Bulk
1284  switch (opt.optMode) {
1285  case RATE_MODE:
1286  case ORATE_MODE:
1287  case GRATE_MODE:
1288  case WRATE_MODE:
1289  case LRATE_MODE:
1290  // diag
1291  myLS.diagVal[wId] += valw;
1292  // off diag
1293  myLS.colId[wId].push_back(k);
1294  myLS.val[wId].push_back(-valw);
1295  // b
1296  myLS.b[wId] -= bw;
1297  break;
1298  case BHP_MODE:
1299  myLS.colId[wId].push_back(k);
1300  myLS.val[wId].push_back(0);
1301  break;
1302  default:
1303  OCP_ABORT("Wrong well option mode!");
1304  }
1305  }
1306 
1307  // Well Self
1308  assert(myLS.val[wId].size() == numPerf);
1309  // the order of perforation is not necessarily in order
1310  switch (opt.optMode) {
1311  case RATE_MODE:
1312  case ORATE_MODE:
1313  case GRATE_MODE:
1314  case WRATE_MODE:
1315  case LRATE_MODE:
1316  // diag
1317  myLS.colId[wId].push_back(wId);
1318  myLS.diagPtr[wId] = numPerf;
1319  myLS.val[wId].push_back(myLS.diagVal[wId]);
1320  // b
1321  myLS.b[wId] += dt * opt.maxRate;
1322  break;
1323  case BHP_MODE:
1324  // diag
1325  myLS.colId[wId].push_back(wId);
1326  myLS.diagPtr[wId] = numPerf;
1327  myLS.val[wId].push_back(dt);
1328  // b
1329  myLS.b[wId] += dt * opt.maxBHP;
1330  // u initial value
1331  myLS.u[wId] = opt.maxBHP;
1332  break;
1333  default:
1334  OCP_ABORT("Wrong well option mode!");
1335  }
1336 }
1337 
1338 void Well::AssembleMatPROD_IMPEC(const Bulk& myBulk, LinearSystem& myLS,
1339  const OCP_DBL& dt) const
1340 {
1341  OCP_FUNCNAME;
1342 
1343  const USI np = myBulk.numPhase;
1344  const USI nc = myBulk.numCom;
1345  const OCP_USI wId = myLS.dim;
1346  // important !
1347  myLS.dim++;
1348 
1349  // Set Prod Weight
1350  if (opt.optMode != BHP_MODE) CalProdWeight(myBulk);
1351 
1352  for (USI p = 0; p < numPerf; p++) {
1353  OCP_USI n = perf[p].location;
1354 
1355  OCP_DBL valb = 0;
1356  OCP_DBL bb = 0;
1357  OCP_DBL valw = 0;
1358  OCP_DBL bw = 0;
1359 
1360  for (USI j = 0; j < np; j++) {
1361  if (!myBulk.phaseExist[n * np + j]) continue;
1362 
1363  OCP_DBL tempb = 0;
1364  OCP_DBL tempw = 0;
1365 
1366  for (USI i = 0; i < nc; i++) {
1367  tempb += myBulk.vfi[n * nc + i] * myBulk.xij[n * np * nc + j * nc + i];
1368  tempw += prodWeight[i] * myBulk.xij[n * np * nc + j * nc + i];
1369  }
1370  OCP_DBL trans = dt * perf[p].transj[j] * myBulk.xi[n * np + j];
1371  valb += tempb * trans;
1372  valw += tempw * trans;
1373 
1374  OCP_DBL dP = dG[p] - myBulk.Pc[n * np + j];
1375  bb += tempb * trans * dP;
1376  bw += tempw * trans * dP;
1377  }
1378 
1379  // Bulk to Well
1380  // diag
1381  USI ptr = myLS.diagPtr[n];
1382  myLS.val[n][ptr] += valb;
1383  // off diag
1384  myLS.colId[n].push_back(wId);
1385  myLS.val[n].push_back(-valb);
1386  // b
1387  myLS.b[n] += bb;
1388 
1389  // Well to Bulk
1390  switch (opt.optMode) {
1391  case RATE_MODE:
1392  case ORATE_MODE:
1393  case GRATE_MODE:
1394  case WRATE_MODE:
1395  case LRATE_MODE:
1396  // diag !!! attention! sign is -
1397  myLS.diagVal[wId] -= valw;
1398  // off diag
1399  myLS.colId[wId].push_back(n);
1400  myLS.val[wId].push_back(valw);
1401  // b
1402  myLS.b[wId] += bw;
1403  break;
1404  case BHP_MODE:
1405  // off diag
1406  myLS.colId[wId].push_back(n);
1407  myLS.val[wId].push_back(0);
1408  break;
1409  default:
1410  OCP_ABORT("Wrong well option mode!");
1411  }
1412  }
1413 
1414  // Well Self
1415  assert(myLS.val[wId].size() == numPerf);
1416  // the order of perforation is not necessarily in order
1417  switch (opt.optMode) {
1418  case RATE_MODE:
1419  case ORATE_MODE:
1420  case GRATE_MODE:
1421  case WRATE_MODE:
1422  case LRATE_MODE:
1423  // diag
1424  myLS.colId[wId].push_back(wId);
1425  myLS.diagPtr[wId] = numPerf;
1426  myLS.val[wId].push_back(myLS.diagVal[wId]);
1427  // b
1428  myLS.b[wId] += dt * opt.maxRate;
1429  break;
1430  case BHP_MODE:
1431  // diag
1432  myLS.colId[wId].push_back(wId);
1433  myLS.diagPtr[wId] = numPerf;
1434  myLS.val[wId].push_back(dt);
1435  // b
1436  myLS.b[wId] += dt * opt.minBHP;
1437  // u initial value
1438  myLS.u[wId] = opt.minBHP;
1439  break;
1440  default:
1441  OCP_ABORT("Wrong well option mode!");
1442  }
1443 }
1444 
1445 void Well::AssembleMatReinjection_IMPEC(const Bulk& myBulk, LinearSystem& myLS,
1446  const OCP_DBL& dt, const vector<Well>& allWell,
1447  const vector<USI>& injId) const
1448 {
1449  // find Open injection well under Rate control
1450  vector<OCP_USI> tarId;
1451  for (USI w = 0; w < injId.size(); w++) {
1452  if (allWell[w].WellState() && allWell[w].opt.optMode != BHP_MODE)
1453  tarId.push_back(allWell[w].wEId + myBulk.numBulk);
1454  }
1455 
1456  USI tlen = tarId.size();
1457  if (tlen > 0) {
1458  const OCP_DBL factor = allWell[injId[0]].opt.factor * dt;
1459  // cout << "Factor(assemble): " << allWell[injId[0]].opt.factor << endl;
1460  const OCP_USI prodId = wEId + myBulk.numBulk;
1461  USI np = myBulk.numPhase;
1462  OCP_USI n, bId;
1463  OCP_DBL tmp, valb;
1464  OCP_DBL valw = 0;
1465  OCP_DBL rhsw = 0;
1466  OCP_USI tar;
1467  USI tarsize;
1468  for (USI p = 0; p < numPerf; p++) {
1469  n = perf[p].location;
1470  valb = 0;
1471 
1472  for (USI j = 0; j < np; j++) {
1473  bId = n * np + j;
1474  if (myBulk.phaseExist[bId]) {
1475  tmp = perf[p].transj[j] * myBulk.xi[bId];
1476  valb += tmp;
1477  rhsw += tmp * (myBulk.Pc[bId] - dG[p]);
1478  }
1479  }
1480  valb *= factor;
1481  // Insert bulk val into equations in ascending order
1482  for (USI t = 0; t < tlen; t++) {
1483  tar = tarId[t];
1484  tarsize = myLS.colId[tar].size();
1485  USI i;
1486  for (i = 0; i < tarsize; i++) {
1487  if (n < myLS.colId[tar][i]) {
1488  // insert
1489  // attention that bulk id is less than well id
1490  myLS.colId[tar].insert(myLS.colId[tar].begin() + i, n);
1491  myLS.val[tar].insert(myLS.val[tar].begin() + i, -valb);
1492  myLS.diagPtr[tar]++;
1493  break;
1494  }
1495  }
1496  }
1497  valw += valb;
1498  }
1499  rhsw *= factor;
1500  // insert prod well var and rhs into equations in ascending order
1501  for (USI t = 0; t < tlen; t++) {
1502  tar = tarId[t];
1503  tarsize = myLS.colId[tar].size();
1504  myLS.b[tar] += rhsw; // rhsw
1505  USI i;
1506  for (i = 0; i < tarsize; i++) {
1507  if (prodId < myLS.colId[tar][i]) {
1508  // insert
1509  myLS.colId[tar].insert(myLS.colId[tar].begin() + i, prodId);
1510  myLS.val[tar].insert(myLS.val[tar].begin() + i, valw);
1511  if (i <= myLS.diagPtr[tar]) myLS.diagPtr[tar]++;
1512  break;
1513  }
1514  }
1515  if (i == tarsize) {
1516  // pushback
1517  myLS.colId[tar].push_back(prodId);
1518  myLS.val[tar].push_back(valw);
1519  }
1520  }
1521  }
1522 }
1523 
1525 // FIM
1527 
1528 void Well::AssembleMatINJ_FIM(const Bulk& myBulk, LinearSystem& myLS,
1529  const OCP_DBL& dt) const
1530 {
1531  OCP_FUNCNAME;
1532 
1533  const OCP_USI wId = myLS.dim;
1534  // important !
1535  myLS.dim++;
1536 
1537  const USI np = myBulk.numPhase;
1538  const USI nc = myBulk.numCom;
1539  const USI ncol = nc + 1;
1540  const USI ncol2 = np * nc + np;
1541  const USI bsize = ncol * ncol;
1542  const USI bsize2 = ncol * ncol2;
1543 
1544  OCP_DBL mu, muP;
1545  OCP_DBL dP;
1546  OCP_DBL transIJ;
1547 
1548  OCP_USI n_np_j;
1549 
1550  vector<OCP_DBL> bmat(bsize, 0);
1551  vector<OCP_DBL> bmat2(bsize, 0);
1552  vector<OCP_DBL> dQdXpB(bsize, 0);
1553  vector<OCP_DBL> dQdXpW(bsize, 0);
1554  vector<OCP_DBL> dQdXsB(bsize2, 0);
1555 
1556  for (USI p = 0; p < numPerf; p++) {
1557  const OCP_USI n = perf[p].location;
1558  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
1559  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
1560  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
1561 
1562  dP = myBulk.P[n] - BHP - dG[p];
1563 
1564  for (USI j = 0; j < np; j++) {
1565  n_np_j = n * np + j;
1566  if (!myBulk.phaseExist[n_np_j]) continue;
1567 
1568  mu = myBulk.mu[n_np_j];
1569  muP = myBulk.muP[n_np_j];
1570 
1571  for (USI i = 0; i < nc; i++) {
1572  // dQ / dP
1573  transIJ = perf[p].transj[j] * perf[p].xi * opt.zi[i];
1574  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
1575  dQdXpW[(i + 1) * ncol] += -transIJ;
1576 
1577  // dQ / dS
1578  for (USI k = 0; k < np; k++) {
1579  dQdXsB[(i + 1) * ncol2 + k] +=
1580  CONV1 * perf[p].WI * perf[p].multiplier * perf[p].xi *
1581  opt.zi[i] * myBulk.dKr_dS[n_np_j * np + k] * dP / mu;
1582  }
1583  // dQ / dxij
1584  for (USI k = 0; k < nc; k++) {
1585  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] +=
1586  -transIJ * dP / mu * myBulk.mux[n_np_j * nc + k];
1587  }
1588  }
1589  }
1590  // Bulk to Well
1591  bmat = dQdXpB;
1592  //myDABpCp1(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data(), &flagB[0], nc - 1);
1593  //myDABpCp(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
1594  //myDABpC(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data());
1595  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], 1,
1596  bmat.data());
1597  Dscalar(bsize, dt, bmat.data());
1598  // Add
1599  USI ptr = myLS.diagPtr[n];
1600  for (USI i = 0; i < bsize; i++) {
1601  myLS.val[n][ptr * bsize + i] += bmat[i];
1602  }
1603  // Insert
1604  bmat = dQdXpW;
1605  Dscalar(bsize, dt, bmat.data());
1606  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
1607  myLS.colId[n].push_back(wId);
1608 
1609  // Well
1610  switch (opt.optMode) {
1611  case RATE_MODE:
1612  case ORATE_MODE:
1613  case GRATE_MODE:
1614  case WRATE_MODE:
1615  case LRATE_MODE:
1616  // Diag
1617  fill(bmat.begin(), bmat.end(), 0.0);
1618  for (USI i = 0; i < nc; i++) {
1619  bmat[0] += dQdXpW[(i + 1) * ncol];
1620  bmat[(i + 1) * ncol + i + 1] = 1;
1621  }
1622  for (USI i = 0; i < bsize; i++) {
1623  myLS.diagVal[wId * bsize + i] += bmat[i];
1624  }
1625 
1626  // OffDiag
1627  bmat = dQdXpB;
1628  //myDABpCp1(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data(), &flagB[0], nc - 1);
1629  //myDABpCp(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
1630  //myDABpC(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data());
1631  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(),
1632  &myBulk.dSec_dPri[n * bsize2], 1, bmat.data());
1633  fill(bmat2.begin(), bmat2.end(), 0.0);
1634  for (USI i = 0; i < nc; i++) {
1635  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
1636  }
1637  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
1638  myLS.colId[wId].push_back(n);
1639  break;
1640 
1641  case BHP_MODE:
1642  // Diag
1643  fill(bmat.begin(), bmat.end(), 0.0);
1644  for (USI i = 0; i < nc + 1; i++) {
1645  bmat[i * ncol + i] = 1;
1646  }
1647  // Add
1648  for (USI i = 0; i < bsize; i++) {
1649  myLS.diagVal[wId * bsize + i] += bmat[i];
1650  }
1651  // OffDiag
1652  fill(bmat.begin(), bmat.end(), 0.0);
1653  // Insert
1654  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
1655  myLS.colId[wId].push_back(n);
1656  // Solution
1657  // myLS.u[wId * ncol] = opt.maxBHP - BHP;
1658  break;
1659 
1660  default:
1661  OCP_ABORT("Wrong Well Opt mode!");
1662  break;
1663  }
1664  }
1665  assert(myLS.val[wId].size() == numPerf * bsize);
1666  // Well self
1667  myLS.colId[wId].push_back(wId);
1668  myLS.diagPtr[wId] = numPerf;
1669  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
1670  myLS.diagVal.data() + wId * bsize + bsize);
1671 }
1672 
1673 void Well::AssembleMatPROD_FIM(const Bulk& myBulk, LinearSystem& myLS,
1674  const OCP_DBL& dt) const
1675 {
1676  OCP_FUNCNAME;
1677 
1678  const OCP_USI wId = myLS.dim;
1679  // important !
1680  myLS.dim++;
1681 
1682  const USI np = myBulk.numPhase;
1683  const USI nc = myBulk.numCom;
1684  const USI ncol = nc + 1;
1685  const USI ncol2 = np * nc + np;
1686  const USI bsize = ncol * ncol;
1687  const USI bsize2 = ncol * ncol2;
1688 
1689  OCP_DBL xij, xi, mu, muP, xiP;
1690  OCP_DBL dP;
1691  OCP_DBL transIJ;
1692  OCP_DBL tmp;
1693 
1694  OCP_USI n_np_j;
1695 
1696  vector<OCP_DBL> bmat(bsize, 0);
1697  vector<OCP_DBL> bmat2(bsize, 0);
1698  vector<OCP_DBL> dQdXpB(bsize, 0);
1699  vector<OCP_DBL> dQdXpW(bsize, 0);
1700  vector<OCP_DBL> dQdXsB(bsize2, 0);
1701 
1702 
1703  // Set Prod Weight
1704  if (opt.optMode != BHP_MODE) CalProdWeight(myBulk);
1705 
1706  for (USI p = 0; p < numPerf; p++) {
1707  const OCP_USI n = perf[p].location;
1708  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
1709  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
1710  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
1711 
1712  for (USI j = 0; j < np; j++) {
1713  n_np_j = n * np + j;
1714  if (!myBulk.phaseExist[n_np_j]) continue;
1715 
1716  dP = myBulk.Pj[n_np_j] - BHP - dG[p];
1717  xi = myBulk.xi[n_np_j];
1718  mu = myBulk.mu[n_np_j];
1719  muP = myBulk.muP[n_np_j];
1720  xiP = myBulk.xiP[n_np_j];
1721 
1722  for (USI i = 0; i < nc; i++) {
1723  xij = myBulk.xij[n_np_j * nc + i];
1724  // dQ / dP
1725  transIJ = perf[p].transj[j] * xi * xij;
1726  dQdXpB[(i + 1) * ncol] +=
1727  transIJ * (1 - dP * muP / mu) + dP * perf[p].transj[j] * xij * xiP;
1728  dQdXpW[(i + 1) * ncol] += -transIJ;
1729 
1730  // if (dQdXpW[ncol] == 0) {
1731  // cout << name << " perf " << p << " " << j << endl;
1732  //}
1733 
1734  // dQ / dS
1735  for (USI k = 0; k < np; k++) {
1736  tmp = CONV1 * perf[p].WI * perf[p].multiplier * dP / mu * xi * xij *
1737  myBulk.dKr_dS[n_np_j * np + k];
1738  // capillary pressure
1739  tmp += transIJ * myBulk.dPcj_dS[n_np_j * np + k];
1740  dQdXsB[(i + 1) * ncol2 + k] += tmp;
1741  }
1742  // dQ / dCij
1743  for (USI k = 0; k < nc; k++) {
1744  tmp = dP * perf[p].transj[j] * xij *
1745  (myBulk.xix[n_np_j * nc + k] -
1746  xi / mu * myBulk.mux[n_np_j * nc + k]);
1747  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] += tmp;
1748  }
1749  dQdXsB[(i + 1) * ncol2 + np + j * nc + i] += perf[p].transj[j] * xi * dP;
1750  }
1751  }
1752  // Bulk to Well
1753  bmat = dQdXpB;
1754  //myDABpCp1(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data(), &flagB[0], nc - 1);
1755  //myDABpCp(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
1756  //myDABpC(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data());
1757  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], 1,
1758  bmat.data());
1759 
1760  Dscalar(bsize, dt, bmat.data());
1761  // Add
1762  USI ptr = myLS.diagPtr[n];
1763  for (USI i = 0; i < bsize; i++) {
1764  myLS.val[n][ptr * bsize + i] += bmat[i];
1765  }
1766  // Insert
1767  bmat = dQdXpW;
1768  Dscalar(bsize, dt, bmat.data());
1769  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
1770  myLS.colId[n].push_back(wId);
1771 
1772  // Well
1773  switch (opt.optMode) {
1774  case RATE_MODE:
1775  case ORATE_MODE:
1776  case GRATE_MODE:
1777  case WRATE_MODE:
1778  case LRATE_MODE:
1779  // Diag
1780  fill(bmat.begin(), bmat.end(), 0.0);
1781  for (USI i = 0; i < nc; i++) {
1782  bmat[0] += dQdXpW[(i + 1) * ncol] * prodWeight[i];
1783  bmat[(i + 1) * ncol + i + 1] = 1;
1784  }
1785  for (USI i = 0; i < bsize; i++) {
1786  myLS.diagVal[wId * bsize + i] += bmat[i];
1787  }
1788 
1789  // OffDiag
1790  bmat = dQdXpB;
1791  //myDABpCp1(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data(), &flagB[0], nc - 1);
1792  //myDABpCp(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
1793  //myDABpC(ncol, ncol2, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data());
1794  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(),
1795  &myBulk.dSec_dPri[n * bsize2], 1, bmat.data());
1796  fill(bmat2.begin(), bmat2.end(), 0.0);
1797  for (USI i = 0; i < nc; i++) {
1798  Daxpy(ncol, prodWeight[i], bmat.data() + (i + 1) * ncol,
1799  bmat2.data());
1800  }
1801  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
1802  myLS.colId[wId].push_back(n);
1803  break;
1804 
1805  case BHP_MODE:
1806  // Diag
1807  fill(bmat.begin(), bmat.end(), 0.0);
1808  for (USI i = 0; i < nc + 1; i++) {
1809  bmat[i * ncol + i] = 1;
1810  }
1811  // Add
1812  for (USI i = 0; i < bsize; i++) {
1813  myLS.diagVal[wId * bsize + i] += bmat[i];
1814  }
1815  // OffDiag
1816  fill(bmat.begin(), bmat.end(), 0.0);
1817  // Insert
1818  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
1819  myLS.colId[wId].push_back(n);
1820  // Solution
1821  // myLS.u[wId * ncol] = opt.minBHP - BHP;
1822  break;
1823 
1824  default:
1825  OCP_ABORT("Wrong Well Opt mode!");
1826  break;
1827  }
1828  }
1829  assert(myLS.val[wId].size() == numPerf * bsize);
1830  // Well self
1831  myLS.colId[wId].push_back(wId);
1832  myLS.diagPtr[wId] = numPerf;
1833  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
1834  myLS.diagVal.data() + wId * bsize + bsize);
1835 
1836  // test
1837  // for (USI i = 0; i < bsize; i++) {
1838  // cout << myLS.diagVal[wId * bsize + i] << endl;
1839  //}
1840 }
1841 
1842 void Well::AssembleMatReinjection_FIM(const Bulk& myBulk, LinearSystem& myLS,
1843  const OCP_DBL& dt, const vector<Well>& allWell,
1844  const vector<USI>& injId) const
1845 {
1846  // find Open injection well under Rate control
1847  vector<OCP_USI> tarId;
1848  for (USI w = 0; w < injId.size(); w++) {
1849  if (allWell[w].WellState() && allWell[w].opt.optMode != BHP_MODE)
1850  tarId.push_back(allWell[w].wEId + myBulk.numBulk);
1851  }
1852 
1853  USI tlen = tarId.size();
1854  if (tlen > 0) {
1855  OCP_USI tar;
1856  USI tarsize;
1857  const OCP_DBL factor = allWell[injId[0]].opt.factor;
1858  const OCP_USI prodId = wEId + myBulk.numBulk;
1859 
1860  // cout << "Factor(assemble): " << factor << endl;
1861 
1862  const USI np = myBulk.numPhase;
1863  const USI nc = myBulk.numCom;
1864  const USI ncol = nc + 1;
1865  const USI ncol2 = np * nc + np;
1866  const USI bsize = ncol * ncol;
1867  const USI bsize2 = ncol * ncol2;
1868 
1869  OCP_DBL xij, xi, mu, muP, xiP, dP, transIJ, tmp;
1870  OCP_USI n_np_j;
1871 
1872  vector<OCP_DBL> bmat(bsize, 0);
1873  vector<OCP_DBL> bmat2(bsize, 0);
1874  vector<OCP_DBL> tmpMat(bsize, 0);
1875  vector<OCP_DBL> dQdXpB(bsize, 0);
1876  vector<OCP_DBL> dQdXpW(bsize, 0);
1877  vector<OCP_DBL> dQdXsB(bsize2, 0);
1878 
1879  for (USI p = 0; p < numPerf; p++) {
1880  OCP_USI n = perf[p].location;
1881  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
1882  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
1883  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
1884 
1885  for (USI j = 0; j < np; j++) {
1886  n_np_j = n * np + j;
1887  if (!myBulk.phaseExist[n_np_j]) continue;
1888 
1889  dP = myBulk.Pj[n_np_j] - BHP - dG[p];
1890  xi = myBulk.xi[n_np_j];
1891  mu = myBulk.mu[n_np_j];
1892  muP = myBulk.muP[n_np_j];
1893  xiP = myBulk.xiP[n_np_j];
1894 
1895  for (USI i = 0; i < nc; i++) {
1896  xij = myBulk.xij[n_np_j * nc + i];
1897  // dQ / dP
1898  transIJ = perf[p].transj[j] * xi * xij;
1899  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu) +
1900  dP * perf[p].transj[j] * xij * xiP;
1901  dQdXpW[(i + 1) * ncol] += -transIJ;
1902 
1903  // dQ / dS
1904  for (USI k = 0; k < np; k++) {
1905  tmp = CONV1 * perf[p].WI * perf[p].multiplier * dP / mu * xi *
1906  xij * myBulk.dKr_dS[n * np * np + j * np + k];
1907  // capillary pressure
1908  tmp += transIJ * myBulk.dPcj_dS[n * np * np + j * np + k];
1909  dQdXsB[(i + 1) * ncol2 + k] += tmp;
1910  }
1911  // dQ / dCij
1912  for (USI k = 0; k < nc; k++) {
1913  tmp = dP * perf[p].transj[j] * xij *
1914  (myBulk.xix[n_np_j * nc + k] -
1915  xi / mu * myBulk.mux[n_np_j * nc + k]);
1916  if (k == i) {
1917  tmp += perf[p].transj[j] * xi * dP;
1918  }
1919  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] += tmp;
1920  }
1921  }
1922  }
1923 
1924  // for Prod Well, be careful!
1925  for (USI i = 0; i < nc; i++) {
1926  // tmpMat[0] -= dQdXpW[(i + 1) * ncol] * factor;
1927  tmpMat[0] += dQdXpW[(i + 1) * ncol] * factor;
1928  }
1929 
1930  // for perf(bulk) of Prod Well
1931  bmat = dQdXpB;
1932  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2],
1933  1, bmat.data());
1934  fill(bmat2.begin(), bmat2.end(), 0.0);
1935  for (USI i = 0; i < nc; i++) {
1936  // becareful '-' before factor
1937  // Daxpy(ncol, -factor, bmat.data() + (i + 1) * ncol, bmat2.data());
1938  Daxpy(ncol, factor, bmat.data() + (i + 1) * ncol, bmat2.data());
1939  }
1940 
1941  // Insert bulk val into equations in ascending order
1942  for (USI t = 0; t < tlen; t++) {
1943  tar = tarId[t];
1944  tarsize = myLS.colId[tar].size();
1945  USI i;
1946  for (i = 0; i < tarsize; i++) {
1947  if (n < myLS.colId[tar][i]) {
1948  // insert
1949  // attention that bulk id is less than well id
1950  myLS.colId[tar].insert(myLS.colId[tar].begin() + i, n);
1951  myLS.val[tar].insert(myLS.val[tar].begin() + i * bsize,
1952  bmat.begin(), bmat.end());
1953  myLS.diagPtr[tar]++;
1954  break;
1955  }
1956  }
1957  }
1958  }
1959  // insert prod well var into equations in ascending order
1960  for (USI t = 0; t < tlen; t++) {
1961  tar = tarId[t];
1962  tarsize = myLS.colId[tar].size();
1963  USI i;
1964  for (i = 0; i < tarsize; i++) {
1965  if (prodId < myLS.colId[tar][i]) {
1966  // insert
1967  myLS.colId[tar].insert(myLS.colId[tar].begin() + i, prodId);
1968  myLS.val[tar].insert(myLS.val[tar].begin() + i * bsize,
1969  tmpMat.begin(), tmpMat.end());
1970  if (i <= myLS.diagPtr[tar]) myLS.diagPtr[tar]++;
1971  break;
1972  }
1973  }
1974  if (i == tarsize) {
1975  // insert into end
1976  myLS.colId[tar].push_back(prodId);
1977  myLS.val[tar].insert(myLS.val[tar].end(), tmpMat.begin(), tmpMat.end());
1978  }
1979  }
1980  }
1981 }
1982 
1983 void Well::CalResFIM(ResFIM& resFIM, const Bulk& myBulk, const OCP_DBL& dt,
1984  const OCP_USI& wId, const vector<Well>& allWell) const
1985 {
1986  OCP_FUNCNAME;
1987 
1988  // Well to Bulk
1989  const USI nc = myBulk.numCom;
1990  const USI len = nc + 1;
1991  OCP_USI k;
1992 
1993  for (USI p = 0; p < numPerf; p++) {
1994  k = perf[p].location;
1995  for (USI i = 0; i < nc; i++) {
1996  resFIM.res[k * len + 1 + i] += perf[p].qi_lbmol[i] * dt;
1997  }
1998  }
1999 
2000  // Well Self
2001  OCP_USI bId = (myBulk.numBulk + wId) * len;
2002  if (opt.type == INJ) {
2003  // Injection
2004  switch (opt.optMode) {
2005  case BHP_MODE:
2006  BHP = opt.maxBHP;
2007  resFIM.res[bId] = BHP - opt.maxBHP;
2008  break;
2009  case RATE_MODE:
2010  case ORATE_MODE:
2011  case GRATE_MODE:
2012  case WRATE_MODE:
2013  case LRATE_MODE:
2014  resFIM.res[bId] = opt.maxRate;
2015  for (USI i = 0; i < nc; i++) {
2016  resFIM.res[bId] += qi_lbmol[i];
2017  }
2018  if (opt.reInj) {
2019  for (auto& w : opt.connWell) {
2020  OCP_DBL tmp = 0;
2021  for (USI i = 0; i < nc; i++) {
2022  tmp += allWell[w].qi_lbmol[i];
2023  // tmp += opt.factor * allWell[w].qi_lbmol[i];
2024  // resFIM.res[bId] += opt.factor * allWell[w].qi_lbmol[i];
2025  }
2026  tmp *= opt.factor;
2027  resFIM.res[bId] += tmp;
2028  // cout << "Temp(INJ): " << tmp / opt.xiINJ / 1000 << endl;
2029  }
2030  // cout << "Factor(res) " << opt.factor << endl;
2031  }
2032  //cout << name << " " << resFIM.res[bId] << " " << opt.maxRate << " " <<
2033  // fabs(resFIM.res[bId] / opt.maxRate) << endl;
2034  resFIM.maxWellRelRes_mol =
2035  max(resFIM.maxWellRelRes_mol, fabs(resFIM.res[bId] / opt.maxRate));
2036  break;
2037  default:
2038  OCP_ABORT("Wrong well opt mode!");
2039  break;
2040  }
2041  } else {
2042  // Production
2043  switch (opt.optMode) {
2044  case BHP_MODE:
2045  BHP = opt.minBHP;
2046  resFIM.res[bId] = BHP - opt.minBHP;
2047  break;
2048  case RATE_MODE:
2049  case ORATE_MODE:
2050  case GRATE_MODE:
2051  case WRATE_MODE:
2052  case LRATE_MODE:
2053  resFIM.res[bId] = -opt.maxRate;
2054  for (USI i = 0; i < nc; i++) {
2055  resFIM.res[bId] += qi_lbmol[i] * prodWeight[i];
2056  }
2057  // cout << "Temp(Prod): " << tmp << endl;
2058  //cout << name << " " << resFIM.res[bId] << " " << opt.maxRate << " "
2059  // << fabs(resFIM.res[bId] / opt.maxRate) << endl;
2060  resFIM.maxWellRelRes_mol =
2061  max(resFIM.maxWellRelRes_mol, fabs(resFIM.res[bId] / opt.maxRate));
2062  break;
2063  default:
2064  OCP_ABORT("Wrong well opt mode!");
2065  break;
2066  }
2067  }
2068 }
2069 
2070 void Well::ShowRes(const OCP_USI& wId, const vector<OCP_DBL>& res, const Bulk& myBulk) const
2071 {
2072  // Well to Bulk
2073  const USI nc = myBulk.numCom;
2074  const USI len = nc + 1;
2075 
2076  OCP_USI bId = (myBulk.numBulk + wId) * len;
2077  cout << name << " " << res[bId] << " ";
2078  // Well Self
2079  if (opt.type == INJ) {
2080  // Injection
2081  switch (opt.optMode) {
2082  case BHP_MODE:
2083  cout << "BHPMode" << endl;
2084  break;
2085  case RATE_MODE:
2086  case ORATE_MODE:
2087  case GRATE_MODE:
2088  case WRATE_MODE:
2089  case LRATE_MODE:
2090  cout << opt.maxRate << " " << fabs(res[bId] / opt.maxRate) << endl;
2091  break;
2092  default:
2093  OCP_ABORT("Wrong well opt mode!");
2094  break;
2095  }
2096  }
2097  else {
2098  // Production
2099  switch (opt.optMode) {
2100  case BHP_MODE:
2101  cout << "BHPMode" << endl;
2102  break;
2103  case RATE_MODE:
2104  case ORATE_MODE:
2105  case GRATE_MODE:
2106  case WRATE_MODE:
2107  case LRATE_MODE:
2108  cout << opt.maxRate << " " << fabs(res[bId] / opt.maxRate) << endl;
2109  break;
2110  default:
2111  OCP_ABORT("Wrong well opt mode!");
2112  break;
2113  }
2114  }
2115 }
2116 
2117 
2119 // FIM(new)
2121 
2122 
2123 void Well::AssembleMatINJ_FIM_new(const Bulk& myBulk, LinearSystem& myLS,
2124  const OCP_DBL& dt) const
2125 {
2126  OCP_FUNCNAME;
2127 
2128  const OCP_USI wId = myLS.dim;
2129  // important !
2130  myLS.dim++;
2131 
2132  const USI np = myBulk.numPhase;
2133  const USI nc = myBulk.numCom;
2134  const USI nch = nc - 1;
2135  const USI ncol = nc + 1;
2136  const USI ncol2 = np * nc + np;
2137  const USI bsize = ncol * ncol;
2138  const USI bsize2 = ncol * ncol2;
2139 
2140  OCP_DBL mu, muP;
2141  OCP_DBL dP;
2142  OCP_DBL transIJ;
2143 
2144  OCP_USI n_np_j;
2145 
2146  vector<OCP_DBL> bmat(bsize, 0);
2147  vector<OCP_DBL> bmat2(bsize, 0);
2148  vector<OCP_DBL> dQdXpB(bsize, 0);
2149  vector<OCP_DBL> dQdXpW(bsize, 0);
2150  vector<OCP_DBL> dQdXsB(bsize2, 0);
2151  vector<bool> phaseExistB(np, false);
2152  vector<USI> pEnumComB(np, 0);
2153  USI ncolB;
2154 
2155 
2156  for (USI p = 0; p < numPerf; p++) {
2157  const OCP_USI n = perf[p].location;
2158  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2159  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2160  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2161 
2162  dP = myBulk.P[n] - BHP - dG[p];
2163 
2164  const USI npB = myBulk.phaseNum[n] + 1; ncolB = npB;
2165  for (USI j = 0; j < np; j++) {
2166  phaseExistB[j] = myBulk.phaseExist[n * np + j];
2167  pEnumComB[j] = myBulk.pEnumCom[n * np + j];
2168  ncolB += pEnumComB[j];
2169  }
2170 
2171 
2172  USI jxB = npB;
2173  for (USI j = 0; j < np; j++) {
2174 
2175  if (!phaseExistB[j]) {
2176  jxB += pEnumComB[j];
2177  continue;
2178  }
2179 
2180  n_np_j = n * np + j;
2181  mu = myBulk.mu[n_np_j];
2182  muP = myBulk.muP[n_np_j];
2183 
2184  for (USI i = 0; i < nc; i++) {
2185  // dQ / dP
2186  transIJ = perf[p].transj[j] * perf[p].xi * opt.zi[i];
2187  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
2188  dQdXpW[(i + 1) * ncol] += -transIJ;
2189 
2190  // dQ / dS
2191  USI j1B = 0;
2192  for (USI j1 = 0; j1 < np; j1++) {
2193  if (phaseExistB[j1]) {
2194  dQdXsB[(i + 1) * ncolB + j1B] +=
2195  CONV1 * perf[p].WI * perf[p].multiplier * perf[p].xi *
2196  opt.zi[i] * myBulk.dKr_dS[n_np_j * np + j1] * dP / mu;
2197  j1B++;
2198  }
2199  }
2200 
2201  // dQ / dxij
2202  for (USI k = 0; k < pEnumComB[j]; k++) {
2203  dQdXsB[(i + 1) * ncolB + jxB + k] +=
2204  -transIJ * dP / mu * myBulk.mux[n_np_j * nc + k];
2205  }
2206  }
2207  jxB += pEnumComB[j];
2208  }
2209  // Bulk to Well
2210  bmat = dQdXpB;
2211  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2212  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2213  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1,
2214  bmat.data());
2215  Dscalar(bsize, dt, bmat.data());
2216  // Add
2217  USI ptr = myLS.diagPtr[n];
2218  for (USI i = 0; i < bsize; i++) {
2219  myLS.val[n][ptr * bsize + i] += bmat[i];
2220  }
2221  // Insert
2222  bmat = dQdXpW;
2223  Dscalar(bsize, dt, bmat.data());
2224  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
2225  myLS.colId[n].push_back(wId);
2226 
2227  // Well
2228  switch (opt.optMode) {
2229  case RATE_MODE:
2230  case ORATE_MODE:
2231  case GRATE_MODE:
2232  case WRATE_MODE:
2233  case LRATE_MODE:
2234  // Diag
2235  fill(bmat.begin(), bmat.end(), 0.0);
2236  for (USI i = 0; i < nc; i++) {
2237  bmat[0] += dQdXpW[(i + 1) * ncol];
2238  bmat[(i + 1) * ncol + i + 1] = 1;
2239  }
2240  for (USI i = 0; i < bsize; i++) {
2241  myLS.diagVal[wId * bsize + i] += bmat[i];
2242  }
2243 
2244  // OffDiag
2245  bmat = dQdXpB;
2246  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2247  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2248  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(),
2249  &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1, bmat.data());
2250  fill(bmat2.begin(), bmat2.end(), 0.0);
2251  for (USI i = 0; i < nc; i++) {
2252  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
2253  }
2254  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
2255  myLS.colId[wId].push_back(n);
2256  break;
2257 
2258  case BHP_MODE:
2259  // Diag
2260  fill(bmat.begin(), bmat.end(), 0.0);
2261  for (USI i = 0; i < nc + 1; i++) {
2262  bmat[i * ncol + i] = 1;
2263  }
2264  // Add
2265  for (USI i = 0; i < bsize; i++) {
2266  myLS.diagVal[wId * bsize + i] += bmat[i];
2267  }
2268  // OffDiag
2269  fill(bmat.begin(), bmat.end(), 0.0);
2270  // Insert
2271  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
2272  myLS.colId[wId].push_back(n);
2273  // Solution
2274  // myLS.u[wId * ncol] = opt.maxBHP - BHP;
2275  break;
2276 
2277  default:
2278  OCP_ABORT("Wrong Well Opt mode!");
2279  break;
2280  }
2281  }
2282  assert(myLS.val[wId].size() == numPerf * bsize);
2283  // Well self
2284  myLS.colId[wId].push_back(wId);
2285  myLS.diagPtr[wId] = numPerf;
2286  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
2287  myLS.diagVal.data() + wId * bsize + bsize);
2288 }
2289 
2290 
2291 void Well::AssembleMatPROD_FIM_new(const Bulk& myBulk, LinearSystem& myLS,
2292  const OCP_DBL& dt) const
2293 {
2294  OCP_FUNCNAME;
2295 
2296  const OCP_USI wId = myLS.dim;
2297  // important !
2298  myLS.dim++;
2299 
2300  const USI np = myBulk.numPhase;
2301  const USI nc = myBulk.numCom;
2302  const USI nch = nc - 1;
2303  const USI ncol = nc + 1;
2304  const USI ncol2 = np * nc + np;
2305  const USI bsize = ncol * ncol;
2306  const USI bsize2 = ncol * ncol2;
2307 
2308  OCP_DBL xij, xi, mu, muP, xiP;
2309  OCP_DBL dP;
2310  OCP_DBL transIJ;
2311  OCP_DBL tmp;
2312 
2313  OCP_USI n_np_j;
2314 
2315  vector<OCP_DBL> bmat(bsize, 0);
2316  vector<OCP_DBL> bmat2(bsize, 0);
2317  vector<OCP_DBL> dQdXpB(bsize, 0);
2318  vector<OCP_DBL> dQdXpW(bsize, 0);
2319  vector<OCP_DBL> dQdXsB(bsize2, 0);
2320  vector<bool> phaseExistB(np, false);
2321  vector<USI> pEnumComB(np, 0);
2322  USI ncolB;
2323 
2324 
2325  // Set Prod Weight
2326  if (opt.optMode != BHP_MODE) CalProdWeight(myBulk);
2327 
2328  for (USI p = 0; p < numPerf; p++) {
2329  OCP_USI n = perf[p].location;
2330  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2331  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2332  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2333 
2334  const USI npB = myBulk.phaseNum[n] + 1; ncolB = npB;
2335  for (USI j = 0; j < np; j++) {
2336  phaseExistB[j] = myBulk.phaseExist[n * np + j];
2337  pEnumComB[j] = myBulk.pEnumCom[n * np + j];
2338  ncolB += pEnumComB[j];
2339  }
2340 
2341 
2342  USI jxB = npB;
2343  for (USI j = 0; j < np; j++) {
2344 
2345  if (!phaseExistB[j]) {
2346  jxB += pEnumComB[j];
2347  continue;
2348  }
2349 
2350  n_np_j = n * np + j;
2351  dP = myBulk.Pj[n_np_j] - BHP - dG[p];
2352  xi = myBulk.xi[n_np_j];
2353  mu = myBulk.mu[n_np_j];
2354  muP = myBulk.muP[n_np_j];
2355  xiP = myBulk.xiP[n_np_j];
2356 
2357 
2358  for (USI i = 0; i < nc; i++) {
2359  xij = myBulk.xij[n_np_j * nc + i];
2360  // dQ / dP
2361  transIJ = perf[p].transj[j] * xi * xij;
2362  dQdXpB[(i + 1) * ncol] +=
2363  transIJ * (1 - dP * muP / mu) + dP * perf[p].transj[j] * xij * xiP;
2364  dQdXpW[(i + 1) * ncol] += -transIJ;
2365 
2366  // dQ / dS
2367  USI j1B = 0;
2368  for (USI j1 = 0; j1 < np; j1++) {
2369  if (phaseExistB[j1]) {
2370  tmp = CONV1 * perf[p].WI * perf[p].multiplier * dP / mu * xi * xij *
2371  myBulk.dKr_dS[n_np_j * np + j1];
2372  // capillary pressure
2373  tmp += transIJ * myBulk.dPcj_dS[n_np_j * np + j1];
2374  dQdXsB[(i + 1) * ncolB + j1B] += tmp;
2375  j1B++;
2376  }
2377  }
2378 
2379  for (USI k = 0; k < pEnumComB[j]; k++) {
2380  tmp = dP * perf[p].transj[j] * xij *
2381  (myBulk.xix[n_np_j * nc + k] -
2382  xi / mu * myBulk.mux[n_np_j * nc + k]);
2383  dQdXsB[(i + 1) * ncolB + jxB + k] += tmp;
2384  }
2385  // WARNING !!!
2386  if (i < pEnumComB[j])
2387  dQdXsB[(i + 1) * ncolB + jxB + i] += perf[p].transj[j] * xi * dP;;
2388  }
2389  jxB += pEnumComB[j];
2390  }
2391  // Bulk to Well
2392  bmat = dQdXpB;
2393  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2394  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2395  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1,
2396  bmat.data());
2397 
2398  Dscalar(bsize, dt, bmat.data());
2399  // Add
2400  USI ptr = myLS.diagPtr[n];
2401  for (USI i = 0; i < bsize; i++) {
2402  myLS.val[n][ptr * bsize + i] += bmat[i];
2403  }
2404  // Insert
2405  bmat = dQdXpW;
2406  Dscalar(bsize, dt, bmat.data());
2407  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
2408  myLS.colId[n].push_back(wId);
2409 
2410  // Well
2411  switch (opt.optMode) {
2412  case RATE_MODE:
2413  case ORATE_MODE:
2414  case GRATE_MODE:
2415  case WRATE_MODE:
2416  case LRATE_MODE:
2417  // Diag
2418  fill(bmat.begin(), bmat.end(), 0.0);
2419  for (USI i = 0; i < nc; i++) {
2420  bmat[0] += dQdXpW[(i + 1) * ncol] * prodWeight[i];
2421  bmat[(i + 1) * ncol + i + 1] = 1;
2422  }
2423  for (USI i = 0; i < bsize; i++) {
2424  myLS.diagVal[wId * bsize + i] += bmat[i];
2425  }
2426 
2427  // OffDiag
2428  bmat = dQdXpB;
2429  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2430  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2431  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(),
2432  &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1, bmat.data());
2433  fill(bmat2.begin(), bmat2.end(), 0.0);
2434  for (USI i = 0; i < nc; i++) {
2435  Daxpy(ncol, prodWeight[i], bmat.data() + (i + 1) * ncol,
2436  bmat2.data());
2437  }
2438  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
2439  myLS.colId[wId].push_back(n);
2440  break;
2441 
2442  case BHP_MODE:
2443  // Diag
2444  fill(bmat.begin(), bmat.end(), 0.0);
2445  for (USI i = 0; i < nc + 1; i++) {
2446  bmat[i * ncol + i] = 1;
2447  }
2448  // Add
2449  for (USI i = 0; i < bsize; i++) {
2450  myLS.diagVal[wId * bsize + i] += bmat[i];
2451  }
2452  // OffDiag
2453  fill(bmat.begin(), bmat.end(), 0.0);
2454  // Insert
2455  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
2456  myLS.colId[wId].push_back(n);
2457  // Solution
2458  // myLS.u[wId * ncol] = opt.minBHP - BHP;
2459  break;
2460 
2461  default:
2462  OCP_ABORT("Wrong Well Opt mode!");
2463  break;
2464  }
2465  }
2466  assert(myLS.val[wId].size() == numPerf * bsize);
2467  // Well self
2468  myLS.colId[wId].push_back(wId);
2469  myLS.diagPtr[wId] = numPerf;
2470  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
2471  myLS.diagVal.data() + wId * bsize + bsize);
2472 }
2473 
2474 
2475 void Well::AssembleMatINJ_FIM_new_n(const Bulk& myBulk, LinearSystem& myLS,
2476  const OCP_DBL& dt) const
2477 {
2478  OCP_FUNCNAME;
2479 
2480  const OCP_USI wId = myLS.dim;
2481  // important !
2482  myLS.dim++;
2483 
2484  const USI np = myBulk.numPhase;
2485  const USI nc = myBulk.numCom;
2486  const USI nch = nc - 1;
2487  const USI ncol = nc + 1;
2488  const USI ncol2 = np * nc + np;
2489  const USI bsize = ncol * ncol;
2490  const USI bsize2 = ncol * ncol2;
2491 
2492  OCP_DBL mu, muP;
2493  OCP_DBL dP;
2494  OCP_DBL transIJ;
2495 
2496  OCP_USI n_np_j;
2497 
2498  vector<OCP_DBL> bmat(bsize, 0);
2499  vector<OCP_DBL> bmat2(bsize, 0);
2500  vector<OCP_DBL> dQdXpB(bsize, 0);
2501  vector<OCP_DBL> dQdXpW(bsize, 0);
2502  vector<OCP_DBL> dQdXsB(bsize2, 0);
2503  vector<bool> phaseExistB(np, false);
2504  vector<USI> pEnumComB(np, 0);
2505  USI ncolB;
2506 
2507 
2508  for (USI p = 0; p < numPerf; p++) {
2509  const OCP_USI n = perf[p].location;
2510  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2511  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2512  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2513 
2514  dP = myBulk.P[n] - BHP - dG[p];
2515 
2516  const USI npB = myBulk.phaseNum[n] + 1; ncolB = npB;
2517  for (USI j = 0; j < np; j++) {
2518  phaseExistB[j] = myBulk.phaseExist[n * np + j];
2519  pEnumComB[j] = myBulk.pEnumCom[n * np + j];
2520  ncolB += pEnumComB[j];
2521  }
2522 
2523 
2524  USI jxB = npB;
2525  for (USI j = 0; j < np; j++) {
2526 
2527  if (!phaseExistB[j]) {
2528  jxB += pEnumComB[j];
2529  continue;
2530  }
2531 
2532  n_np_j = n * np + j;
2533  mu = myBulk.mu[n_np_j];
2534  muP = myBulk.muP[n_np_j];
2535 
2536  for (USI i = 0; i < nc; i++) {
2537  // dQ / dP
2538  transIJ = perf[p].transj[j] * perf[p].xi * opt.zi[i];
2539  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
2540  dQdXpW[(i + 1) * ncol] += -transIJ;
2541 
2542  // dQ / dS
2543  USI j1B = 0;
2544  for (USI j1 = 0; j1 < np; j1++) {
2545  if (phaseExistB[j1]) {
2546  dQdXsB[(i + 1) * ncolB + j1B] +=
2547  CONV1 * perf[p].WI * perf[p].multiplier * perf[p].xi *
2548  opt.zi[i] * myBulk.dKr_dS[n_np_j * np + j1] * dP / mu;
2549  j1B++;
2550  }
2551  }
2552 
2553  // dQ / dxij
2554  for (USI k = 0; k < pEnumComB[j]; k++) {
2555  dQdXsB[(i + 1) * ncolB + jxB + k] +=
2556  -transIJ * dP / mu * myBulk.mux[n_np_j * nc + k];
2557  }
2558  }
2559  jxB += pEnumComB[j];
2560  }
2561 
2562  // Assemble rhs
2563  // Well
2564  if (npB > 2) {
2565  DaAxpby(ncol, ncolB, -1.0, dQdXsB.data(),
2566  &myBulk.res_n[myBulk.resIndex[n]], 1.0, &myLS.b[n * ncol]);
2567  }
2568 
2569 
2570  // Bulk to Well
2571  bmat = dQdXpB;
2572  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2573  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2574  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1,
2575  bmat.data());
2576  Dscalar(bsize, dt, bmat.data());
2577  // Add
2578  USI ptr = myLS.diagPtr[n];
2579  for (USI i = 0; i < bsize; i++) {
2580  myLS.val[n][ptr * bsize + i] += bmat[i];
2581  }
2582  // Insert
2583  bmat = dQdXpW;
2584  Dscalar(bsize, dt, bmat.data());
2585  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
2586  myLS.colId[n].push_back(wId);
2587 
2588  // Well
2589  switch (opt.optMode) {
2590  case RATE_MODE:
2591  case ORATE_MODE:
2592  case GRATE_MODE:
2593  case WRATE_MODE:
2594  case LRATE_MODE:
2595  // Diag
2596  fill(bmat.begin(), bmat.end(), 0.0);
2597  for (USI i = 0; i < nc; i++) {
2598  bmat[0] += dQdXpW[(i + 1) * ncol];
2599  bmat[(i + 1) * ncol + i + 1] = 1;
2600  }
2601  for (USI i = 0; i < bsize; i++) {
2602  myLS.diagVal[wId * bsize + i] += bmat[i];
2603  }
2604 
2605  // OffDiag
2606  bmat = dQdXpB;
2607  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2608  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2609  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(),
2610  &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1, bmat.data());
2611  fill(bmat2.begin(), bmat2.end(), 0.0);
2612  for (USI i = 0; i < nc; i++) {
2613  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
2614  }
2615  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
2616  myLS.colId[wId].push_back(n);
2617  break;
2618 
2619  case BHP_MODE:
2620  // Diag
2621  fill(bmat.begin(), bmat.end(), 0.0);
2622  for (USI i = 0; i < nc + 1; i++) {
2623  bmat[i * ncol + i] = 1;
2624  }
2625  // Add
2626  for (USI i = 0; i < bsize; i++) {
2627  myLS.diagVal[wId * bsize + i] += bmat[i];
2628  }
2629  // OffDiag
2630  fill(bmat.begin(), bmat.end(), 0.0);
2631  // Insert
2632  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
2633  myLS.colId[wId].push_back(n);
2634  // Solution
2635  // myLS.u[wId * ncol] = opt.maxBHP - BHP;
2636  break;
2637 
2638  default:
2639  OCP_ABORT("Wrong Well Opt mode!");
2640  break;
2641  }
2642  }
2643  assert(myLS.val[wId].size() == numPerf * bsize);
2644  // Well self
2645  myLS.colId[wId].push_back(wId);
2646  myLS.diagPtr[wId] = numPerf;
2647  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
2648  myLS.diagVal.data() + wId * bsize + bsize);
2649 }
2650 
2651 
2652 void Well::AssembleMatPROD_FIM_new_n(const Bulk& myBulk, LinearSystem& myLS,
2653  const OCP_DBL& dt) const
2654 {
2655  OCP_FUNCNAME;
2656 
2657  const OCP_USI wId = myLS.dim;
2658  // important !
2659  myLS.dim++;
2660 
2661  const USI np = myBulk.numPhase;
2662  const USI nc = myBulk.numCom;
2663  const USI nch = nc - 1;
2664  const USI ncol = nc + 1;
2665  const USI ncol2 = np * nc + np;
2666  const USI bsize = ncol * ncol;
2667  const USI bsize2 = ncol * ncol2;
2668 
2669  OCP_DBL xij, xi, mu, muP, xiP;
2670  OCP_DBL dP;
2671  OCP_DBL transIJ;
2672  OCP_DBL tmp;
2673 
2674  OCP_USI n_np_j;
2675 
2676  vector<OCP_DBL> bmat(bsize, 0);
2677  vector<OCP_DBL> bmat2(bsize, 0);
2678  vector<OCP_DBL> dQdXpB(bsize, 0);
2679  vector<OCP_DBL> dQdXpW(bsize, 0);
2680  vector<OCP_DBL> dQdXsB(bsize2, 0);
2681  vector<bool> phaseExistB(np, false);
2682  vector<USI> pEnumComB(np, 0);
2683  USI ncolB;
2684 
2685 
2686  // Set Prod Weight
2687  if (opt.optMode != BHP_MODE) CalProdWeight(myBulk);
2688 
2689  for (USI p = 0; p < numPerf; p++) {
2690  OCP_USI n = perf[p].location;
2691  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2692  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2693  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2694 
2695  const USI npB = myBulk.phaseNum[n] + 1; ncolB = npB;
2696  for (USI j = 0; j < np; j++) {
2697  phaseExistB[j] = myBulk.phaseExist[n * np + j];
2698  pEnumComB[j] = myBulk.pEnumCom[n * np + j];
2699  ncolB += pEnumComB[j];
2700  }
2701 
2702 
2703  USI jxB = npB;
2704  for (USI j = 0; j < np; j++) {
2705 
2706  if (!phaseExistB[j]) {
2707  jxB += pEnumComB[j];
2708  continue;
2709  }
2710 
2711  n_np_j = n * np + j;
2712  dP = myBulk.Pj[n_np_j] - BHP - dG[p];
2713  xi = myBulk.xi[n_np_j];
2714  mu = myBulk.mu[n_np_j];
2715  muP = myBulk.muP[n_np_j];
2716  xiP = myBulk.xiP[n_np_j];
2717 
2718 
2719  for (USI i = 0; i < nc; i++) {
2720  xij = myBulk.xij[n_np_j * nc + i];
2721  // dQ / dP
2722  transIJ = perf[p].transj[j] * xi * xij;
2723  dQdXpB[(i + 1) * ncol] +=
2724  transIJ * (1 - dP * muP / mu) + dP * perf[p].transj[j] * xij * xiP;
2725  dQdXpW[(i + 1) * ncol] += -transIJ;
2726 
2727  // dQ / dS
2728  USI j1B = 0;
2729  for (USI j1 = 0; j1 < np; j1++) {
2730  if (phaseExistB[j1]) {
2731  tmp = CONV1 * perf[p].WI * perf[p].multiplier * dP / mu * xi * xij *
2732  myBulk.dKr_dS[n_np_j * np + j1];
2733  // capillary pressure
2734  tmp += transIJ * myBulk.dPcj_dS[n_np_j * np + j1];
2735  dQdXsB[(i + 1) * ncolB + j1B] += tmp;
2736  j1B++;
2737  }
2738  }
2739 
2740  for (USI k = 0; k < pEnumComB[j]; k++) {
2741  tmp = dP * perf[p].transj[j] * xij *
2742  (myBulk.xix[n_np_j * nc + k] -
2743  xi / mu * myBulk.mux[n_np_j * nc + k]);
2744  tmp -= transIJ * dP / myBulk.nj[n_np_j];
2745  dQdXsB[(i + 1) * ncolB + jxB + k] += tmp;
2746  }
2747  // WARNING !!!
2748  if (i < pEnumComB[j])
2749  dQdXsB[(i + 1) * ncolB + jxB + i] += perf[p].transj[j] * xi * dP / myBulk.nj[n_np_j];
2750  }
2751  jxB += pEnumComB[j];
2752  }
2753 
2754  // Assemble rhs
2755  // Well
2756  if (npB > 2) {
2757  DaAxpby(ncol, ncolB, -1.0, dQdXsB.data(),
2758  &myBulk.res_n[myBulk.resIndex[n]], 1.0, &myLS.b[n * ncol]);
2759  }
2760 
2761 
2762  // Bulk to Well
2763  bmat = dQdXpB;
2764  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2765  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2766  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1,
2767  bmat.data());
2768 
2769  Dscalar(bsize, dt, bmat.data());
2770  // Add
2771  USI ptr = myLS.diagPtr[n];
2772  for (USI i = 0; i < bsize; i++) {
2773  myLS.val[n][ptr * bsize + i] += bmat[i];
2774  }
2775  // Insert
2776  bmat = dQdXpW;
2777  Dscalar(bsize, dt, bmat.data());
2778  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
2779  myLS.colId[n].push_back(wId);
2780 
2781  // Well
2782  switch (opt.optMode) {
2783  case RATE_MODE:
2784  case ORATE_MODE:
2785  case GRATE_MODE:
2786  case WRATE_MODE:
2787  case LRATE_MODE:
2788  // Diag
2789  fill(bmat.begin(), bmat.end(), 0.0);
2790  for (USI i = 0; i < nc; i++) {
2791  bmat[0] += dQdXpW[(i + 1) * ncol] * prodWeight[i];
2792  bmat[(i + 1) * ncol + i + 1] = 1;
2793  }
2794  for (USI i = 0; i < bsize; i++) {
2795  myLS.diagVal[wId * bsize + i] += bmat[i];
2796  }
2797 
2798  // OffDiag
2799  bmat = dQdXpB;
2800  //myDABpCp2(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[n * bsize2], bmat.data(), &flagB[0], nc - 1);
2801  //myDABpC(ncol, ncolB, ncol, dQdXsB.data(), &myBulk.dSec_dPri[myBulk.dSdPindex[n]], bmat.data());
2802  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(),
2803  &myBulk.dSec_dPri[myBulk.dSdPindex[n]], 1, bmat.data());
2804  fill(bmat2.begin(), bmat2.end(), 0.0);
2805  for (USI i = 0; i < nc; i++) {
2806  Daxpy(ncol, prodWeight[i], bmat.data() + (i + 1) * ncol,
2807  bmat2.data());
2808  }
2809  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
2810  myLS.colId[wId].push_back(n);
2811  break;
2812 
2813  case BHP_MODE:
2814  // Diag
2815  fill(bmat.begin(), bmat.end(), 0.0);
2816  for (USI i = 0; i < nc + 1; i++) {
2817  bmat[i * ncol + i] = 1;
2818  }
2819  // Add
2820  for (USI i = 0; i < bsize; i++) {
2821  myLS.diagVal[wId * bsize + i] += bmat[i];
2822  }
2823  // OffDiag
2824  fill(bmat.begin(), bmat.end(), 0.0);
2825  // Insert
2826  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
2827  myLS.colId[wId].push_back(n);
2828  // Solution
2829  // myLS.u[wId * ncol] = opt.minBHP - BHP;
2830  break;
2831 
2832  default:
2833  OCP_ABORT("Wrong Well Opt mode!");
2834  break;
2835  }
2836  }
2837  assert(myLS.val[wId].size() == numPerf * bsize);
2838  // Well self
2839  myLS.colId[wId].push_back(wId);
2840  myLS.diagPtr[wId] = numPerf;
2841  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
2842  myLS.diagVal.data() + wId * bsize + bsize);
2843 }
2844 
2845 
2846 void Well::AssembleMatINJ_AIMt(const Bulk& myBulk, LinearSystem& myLS,
2847  const OCP_DBL& dt) const
2848 {
2849  OCP_FUNCNAME;
2850 
2851  const OCP_USI wId = myLS.dim;
2852  // important !
2853  myLS.dim++;
2854 
2855  const USI np = myBulk.numPhase;
2856  const USI nc = myBulk.numCom;
2857  const USI ncol = nc + 1;
2858  const USI ncol2 = np * nc + np;
2859  const USI bsize = ncol * ncol;
2860  const USI bsize2 = ncol * ncol2;
2861 
2862  OCP_DBL mu, muP;
2863  OCP_DBL dP;
2864  OCP_DBL transIJ;
2865 
2866  OCP_USI n_np_j;
2867  OCP_USI pIde, pIde_np_j;
2868 
2869  vector<OCP_DBL> bmat(bsize, 0);
2870  vector<OCP_DBL> bmat2(bsize, 0);
2871  vector<OCP_DBL> dQdXpB(bsize, 0);
2872  vector<OCP_DBL> dQdXpW(bsize, 0);
2873  vector<OCP_DBL> dQdXsB(bsize2, 0);
2874 
2875  for (USI p = 0; p < numPerf; p++) {
2876  OCP_USI n = perf[p].location;
2877  pIde = myBulk.map_Bulk2FIM[n];
2878 
2879  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2880  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2881  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2882 
2883  dP = myBulk.P[n] - BHP - dG[p];
2884 
2885  for (USI j = 0; j < np; j++) {
2886  n_np_j = n * np + j;
2887  if (!myBulk.phaseExist[n_np_j]) continue;
2888 
2889  pIde_np_j = pIde * np + j;
2890  mu = myBulk.mu[n_np_j];
2891  muP = myBulk.muP[pIde_np_j];
2892 
2893  for (USI i = 0; i < nc; i++) {
2894  // dQ / dP
2895  transIJ = perf[p].transj[j] * perf[p].xi * opt.zi[i];
2896  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
2897  dQdXpW[(i + 1) * ncol] += -transIJ;
2898 
2899  // dQ / dS
2900  for (USI k = 0; k < np; k++) {
2901  dQdXsB[(i + 1) * ncol2 + k] +=
2902  CONV1 * perf[p].WI * perf[p].multiplier * perf[p].xi *
2903  opt.zi[i] * myBulk.dKr_dS[pIde * np * np + j * np + k] * dP / mu;
2904  }
2905  // dQ / dxij
2906  for (USI k = 0; k < nc; k++) {
2907  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] +=
2908  -transIJ * dP / mu * myBulk.mux[pIde_np_j * nc + k];
2909  }
2910  }
2911  }
2912  // Bulk to Well
2913  bmat = dQdXpB;
2914  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[pIde * bsize2], 1,
2915  bmat.data());
2916  Dscalar(bsize, dt, bmat.data());
2917  // Add
2918  USI ptr = myLS.diagPtr[pIde];
2919  for (USI i = 0; i < bsize; i++) {
2920  myLS.val[pIde][ptr * bsize + i] += bmat[i];
2921  }
2922  // Insert
2923  bmat = dQdXpW;
2924  Dscalar(bsize, dt, bmat.data());
2925  myLS.val[pIde].insert(myLS.val[pIde].end(), bmat.begin(), bmat.end());
2926  myLS.colId[pIde].push_back(wId);
2927 
2928  // Well
2929  switch (opt.optMode) {
2930  case RATE_MODE:
2931  case ORATE_MODE:
2932  case GRATE_MODE:
2933  case WRATE_MODE:
2934  case LRATE_MODE:
2935  // Diag
2936  fill(bmat.begin(), bmat.end(), 0.0);
2937  for (USI i = 0; i < nc; i++) {
2938  bmat[0] += dQdXpW[(i + 1) * ncol];
2939  bmat[(i + 1) * ncol + i + 1] = 1;
2940  }
2941  for (USI i = 0; i < bsize; i++) {
2942  myLS.diagVal[wId * bsize + i] += bmat[i];
2943  }
2944 
2945  // OffDiag
2946  bmat = dQdXpB;
2947  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(),
2948  &myBulk.dSec_dPri[pIde * bsize2], 1, bmat.data());
2949  fill(bmat2.begin(), bmat2.end(), 0.0);
2950  for (USI i = 0; i < nc; i++) {
2951  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
2952  }
2953  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
2954  myLS.colId[wId].push_back(pIde);
2955  break;
2956 
2957  case BHP_MODE:
2958  // Diag
2959  fill(bmat.begin(), bmat.end(), 0.0);
2960  for (USI i = 0; i < nc + 1; i++) {
2961  bmat[i * ncol + i] = 1;
2962  }
2963  // Add
2964  for (USI i = 0; i < bsize; i++) {
2965  myLS.diagVal[wId * bsize + i] += bmat[i];
2966  }
2967  // OffDiag
2968  fill(bmat.begin(), bmat.end(), 0.0);
2969  // Insert
2970  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
2971  myLS.colId[wId].push_back(pIde);
2972  // Solution
2973  // myLS.u[wId * ncol] = opt.maxBHP - BHP;
2974  break;
2975 
2976  default:
2977  OCP_ABORT("Wrong Well Opt mode!");
2978  break;
2979  }
2980  }
2981  assert(myLS.val[wId].size() == numPerf * bsize);
2982  // Well self
2983  myLS.colId[wId].push_back(wId);
2984  myLS.diagPtr[wId] = numPerf;
2985  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
2986  myLS.diagVal.data() + wId * bsize + bsize);
2987 }
2988 
2989 
2990 void Well::AssembleMatPROD_AIMt(const Bulk& myBulk, LinearSystem& myLS,
2991  const OCP_DBL& dt) const
2992 {
2993  OCP_FUNCNAME;
2994 
2995  const OCP_USI wId = myLS.dim;
2996  // important !
2997  myLS.dim++;
2998 
2999  const USI np = myBulk.numPhase;
3000  const USI nc = myBulk.numCom;
3001  const USI ncol = nc + 1;
3002  const USI ncol2 = np * nc + np;
3003  const USI bsize = ncol * ncol;
3004  const USI bsize2 = ncol * ncol2;
3005 
3006  OCP_DBL xij, xi, mu, muP, xiP;
3007  OCP_DBL dP;
3008  OCP_DBL transIJ;
3009  OCP_DBL tmp;
3010 
3011  OCP_USI n_np_j;
3012  OCP_USI pIde, pIde_np_j;
3013 
3014  vector<OCP_DBL> bmat(bsize, 0);
3015  vector<OCP_DBL> bmat2(bsize, 0);
3016  vector<OCP_DBL> dQdXpB(bsize, 0);
3017  vector<OCP_DBL> dQdXpW(bsize, 0);
3018  vector<OCP_DBL> dQdXsB(bsize2, 0);
3019 
3020  // Set Prod Weight
3021  if (opt.optMode != BHP_MODE) CalProdWeight(myBulk);
3022 
3023  for (USI p = 0; p < numPerf; p++) {
3024  OCP_USI n = perf[p].location;
3025  pIde = myBulk.map_Bulk2FIM[n];
3026  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
3027  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
3028  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
3029 
3030  for (USI j = 0; j < np; j++) {
3031  n_np_j = n * np + j;
3032  if (!myBulk.phaseExist[n_np_j]) continue;
3033 
3034  pIde_np_j = pIde * np + j;
3035  dP = myBulk.Pj[n_np_j] - BHP - dG[p];
3036  xi = myBulk.xi[n_np_j];
3037  mu = myBulk.mu[n_np_j];
3038  muP = myBulk.muP[pIde_np_j];
3039  xiP = myBulk.xiP[pIde_np_j];
3040 
3041  for (USI i = 0; i < nc; i++) {
3042  xij = myBulk.xij[n_np_j * nc + i];
3043  // dQ / dP
3044  transIJ = perf[p].transj[j] * xi * xij;
3045  dQdXpB[(i + 1) * ncol] +=
3046  transIJ * (1 - dP * muP / mu) + dP * perf[p].transj[j] * xij * xiP;
3047  dQdXpW[(i + 1) * ncol] += -transIJ;
3048 
3049  // dQ / dS
3050  for (USI k = 0; k < np; k++) {
3051  tmp = CONV1 * perf[p].WI * perf[p].multiplier * dP / mu * xi * xij *
3052  myBulk.dKr_dS[pIde * np * np + j * np + k];
3053  // capillary pressure
3054  tmp += transIJ * myBulk.dPcj_dS[pIde * np * np + j * np + k];
3055  dQdXsB[(i + 1) * ncol2 + k] += tmp;
3056  }
3057  // dQ / dCij
3058  for (USI k = 0; k < nc; k++) {
3059  tmp = dP * perf[p].transj[j] * xij *
3060  (myBulk.xix[pIde_np_j * nc + k] -
3061  xi / mu * myBulk.mux[pIde_np_j * nc + k]);
3062  if (k == i) {
3063  tmp += perf[p].transj[j] * xi * dP;
3064  }
3065  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] += tmp;
3066  }
3067  }
3068  }
3069  // Bulk to Well
3070  bmat = dQdXpB;
3071  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[pIde * bsize2], 1,
3072  bmat.data());
3073 
3074  // for (USI i = 0; i < bsize; i++) {
3075  // cout << fixed << setprecision(6) << bmat[i] << endl;
3076  //}
3077 
3078  Dscalar(bsize, dt, bmat.data());
3079  // Add
3080  USI ptr = myLS.diagPtr[pIde];
3081  for (USI i = 0; i < bsize; i++) {
3082  myLS.val[pIde][ptr * bsize + i] += bmat[i];
3083  }
3084  // Insert
3085  bmat = dQdXpW;
3086  Dscalar(bsize, dt, bmat.data());
3087  myLS.val[pIde].insert(myLS.val[pIde].end(), bmat.begin(), bmat.end());
3088  myLS.colId[pIde].push_back(wId);
3089 
3090  // Well
3091  switch (opt.optMode) {
3092  case RATE_MODE:
3093  case ORATE_MODE:
3094  case GRATE_MODE:
3095  case WRATE_MODE:
3096  case LRATE_MODE:
3097  // Diag
3098  fill(bmat.begin(), bmat.end(), 0.0);
3099  for (USI i = 0; i < nc; i++) {
3100  bmat[0] += dQdXpW[(i + 1) * ncol] * prodWeight[i];
3101  bmat[(i + 1) * ncol + i + 1] = 1;
3102  }
3103  for (USI i = 0; i < bsize; i++) {
3104  myLS.diagVal[wId * bsize + i] += bmat[i];
3105  }
3106 
3107  // OffDiag
3108  bmat = dQdXpB;
3109  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(),
3110  &myBulk.dSec_dPri[pIde * bsize2], 1, bmat.data());
3111  fill(bmat2.begin(), bmat2.end(), 0.0);
3112  for (USI i = 0; i < nc; i++) {
3113  Daxpy(ncol, prodWeight[i], bmat.data() + (i + 1) * ncol,
3114  bmat2.data());
3115  }
3116  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
3117  myLS.colId[wId].push_back(pIde);
3118  break;
3119 
3120  case BHP_MODE:
3121  // Diag
3122  fill(bmat.begin(), bmat.end(), 0.0);
3123  for (USI i = 0; i < nc + 1; i++) {
3124  bmat[i * ncol + i] = 1;
3125  }
3126  // Add
3127  for (USI i = 0; i < bsize; i++) {
3128  myLS.diagVal[wId * bsize + i] += bmat[i];
3129  }
3130  // OffDiag
3131  fill(bmat.begin(), bmat.end(), 0.0);
3132  // Insert
3133  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
3134  myLS.colId[wId].push_back(pIde);
3135  // Solution
3136  // myLS.u[wId * ncol] = opt.minBHP - BHP;
3137  break;
3138 
3139  default:
3140  OCP_ABORT("Wrong Well Opt mode!");
3141  break;
3142  }
3143  }
3144  assert(myLS.val[wId].size() == numPerf * bsize);
3145  // Well self
3146  myLS.colId[wId].push_back(wId);
3147  myLS.diagPtr[wId] = numPerf;
3148  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
3149  myLS.diagVal.data() + wId * bsize + bsize);
3150 
3151  // test
3152  // for (USI i = 0; i < bsize; i++) {
3153  // cout << myLS.diagVal[wId * bsize + i] << endl;
3154  //}
3155 }
3156 
3157 void Well::CalResAIMt(ResFIM& resFIM, const Bulk& myBulk, const OCP_DBL& dt,
3158  const OCP_USI& wId, const vector<Well>& allWell) const
3159 {
3160  OCP_FUNCNAME;
3161 
3162  // Well to Bulk
3163  const USI nc = myBulk.numCom;
3164  const USI len = nc + 1;
3165  OCP_USI k, bIde;
3166 
3167  for (USI p = 0; p < numPerf; p++) {
3168  k = perf[p].location;
3169  bIde = myBulk.map_Bulk2FIM[k];
3170  for (USI i = 0; i < nc; i++) {
3171  resFIM.res[bIde * len + 1 + i] += perf[p].qi_lbmol[i] * dt;
3172  }
3173  }
3174 
3175  OCP_USI bId = (myBulk.numFIMBulk + wId) * len;
3176  // Well Self
3177  if (opt.type == INJ) {
3178  // Injection
3179  switch (opt.optMode) {
3180  case BHP_MODE:
3181  BHP = opt.maxBHP;
3182  resFIM.res[bId] = BHP - opt.maxBHP;
3183  break;
3184  case RATE_MODE:
3185  case ORATE_MODE:
3186  case GRATE_MODE:
3187  case WRATE_MODE:
3188  case LRATE_MODE:
3189  resFIM.res[bId] = opt.maxRate;
3190  for (USI i = 0; i < nc; i++) {
3191  resFIM.res[bId] += qi_lbmol[i];
3192  }
3193  if (opt.reInj) {
3194  for (auto& w : opt.connWell) {
3195  OCP_DBL tmp = 0;
3196  for (USI i = 0; i < nc; i++) {
3197  tmp += allWell[w].qi_lbmol[i];
3198  // tmp += opt.factor * allWell[w].qi_lbmol[i];
3199  // resFIM.res[bId] += opt.factor * allWell[w].qi_lbmol[i];
3200  }
3201  tmp *= opt.factor;
3202  resFIM.res[bId] += tmp;
3203  // cout << "Temp(INJ): " << tmp / opt.xiINJ / 1000 << endl;
3204  }
3205  // cout << "Factor(res) " << opt.factor << endl;
3206  }
3207  // cout << name << " " << resFIM.res[bId] << " " << opt.maxRate << "
3208  // " << fabs(resFIM.res[bId] / opt.maxRate) << endl;
3209  resFIM.maxRelRes_v =
3210  max(resFIM.maxRelRes_v, fabs(resFIM.res[bId] / opt.maxRate));
3211  break;
3212  default:
3213  OCP_ABORT("Wrong well opt mode!");
3214  break;
3215  }
3216  }
3217  else {
3218  // Production
3219  switch (opt.optMode) {
3220  case BHP_MODE:
3221  BHP = opt.minBHP;
3222  resFIM.res[bId] = BHP - opt.minBHP;
3223  break;
3224  case RATE_MODE:
3225  case ORATE_MODE:
3226  case GRATE_MODE:
3227  case WRATE_MODE:
3228  case LRATE_MODE:
3229  resFIM.res[bId] = -opt.maxRate;
3230  for (USI i = 0; i < nc; i++) {
3231  resFIM.res[bId] += qi_lbmol[i] * prodWeight[i];
3232  }
3233  // cout << "Temp(Prod): " << tmp << endl;
3234  // cout << name << " " << resFIM.res[bId] << " " << opt.maxRate << "
3235  // " << fabs(resFIM.res[bId] / opt.maxRate) << endl;
3236  resFIM.maxRelRes_v =
3237  max(resFIM.maxRelRes_v, fabs(resFIM.res[bId] / opt.maxRate));
3238  break;
3239  default:
3240  OCP_ABORT("Wrong well opt mode!");
3241  break;
3242  }
3243  }
3244 }
3245 
3246 
3247 void Well::AssembleMatINJ_AIMs(const Bulk& myBulk, LinearSystem& myLS,
3248  const OCP_DBL& dt) const
3249 {
3250  OCP_FUNCNAME;
3251 
3252  const OCP_USI wId = myLS.dim;
3253  // important !
3254  myLS.dim++;
3255 
3256  const USI np = myBulk.numPhase;
3257  const USI nc = myBulk.numCom;
3258  const USI ncol = nc + 1;
3259  const USI ncol2 = np * nc + np;
3260  const USI bsize = ncol * ncol;
3261  const USI bsize2 = ncol * ncol2;
3262 
3263  OCP_DBL mu, muP;
3264  OCP_DBL dP;
3265  OCP_DBL transIJ;
3266 
3267  OCP_USI n_np_j;
3268 
3269  vector<OCP_DBL> bmat(bsize, 0);
3270  vector<OCP_DBL> bmat2(bsize, 0);
3271  vector<OCP_DBL> dQdXpB(bsize, 0);
3272  vector<OCP_DBL> dQdXpW(bsize, 0);
3273  vector<OCP_DBL> dQdXsB(bsize2, 0);
3274 
3275  for (USI p = 0; p < numPerf; p++) {
3276  OCP_USI n = perf[p].location;
3277  OCP_USI pId = myBulk.map_Bulk2FIM[n];
3278  OCP_ASSERT(pId < myBulk.numFIMBulk, "wrong FIM Bulk");
3279 
3280  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
3281  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
3282  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
3283 
3284  dP = myBulk.P[n] - BHP - dG[p];
3285 
3286  for (USI j = 0; j < np; j++) {
3287  n_np_j = n * np + j;
3288  if (!myBulk.phaseExist[n_np_j]) continue;
3289 
3290  mu = myBulk.mu[n_np_j];
3291  muP = myBulk.muP[n_np_j];
3292 
3293  for (USI i = 0; i < nc; i++) {
3294  // dQ / dP
3295  transIJ = perf[p].transj[j] * perf[p].xi * opt.zi[i];
3296  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
3297  dQdXpW[(i + 1) * ncol] += -transIJ;
3298 
3299  // dQ / dS
3300  for (USI k = 0; k < np; k++) {
3301  dQdXsB[(i + 1) * ncol2 + k] +=
3302  CONV1 * perf[p].WI * perf[p].multiplier * perf[p].xi *
3303  opt.zi[i] * myBulk.dKr_dS[pId * np * np + j * np + k] * dP / mu;
3304  }
3305  // dQ / dxij
3306  for (USI k = 0; k < nc; k++) {
3307  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] +=
3308  -transIJ * dP / mu * myBulk.mux[n_np_j * nc + k];
3309  }
3310  }
3311  }
3312  // Bulk to Well
3313  bmat = dQdXpB;
3314  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[pId * bsize2], 1,
3315  bmat.data());
3316  Dscalar(bsize, dt, bmat.data());
3317  // Add
3318  USI ptr = myLS.diagPtr[n];
3319  for (USI i = 0; i < bsize; i++) {
3320  myLS.val[n][ptr * bsize + i] += bmat[i];
3321  }
3322  // Insert
3323  bmat = dQdXpW;
3324  Dscalar(bsize, dt, bmat.data());
3325  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
3326  myLS.colId[n].push_back(wId);
3327 
3328  // Well
3329  switch (opt.optMode) {
3330  case RATE_MODE:
3331  case ORATE_MODE:
3332  case GRATE_MODE:
3333  case WRATE_MODE:
3334  case LRATE_MODE:
3335  // Diag
3336  fill(bmat.begin(), bmat.end(), 0.0);
3337  for (USI i = 0; i < nc; i++) {
3338  bmat[0] += dQdXpW[(i + 1) * ncol];
3339  bmat[(i + 1) * ncol + i + 1] = 1;
3340  }
3341  for (USI i = 0; i < bsize; i++) {
3342  myLS.diagVal[wId * bsize + i] += bmat[i];
3343  }
3344 
3345  // OffDiag
3346  bmat = dQdXpB;
3347  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(),
3348  &myBulk.dSec_dPri[pId * bsize2], 1, bmat.data());
3349  fill(bmat2.begin(), bmat2.end(), 0.0);
3350  for (USI i = 0; i < nc; i++) {
3351  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
3352  }
3353  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
3354  myLS.colId[wId].push_back(n);
3355  break;
3356 
3357  case BHP_MODE:
3358  // Diag
3359  fill(bmat.begin(), bmat.end(), 0.0);
3360  for (USI i = 0; i < nc + 1; i++) {
3361  bmat[i * ncol + i] = 1;
3362  }
3363  // Add
3364  for (USI i = 0; i < bsize; i++) {
3365  myLS.diagVal[wId * bsize + i] += bmat[i];
3366  }
3367  // OffDiag
3368  fill(bmat.begin(), bmat.end(), 0.0);
3369  // Insert
3370  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
3371  myLS.colId[wId].push_back(n);
3372  // Solution
3373  // myLS.u[wId * ncol] = opt.maxBHP - BHP;
3374  break;
3375 
3376  default:
3377  OCP_ABORT("Wrong Well Opt mode!");
3378  break;
3379  }
3380  }
3381  assert(myLS.val[wId].size() == numPerf * bsize);
3382  // Well self
3383  myLS.colId[wId].push_back(wId);
3384  myLS.diagPtr[wId] = numPerf;
3385  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
3386  myLS.diagVal.data() + wId * bsize + bsize);
3387 }
3388 
3389 
3390 void Well::AssembleMatPROD_AIMs(const Bulk& myBulk, LinearSystem& myLS,
3391  const OCP_DBL& dt) const
3392 {
3393  OCP_FUNCNAME;
3394 
3395  const OCP_USI wId = myLS.dim;
3396  // important !
3397  myLS.dim++;
3398 
3399  const USI np = myBulk.numPhase;
3400  const USI nc = myBulk.numCom;
3401  const USI ncol = nc + 1;
3402  const USI ncol2 = np * nc + np;
3403  const USI bsize = ncol * ncol;
3404  const USI bsize2 = ncol * ncol2;
3405 
3406  OCP_DBL xij, xi, mu, muP, xiP;
3407  OCP_DBL dP;
3408  OCP_DBL transIJ;
3409  OCP_DBL tmp;
3410 
3411  OCP_USI n_np_j;
3412 
3413  vector<OCP_DBL> bmat(bsize, 0);
3414  vector<OCP_DBL> bmat2(bsize, 0);
3415  vector<OCP_DBL> dQdXpB(bsize, 0);
3416  vector<OCP_DBL> dQdXpW(bsize, 0);
3417  vector<OCP_DBL> dQdXsB(bsize2, 0);
3418 
3419  // Set Prod Weight
3420  if (opt.optMode != BHP_MODE) CalProdWeight(myBulk);
3421 
3422  for (USI p = 0; p < numPerf; p++) {
3423  OCP_USI n = perf[p].location;
3424  OCP_USI pId = myBulk.map_Bulk2FIM[n];
3425  OCP_ASSERT(pId < myBulk.numFIMBulk, "wrong FIM Bulk");
3426 
3427  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
3428  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
3429  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
3430 
3431  for (USI j = 0; j < np; j++) {
3432  n_np_j = n * np + j;
3433  if (!myBulk.phaseExist[n_np_j]) continue;
3434 
3435  dP = myBulk.Pj[n_np_j] - BHP - dG[p];
3436  xi = myBulk.xi[n_np_j];
3437  mu = myBulk.mu[n_np_j];
3438  muP = myBulk.muP[n_np_j];
3439  xiP = myBulk.xiP[n_np_j];
3440 
3441  for (USI i = 0; i < nc; i++) {
3442  xij = myBulk.xij[n_np_j * nc + i];
3443  // dQ / dP
3444  transIJ = perf[p].transj[j] * xi * xij;
3445  dQdXpB[(i + 1) * ncol] +=
3446  transIJ * (1 - dP * muP / mu) + dP * perf[p].transj[j] * xij * xiP;
3447  dQdXpW[(i + 1) * ncol] += -transIJ;
3448 
3449  // if (dQdXpW[ncol] == 0) {
3450  // cout << name << " perf " << p << " " << j << endl;
3451  //}
3452 
3453  // dQ / dS
3454  for (USI k = 0; k < np; k++) {
3455  tmp = CONV1 * perf[p].WI * perf[p].multiplier * dP / mu * xi * xij *
3456  myBulk.dKr_dS[pId * np * np + j * np + k];
3457  // capillary pressure
3458  tmp += transIJ * myBulk.dPcj_dS[pId * np * np + j * np + k];
3459  dQdXsB[(i + 1) * ncol2 + k] += tmp;
3460  }
3461  // dQ / dCij
3462  for (USI k = 0; k < nc; k++) {
3463  tmp = dP * perf[p].transj[j] * xij *
3464  (myBulk.xix[n_np_j * nc + k] -
3465  xi / mu * myBulk.mux[n_np_j * nc + k]);
3466  if (k == i) {
3467  tmp += perf[p].transj[j] * xi * dP;
3468  }
3469  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] += tmp;
3470  }
3471  }
3472  }
3473  // Bulk to Well
3474  bmat = dQdXpB;
3475  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &myBulk.dSec_dPri[pId * bsize2], 1,
3476  bmat.data());
3477 
3478  Dscalar(bsize, dt, bmat.data());
3479  // Add
3480  USI ptr = myLS.diagPtr[n];
3481  for (USI i = 0; i < bsize; i++) {
3482  myLS.val[n][ptr * bsize + i] += bmat[i];
3483  }
3484  // Insert
3485  bmat = dQdXpW;
3486  Dscalar(bsize, dt, bmat.data());
3487  myLS.val[n].insert(myLS.val[n].end(), bmat.begin(), bmat.end());
3488  myLS.colId[n].push_back(wId);
3489 
3490  // Well
3491  switch (opt.optMode) {
3492  case RATE_MODE:
3493  case ORATE_MODE:
3494  case GRATE_MODE:
3495  case WRATE_MODE:
3496  case LRATE_MODE:
3497  // Diag
3498  fill(bmat.begin(), bmat.end(), 0.0);
3499  for (USI i = 0; i < nc; i++) {
3500  bmat[0] += dQdXpW[(i + 1) * ncol] * prodWeight[i];
3501  bmat[(i + 1) * ncol + i + 1] = 1;
3502  }
3503  for (USI i = 0; i < bsize; i++) {
3504  myLS.diagVal[wId * bsize + i] += bmat[i];
3505  }
3506 
3507  // OffDiag
3508  bmat = dQdXpB;
3509  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(),
3510  &myBulk.dSec_dPri[pId * bsize2], 1, bmat.data());
3511  fill(bmat2.begin(), bmat2.end(), 0.0);
3512  for (USI i = 0; i < nc; i++) {
3513  Daxpy(ncol, prodWeight[i], bmat.data() + (i + 1) * ncol,
3514  bmat2.data());
3515  }
3516  myLS.val[wId].insert(myLS.val[wId].end(), bmat2.begin(), bmat2.end());
3517  myLS.colId[wId].push_back(n);
3518  break;
3519 
3520  case BHP_MODE:
3521  // Diag
3522  fill(bmat.begin(), bmat.end(), 0.0);
3523  for (USI i = 0; i < nc + 1; i++) {
3524  bmat[i * ncol + i] = 1;
3525  }
3526  // Add
3527  for (USI i = 0; i < bsize; i++) {
3528  myLS.diagVal[wId * bsize + i] += bmat[i];
3529  }
3530  // OffDiag
3531  fill(bmat.begin(), bmat.end(), 0.0);
3532  // Insert
3533  myLS.val[wId].insert(myLS.val[wId].end(), bmat.begin(), bmat.end());
3534  myLS.colId[wId].push_back(n);
3535  // Solution
3536  // myLS.u[wId * ncol] = opt.minBHP - BHP;
3537  break;
3538 
3539  default:
3540  OCP_ABORT("Wrong Well Opt mode!");
3541  break;
3542  }
3543  }
3544  assert(myLS.val[wId].size() == numPerf * bsize);
3545  // Well self
3546  myLS.colId[wId].push_back(wId);
3547  myLS.diagPtr[wId] = numPerf;
3548  myLS.val[wId].insert(myLS.val[wId].end(), myLS.diagVal.data() + wId * bsize,
3549  myLS.diagVal.data() + wId * bsize + bsize);
3550 
3551  // test
3552  // for (USI i = 0; i < bsize; i++) {
3553  // cout << myLS.diagVal[wId * bsize + i] << endl;
3554  //}
3555 }
3556 
3557 /*----------------------------------------------------------------------------*/
3558 /* Brief Change History of This File */
3559 /*----------------------------------------------------------------------------*/
3560 /* Author Date Actions */
3561 /*----------------------------------------------------------------------------*/
3562 /* Shizhe Li Oct/01/2021 Create file */
3563 /* Chensong Zhang Oct/15/2021 Format file */
3564 /*----------------------------------------------------------------------------*/
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:76
Dscalar
void Dscalar(const int &n, const double &alpha, double *x)
Scales a vector by a constant.
Definition: DenseMat.cpp:62
DaAxpby
void DaAxpby(const int &m, const int &n, const double &a, const double *A, const double *x, const double &b, double *y)
Computes y = a A x + b y.
Definition: DenseMat.cpp:173
Daxpy
void Daxpy(const int &n, const double &alpha, const double *x, double *y)
Constant times a vector plus a vector.
Definition: DenseMat.cpp:69
Dnorm1
double Dnorm1(const int &N, double *x)
Computes the L1-norm of a vector.
Definition: DenseMat.cpp:50
X_DIRECTION
const USI X_DIRECTION
x-direction
Definition: OCPConst.hpp:123
BLKOIL
const USI BLKOIL
Mixture model = black-oil.
Definition: OCPConst.hpp:88
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:36
EOS_PVTW
const USI EOS_PVTW
Mixture model = equation-of-state.
Definition: OCPConst.hpp:89
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
GAS
const USI GAS
Fluid type = gas.
Definition: OCPConst.hpp:83
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
CONV1
const OCP_DBL CONV1
1 bbl = CONV1 ft3
Definition: OCPConst.hpp:59
GRAVITY_FACTOR
const OCP_DBL GRAVITY_FACTOR
0.00694444 ft2 psi / lb
Definition: OCPConst.hpp:52
LRATE_MODE
const USI LRATE_MODE
Well option = fixed fluid rate???
Definition: OCPConst.hpp:119
WRATE_MODE
const USI WRATE_MODE
Well option = fixed water rate.
Definition: OCPConst.hpp:118
WATER
const USI WATER
Fluid type = water.
Definition: OCPConst.hpp:84
PHASE_ODGW
const USI PHASE_ODGW
Phase type = oil-dry gas-water.
Definition: OCPConst.hpp:100
Y_DIRECTION
const USI Y_DIRECTION
y-direction
Definition: OCPConst.hpp:124
ORATE_MODE
const USI ORATE_MODE
Well option = fixed oil rate.
Definition: OCPConst.hpp:116
OIL
const USI OIL
Fluid type = oil.
Definition: OCPConst.hpp:82
CLOSE
const bool CLOSE
Well type = closed.
Definition: OCPConst.hpp:112
PHASE_W
const USI PHASE_W
Phase type = water only.
Definition: OCPConst.hpp:96
RATE_MODE
const USI RATE_MODE
Well option = fixed total rate???
Definition: OCPConst.hpp:115
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:24
OPEN
const bool OPEN
Well type = open.
Definition: OCPConst.hpp:111
Z_DIRECTION
const USI Z_DIRECTION
z-direction
Definition: OCPConst.hpp:125
PHASE_OW
const USI PHASE_OW
Phase type = oil-water.
Definition: OCPConst.hpp:98
PROD
const USI PROD
Well type = producer.
Definition: OCPConst.hpp:108
TEMPERATURE_STD
const OCP_DBL TEMPERATURE_STD
Standard temperature.
Definition: OCPConst.hpp:56
PHASE_DOGW
const USI PHASE_DOGW
Phase type = dead oil-gas-water.
Definition: OCPConst.hpp:101
PI
const OCP_DBL PI
Pi.
Definition: OCPConst.hpp:37
OCP_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:25
INJ
const USI INJ
Well type = injector.
Definition: OCPConst.hpp:107
BHP_MODE
const USI BHP_MODE
Well option = fixed bottom-hole-pressure.
Definition: OCPConst.hpp:120
GRATE_MODE
const USI GRATE_MODE
Well option = fixed gas rate.
Definition: OCPConst.hpp:117
PRESSURE_STD
const OCP_DBL PRESSURE_STD
14.6959 psia = 1 atm
Definition: OCPConst.hpp:55
CONV2
const OCP_DBL CONV2
Darcy constant in field unit.
Definition: OCPConst.hpp:60
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_ASSERT
#define OCP_ASSERT(cond, msg)
Assert condition and log user messages in DEBUG mode.
Definition: UtilError.hpp:58
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:58
Grid
Basic information of computational grid, including the rock properties.
Definition: Grid.hpp:76
LinearSystem
Linear solvers for discrete systems.
Definition: LinearSystem.hpp:31
WellOptParam::maxBHP
OCP_DBL maxBHP
Maximum allowable pressure in the injection well.
Definition: ParamWell.hpp:39
WellOptParam::type
string type
Type of well, injection or production?
Definition: ParamWell.hpp:33
WellOptParam::fluidType
string fluidType
Type of fluid into the injection well. (injection well only)
Definition: ParamWell.hpp:34
WellOptParam::maxRate
OCP_DBL maxRate
Maximum allowable flow rate into/out the well.
Definition: ParamWell.hpp:38
WellOptParam::optMode
string optMode
Mode of well, Rate or BHP?
Definition: ParamWell.hpp:36
WellOptParam::minBHP
OCP_DBL minBHP
Minimum allowable pressure in the production well.
Definition: ParamWell.hpp:40
WellOptParam::state
string state
State of well, open or close?
Definition: ParamWell.hpp:35
WellOpt
Definition: Well.hpp:49
WellOpt::WellOpt
WellOpt()=default
Default constructor.
WellOpt::operator!=
bool operator!=(const WellOpt &Opt) const
overload inequality
Definition: Well.cpp:63
WellParam
TODO: Add Doxygen.
Definition: ParamWell.hpp:61
WellParam::direction
vector< string > direction
Direction of perforations.
Definition: ParamWell.hpp:80
WellParam::diameter
vector< OCP_DBL > diameter
Diameter of perforations.
Definition: ParamWell.hpp:77
WellParam::skinFactor
vector< OCP_DBL > skinFactor
Skin factor.
Definition: ParamWell.hpp:79
WellParam::I_perf
vector< USI > I_perf
I-index of perforation in grid.
Definition: ParamWell.hpp:73
WellParam::WI
vector< OCP_DBL > WI
Transmissiblity connection factor.
Definition: ParamWell.hpp:76
WellParam::J_perf
vector< USI > J_perf
J-index of perforation in grid.
Definition: ParamWell.hpp:74
WellParam::K_perf
vector< USI > K_perf
K-index of perforation in grid.
Definition: ParamWell.hpp:75
Well::AssembleMatPROD_AIMt
void AssembleMatPROD_AIMt(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for AIMt, parts related to Production well are included.
Definition: Well.cpp:2990
Well::AssembleMatPROD_IMPEC
void AssembleMatPROD_IMPEC(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for IMPEC, parts related to production well are included.
Definition: Well.cpp:1338
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:1445
Well::SetBHP
void SetBHP()
Set BHP if opt mode is BHPMode.
Definition: Well.cpp:949
Well::Setup
void Setup(const Grid &myGrid, const Bulk &myBulk, const vector< SolventINJ > &sols)
Setup the well after Grid and Bulk finish setupping.
Definition: Well.cpp:105
Well::CheckOptMode
void CheckOptMode(const Bulk &myBulk)
Check if well operation mode would be changed.
Definition: Well.cpp:971
Well::CalProdQjCOMP
void CalProdQjCOMP(const Bulk &myBulk)
Calculate flow rate of moles of phase for production well in Compositional Model.
Definition: Well.cpp:548
Well::CalWI_Peaceman_Vertical
void CalWI_Peaceman_Vertical(const Bulk &myBulk)
Calculate Well Index with Peaceman model for vertical well.
Definition: Well.cpp:240
Well::CheckP
OCP_INT CheckP(const Bulk &myBulk)
Check if abnormal Pressure occurs.
Definition: Well.cpp:1045
Well::AssembleMatINJ_AIMt
void AssembleMatINJ_AIMt(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for AIMt, parts related to Injection well are included.
Definition: Well.cpp:2846
Well::CalInjdG
void CalInjdG(const Bulk &myBulk)
Calculate pressure difference between well and perforations for Injection.
Definition: Well.cpp:608
Well::AssembleMatINJ_FIM
void AssembleMatINJ_FIM(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for FIM, parts related to Injection well are included.
Definition: Well.cpp:1528
Well::AssembleMatPROD_FIM_new
void AssembleMatPROD_FIM_new(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for FIM, parts related to Production well are included.
Definition: Well.cpp:2291
Well::CalProdQj
void CalProdQj(const Bulk &myBulk, const OCP_DBL &dt)
Calculate flow rate of moles of phase for production well.
Definition: Well.cpp:535
Well::AssembleMatINJ_FIM_new
void AssembleMatINJ_FIM_new(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for FIM, parts related to Injection well are included.
Definition: Well.cpp:2123
Well::SetupWellBulk
void SetupWellBulk(Bulk &myBulk) const
Setup bulks which are penetrated by wells.
Definition: Well.cpp:1176
Well::MassConserveIMPEC
void MassConserveIMPEC(Bulk &myBulk, const OCP_DBL &dt) const
Update moles of components in those bulks who connects to the well.
Definition: Well.cpp:1235
Well::CalReInjFluid
void CalReInjFluid(const Bulk &myBulk, vector< OCP_DBL > &myZi)
Calculate the contribution of production well to reinjection defaulted.
Definition: Well.cpp:930
Well::InputPerfo
void InputPerfo(const WellParam &well)
Input the param of perforations.
Definition: Well.cpp:79
Well::AssembleMatPROD_FIM
void AssembleMatPROD_FIM(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for FIM, parts related to Production well are included.
Definition: Well.cpp:1673
Well::AssembleMatINJ_AIMs
void AssembleMatINJ_AIMs(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for AIMs, parts related to Injection well are included.
Definition: Well.cpp:3247
Well::CalProdQiBO
void CalProdQiBO(const Bulk &myBulk)
Definition: Well.cpp:557
Well::WellState
bool WellState() const
Return the state of the well, Open or Close.
Definition: Well.hpp:172
Well::AssembleMatPROD_AIMs
void AssembleMatPROD_AIMs(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for AIMs, parts related to Production well are included.
Definition: Well.cpp:3390
Well::CalTrans
void CalTrans(const Bulk &myBulk)
Calculate transmissibility for each phase in perforations.
Definition: Well.cpp:311
Well::CalFlux
void CalFlux(const Bulk &myBulk, const bool flag=false)
Calculate the flux for each perforations.
Definition: Well.cpp:373
Well::CalProddG
void CalProddG(const Bulk &myBulk)
Calculate pressure difference between well and perforations for Prodcution.
Definition: Well.cpp:678
Well::AllocateMat
void AllocateMat(LinearSystem &myLS) const
Allocate memory for matrix.
Definition: Well.cpp:1167
Well::CaldG
void CaldG(const Bulk &myBulk)
Calculate pressure difference between well and perforations.
Definition: Well.cpp:577
Well::AssembleMatINJ_IMPEC
void AssembleMatINJ_IMPEC(const Bulk &myBulk, LinearSystem &myLS, const OCP_DBL &dt) const
Assemble matrix for IMPEC, parts related to injection well are included.
Definition: Well.cpp:1249
Well::CheckCrossFlow
OCP_INT CheckCrossFlow(const Bulk &myBulk)
Check if crossflow happens.
Definition: Well.cpp:1091
Well::SmoothdG
void SmoothdG()
Try to smooth the dG by average it with dG at last time step.
Definition: Well.cpp:959
Well::CalResAIMt
void CalResAIMt(ResFIM &resFIM, const Bulk &myBulk, const OCP_DBL &dt, const OCP_USI &wId, const vector< Well > &allWell) const
Calculate Resiual and relative Resiual for local FIM.
Definition: Well.cpp:3157
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:1842
Well::ShowPerfStatus
void ShowPerfStatus(const Bulk &myBulk) const
Display operation mode of well and state of perforations.
Definition: Well.cpp:1184
Well::InitBHP
void InitBHP(const Bulk &myBulk)
Initialize the Well BHP.
Definition: Well.cpp:238
Well::CalInjQi
void CalInjQi(const Bulk &myBulk, const OCP_DBL &dt)
Calculate flow rate of moles of components for injection well.
Definition: Well.cpp:511
Well::CalProdRate
OCP_DBL CalProdRate(const Bulk &myBulk, const bool &minBHP)
calculate flow rate of moles of components for production well with minBHP
Definition: Well.cpp:481
Well::CalProdWeight
void CalProdWeight(const Bulk &myBulk) const
Calculate the Prodweight.
Definition: Well.cpp:833
Well::CalCFL
void CalCFL(const Bulk &myBulk, const OCP_DBL &dt) const
Calculate the CFL number, only parts related to wells are considered.
Definition: Well.cpp:1217
Well::CalInjRate
OCP_DBL CalInjRate(const Bulk &myBulk, const bool &maxBHP)
calculate flow rate of moles of components for injection well with maxBHP
Definition: Well.cpp:460
Well::CalResFIM
void CalResFIM(ResFIM &resFIM, const Bulk &myBulk, const OCP_DBL &dt, const OCP_USI &wId, const vector< Well > &allWell) const
Calculate Resiual and relative Resiual for FIM.
Definition: Well.cpp:1983
Generated on Sat Oct 15 2022 00:08:52 for OpenCAEPoro by doxygen 1.9.1