OpenCAEPoro  v0.5.0
A simulator for multicomponent porous media flow
OCPFluidMethod.cpp
Go to the documentation of this file.
1 
12 #include "OCPFluidMethod.hpp"
13 
15 // IsothermalMethod
17 
18 void IsothermalMethod::InitRock(Bulk& bk) const
19 {
20  for (OCP_USI n = 0; n < bk.numBulk; n++) {
21  bk.poroInit[n] *= bk.ntg[n];
22  }
23 }
24 
25 void IsothermalMethod::CalRock(Bulk& bk) const
26 {
27  for (OCP_USI n = 0; n < bk.numBulk; n++) {
28  bk.rock[bk.ROCKNUM[n]]->CalPoro(bk.P[n], bk.poroInit[n], bk.poro[n],
29  bk.poroP[n]);
30  bk.rockVp[n] = bk.v[n] * bk.poro[n];
31  }
32 }
33 
35 // IsoT_IMPEC
37 
38 void IsoT_IMPEC::Setup(Reservoir& rs, LinearSystem& ls, const OCPControl& ctrl)
39 {
40  // Allocate Memory of auxiliary variables for IMPEC
41  AllocateReservoir(rs);
42  // Allocate Memory of Matrix for IMPEC
43  AllocateLinearSystem(ls, rs, ctrl);
44 }
45 
47 void IsoT_IMPEC::InitReservoir(Reservoir& rs) const
48 {
49  rs.bulk.InitPTSw(50);
50 
51  InitRock(rs.bulk);
52  CalRock(rs.bulk);
53 
54  InitFlash(rs.bulk);
55  CalKrPc(rs.bulk);
56 
57  CalBulkFlux(rs);
58 
59  rs.allWells.InitBHP(rs.bulk);
60 
61  UpdateLastTimeStep(rs);
62 }
63 
64 void IsoT_IMPEC::Prepare(Reservoir& rs, OCPControl& ctrl)
65 {
66  rs.allWells.PrepareWell(rs.bulk);
67  rs.CalCFL(ctrl.GetCurDt());
68  ctrl.Check(rs, {"CFL"});
69 }
70 
71 void IsoT_IMPEC::AssembleMat(LinearSystem& ls,
72  const Reservoir& rs,
73  const OCP_DBL& dt) const
74 {
75  AssembleMatBulks(ls, rs, dt);
76  AssembleMatWells(ls, rs, dt);
77 }
78 
79 void IsoT_IMPEC::SolveLinearSystem(LinearSystem& ls, Reservoir& rs, OCPControl& ctrl)
80 {
81 #ifdef DEBUG
82  ls.CheckEquation();
83 #endif // DEBUG
84 
85  ls.AssembleMatLinearSolver();
86 
87 #ifdef DEBUG
88  // ls.OutputLinearSystem("testA_IMPEC.out", "testb_IMPEC.out");
89 #endif // DEBUG
90 
91  GetWallTime Timer;
92  Timer.Start();
93  int status = ls.Solve();
94  if (status < 0) {
95  status = ls.GetNumIters();
96  }
97 
98  ctrl.RecordTimeLS(Timer.Stop() / 1000);
99  ctrl.UpdateIterLS(status);
100  ctrl.UpdateIterNR();
101 
102 #ifdef DEBUG
103  // ls.OutputSolution("testx_IMPEC.out");
104 #endif // DEBUG
105 
106  // rs.GetSolutionIMPEC(ls.GetSolution());
107  GetSolution(rs, ls.GetSolution());
108  ls.ClearData();
109 }
110 
111 OCP_BOOL IsoT_IMPEC::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
112 {
113  OCP_DBL& dt = ctrl.current_dt;
114 
115  // First check : Pressure check
116  if (!ctrl.Check(rs, {"BulkP", "WellP"})) {
117  return OCP_FALSE;
118  }
119 
120  // Calculate Flux between bulks and between bulks and wells
121  CalFlux(rs);
122 
123  // Second check : CFL check
124  rs.CalCFL(dt);
125  if (!ctrl.Check(rs, {"CFL"})) {
126  ResetToLastTimeStep01(rs, ctrl);
127  cout << "CFL is too big" << endl;
128  return OCP_FALSE;
129  }
130 
131  MassConserve(rs, dt);
132 
133  // Third check: Ni check
134  if (!ctrl.Check(rs, {"BulkNi"})) {
135  ResetToLastTimeStep02(rs, ctrl);
136  return OCP_FALSE;
137  }
138 
139  CalRock(rs.bulk);
140  CalFlash(rs.bulk);
141 
142  // Fouth check: Volume error check
143  if (!ctrl.Check(rs, {"BulkVe"})) {
144  ResetToLastTimeStep03(rs, ctrl);
145  return OCP_FALSE;
146  }
147 
148  CalKrPc(rs.bulk);
149  CalBulkFlux(rs);
150 
151  return OCP_TRUE;
152 }
153 
154 OCP_BOOL IsoT_IMPEC::FinishNR(const Reservoir& rs) { return OCP_TRUE; }
155 
156 void IsoT_IMPEC::FinishStep(Reservoir& rs, OCPControl& ctrl)
157 {
158  rs.CalIPRT(ctrl.GetCurDt());
159  rs.CalMaxChange();
160  UpdateLastTimeStep(rs);
161  // ctrl.CalNextTstepIMPEC(rs);
162  ctrl.CalNextTimeStep(rs, {"dP", "dN", "dS", "eV"});
163 }
164 
165 void IsoT_IMPEC::AllocateReservoir(Reservoir& rs)
166 {
167  Bulk& bk = rs.bulk;
168  const OCP_USI nb = bk.numBulk;
169  const USI np = bk.numPhase;
170  const USI nc = bk.numCom;
171 
172  // Rock
173  bk.poro.resize(nb);
174  bk.rockVp.resize(nb);
175 
176  bk.lporo.resize(nb);
177  bk.lrockVp.resize(nb);
178 
179  // derivatives
180  bk.poroP.resize(nb);
181  bk.lporoP.resize(nb);
182 
183  // Fluid
184  bk.phaseNum.resize(nb);
185  bk.Nt.resize(nb);
186  bk.Ni.resize(nb * nc);
187  bk.vf.resize(nb);
188  bk.T.resize(nb);
189  bk.P.resize(nb);
190  bk.Pb.resize(nb);
191  bk.Pj.resize(nb * np);
192  bk.Pc.resize(nb * np);
193  bk.phaseExist.resize(nb * np);
194  bk.S.resize(nb * np);
195  bk.vj.resize(nb * np);
196  bk.xij.resize(nb * np * nc);
197  bk.rho.resize(nb * np);
198  bk.xi.resize(nb * np);
199  bk.mu.resize(nb * np);
200  bk.kr.resize(nb * np);
201 
202  bk.lphaseNum.resize(nb);
203  bk.lNt.resize(nb);
204  bk.lNi.resize(nb * nc);
205  bk.lvf.resize(nb);
206  bk.lT.resize(nb);
207  bk.lP.resize(nb);
208  bk.lPj.resize(nb * np);
209  bk.lPc.resize(nb * np);
210  bk.lphaseExist.resize(nb * np);
211  bk.lS.resize(nb * np);
212  bk.vj.resize(nb * np);
213  bk.lxij.resize(nb * np * nc);
214  bk.lrho.resize(nb * np);
215  bk.lxi.resize(nb * np);
216  bk.lmu.resize(nb * np);
217  bk.lkr.resize(nb * np);
218 
219  // derivatives
220  bk.vfP.resize(nb);
221  bk.vfi.resize(nb * nc);
222 
223  bk.lvfP.resize(nb);
224  bk.lvfi.resize(nb * nc);
225 
226  // others
227  bk.cfl.resize(nb * np);
228 
229  BulkConn& conn = rs.conn;
230 
231  conn.upblock.resize(conn.numConn * np);
232  conn.upblock_Rho.resize(conn.numConn * np);
233  conn.upblock_Trans.resize(conn.numConn * np);
234  conn.upblock_Velocity.resize(conn.numConn * np);
235 
236  conn.lupblock.resize(conn.numConn * np);
237  conn.lupblock_Rho.resize(conn.numConn * np);
238  conn.lupblock_Trans.resize(conn.numConn * np);
239  conn.lupblock_Velocity.resize(conn.numConn * np);
240 }
241 
242 void IsoT_IMPEC::AllocateLinearSystem(LinearSystem& ls,
243  const Reservoir& rs,
244  const OCPControl& ctrl)
245 {
246  ls.AllocateRowMem(rs.GetBulkNum() + rs.GetWellNum(), 1);
247  ls.AllocateColMem(rs.conn.GetNeighborNum(), rs.allWells.GetWell2Bulk());
248  ls.SetupLinearSolver(SCALARFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
249 }
250 
251 void IsoT_IMPEC::InitFlash(Bulk& bk) const
252 {
253  for (OCP_USI n = 0; n < bk.numBulk; n++) {
254  bk.flashCal[bk.PVTNUM[n]]->InitFlashIMPEC(bk.P[n], bk.Pb[n], bk.T[n],
255  &bk.S[n * bk.numPhase], bk.rockVp[n],
256  bk.Ni.data() + n * bk.numCom, n);
257  for (USI i = 0; i < bk.numCom; i++) {
258  bk.Ni[n * bk.numCom + i] = bk.flashCal[bk.PVTNUM[n]]->GetNi(i);
259  }
260  PassFlashValue(bk, n);
261  }
262 }
263 
264 void IsoT_IMPEC::CalFlash(Bulk& bk)
265 {
266  for (OCP_USI n = 0; n < bk.numBulk; n++) {
267 
268  bk.flashCal[bk.PVTNUM[n]]->FlashIMPEC(bk.P[n], bk.T[n], &bk.Ni[n * bk.numCom],
269  bk.phaseNum[n],
270  &bk.xij[n * bk.numPhase * bk.numCom], n);
271  PassFlashValue(bk, n);
272  }
273 }
274 
275 void IsoT_IMPEC::PassFlashValue(Bulk& bk, const OCP_USI& n) const
276 {
277  const USI np = bk.numPhase;
278  const USI nc = bk.numCom;
279  const OCP_USI bIdp = n * np;
280  const USI pvtnum = bk.PVTNUM[n];
281 
282  bk.phaseNum[n] = 0;
283  bk.Nt[n] = bk.flashCal[pvtnum]->GetNt();
284  bk.vf[n] = bk.flashCal[pvtnum]->GetVf();
285 
286  for (USI j = 0; j < np; j++) {
287  // Important! Saturation must be passed no matter if the phase exists. This is
288  // because it will be used to calculate relative permeability and capillary
289  // pressure at each time step. Make sure that all saturations are updated at
290  // each step!
291  bk.phaseExist[bIdp + j] = bk.flashCal[pvtnum]->GetPhaseExist(j);
292  bk.S[bIdp + j] = bk.flashCal[pvtnum]->GetS(j);
293  if (bk.phaseExist[bIdp + j]) {
294  bk.phaseNum[n]++;
295  for (USI i = 0; i < nc; i++) {
296  bk.xij[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetXij(j, i);
297  }
298  bk.vj[bIdp + j] = bk.flashCal[pvtnum]->GetVj(j);
299  bk.rho[bIdp + j] = bk.flashCal[pvtnum]->GetRho(j);
300  bk.xi[bIdp + j] = bk.flashCal[pvtnum]->GetXi(j);
301  bk.mu[bIdp + j] = bk.flashCal[pvtnum]->GetMu(j);
302  }
303  }
304 
305  bk.vfP[n] = bk.flashCal[pvtnum]->GetVfP();
306  for (USI i = 0; i < nc; i++) {
307  bk.vfi[n * nc + i] = bk.flashCal[pvtnum]->GetVfi(i);
308  }
309 }
310 
311 void IsoT_IMPEC::CalKrPc(Bulk& bk) const
312 {
313  for (OCP_USI n = 0; n < bk.numBulk; n++) {
314  OCP_USI bId = n * bk.numPhase;
315  bk.flow[bk.SATNUM[n]]->CalKrPc(&bk.S[bId], &bk.kr[bId], &bk.Pc[bId], n);
316  for (USI j = 0; j < bk.numPhase; j++)
317  bk.Pj[n * bk.numPhase + j] = bk.P[n] + bk.Pc[n * bk.numPhase + j];
318  }
319 }
320 
321 void IsoT_IMPEC::CalFlux(Reservoir& rs) const
322 {
323  CalBulkFlux(rs);
324  rs.allWells.CalFlux(rs.bulk);
325 }
326 
327 void IsoT_IMPEC::CalBulkFlux(Reservoir& rs) const
328 {
329  const Bulk& bk = rs.bulk;
330  BulkConn& conn = rs.conn;
331  const USI np = bk.numPhase;
332 
333  // calculate a step flux using iteratorConn
334  OCP_USI bId, eId, uId;
335  OCP_USI bId_np_j, eId_np_j;
336  OCP_BOOL exbegin, exend, exup;
337  OCP_DBL rho, dP, Akd;
338 
339  for (OCP_USI c = 0; c < conn.numConn; c++) {
340  bId = conn.iteratorConn[c].BId();
341  eId = conn.iteratorConn[c].EId();
342  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
343 
344  for (USI j = 0; j < np; j++) {
345  bId_np_j = bId * np + j;
346  eId_np_j = eId * np + j;
347 
348  exbegin = bk.phaseExist[bId_np_j];
349  exend = bk.phaseExist[eId_np_j];
350 
351  if ((exbegin) && (exend)) {
352  rho = (bk.rho[bId_np_j] + bk.rho[eId_np_j]) / 2;
353  } else if (exbegin && (!exend)) {
354  rho = bk.rho[bId_np_j];
355  } else if ((!exbegin) && (exend)) {
356  rho = bk.rho[eId_np_j];
357  } else {
358  conn.upblock[c * np + j] = bId;
359  continue;
360  }
361 
362  dP = (bk.Pj[bId_np_j] - GRAVITY_FACTOR * rho * bk.depth[bId]) -
363  (bk.Pj[eId_np_j] - GRAVITY_FACTOR * rho * bk.depth[eId]);
364  if (dP < 0) {
365  uId = eId;
366  exup = exend;
367  } else {
368  uId = bId;
369  exup = exbegin;
370  }
371 
372  conn.upblock_Rho[c * np + j] = rho;
373  conn.upblock[c * np + j] = uId;
374 
375  if (exup) {
376  conn.upblock_Trans[c * np + j] =
377  Akd * bk.kr[uId * np + j] / bk.mu[uId * np + j];
378  conn.upblock_Velocity[c * np + j] = conn.upblock_Trans[c * np + j] * dP;
379  } else {
380  conn.upblock_Trans[c * np + j] = 0;
381  conn.upblock_Velocity[c * np + j] = 0;
382  }
383  }
384  }
385 }
386 
387 void IsoT_IMPEC::MassConserve(Reservoir& rs, const OCP_DBL& dt) const
388 {
389 
390  // Bulk to Bulk
391  Bulk& bk = rs.bulk;
392  const BulkConn& conn = rs.conn;
393 
394  const USI np = bk.numPhase;
395  const USI nc = bk.numCom;
396 
397  OCP_USI bId, eId, uId;
398  OCP_USI uId_np_j;
399  OCP_DBL phaseVelocity, dNi;
400 
401  for (OCP_USI c = 0; c < conn.numConn; c++) {
402  bId = conn.iteratorConn[c].BId();
403  eId = conn.iteratorConn[c].EId();
404 
405  for (USI j = 0; j < np; j++) {
406  uId = conn.upblock[c * np + j];
407  uId_np_j = uId * np + j;
408 
409  if (!bk.phaseExist[uId_np_j]) continue;
410 
411  phaseVelocity = conn.upblock_Velocity[c * np + j];
412  for (USI i = 0; i < nc; i++) {
413  dNi = dt * phaseVelocity * bk.xi[uId_np_j] * bk.xij[uId_np_j * nc + i];
414  bk.Ni[eId * nc + i] += dNi;
415  bk.Ni[bId * nc + i] -= dNi;
416  }
417  }
418  }
419 
420  // Well to Bulk
421  for (auto& wl : rs.allWells.wells) {
422  if (wl.IsOpen()) {
423  for (USI p = 0; p < wl.PerfNum(); p++) {
424  OCP_USI k = wl.PerfLocation(p);
425  for (USI i = 0; i < nc; i++) {
426  bk.Ni[k * nc + i] -= wl.PerfQi_lbmol(p, i) * dt;
427  }
428  }
429  }
430  }
431 }
432 
433 void IsoT_IMPEC::AssembleMatBulks(LinearSystem& ls,
434  const Reservoir& rs,
435  const OCP_DBL& dt) const
436 {
437 
438  const Bulk& bk = rs.bulk;
439  const BulkConn& conn = rs.conn;
440  const OCP_USI nb = bk.numBulk;
441  const USI np = bk.numPhase;
442  const USI nc = bk.numCom;
443 
444  ls.AddDim(nb);
445 
446  // accumulate term
447  OCP_DBL Vpp, Vp, vf, vfP, P;
448  for (OCP_USI n = 0; n < nb; n++) {
449  vf = bk.vf[n];
450  vfP = bk.vfP[n];
451  P = bk.lP[n];
452  Vpp = bk.v[n] * bk.poroP[n];
453  Vp = bk.rockVp[n];
454 
455  ls.NewDiag(n, Vpp - vfP);
456  ls.AddRhs(n, (Vpp - vfP) * P + dt * (vf - Vp));
457  }
458 
459  // flux term
460  OCP_USI bId, eId, uId_np_j;
461  OCP_DBL valupi, valdowni, valup, rhsup, valdown, rhsdown;
462  OCP_DBL dD, tmp;
463 
464  // Be careful when first bulk has no neighbors!
465  for (OCP_USI c = 0; c < conn.numConn; c++) {
466  bId = conn.iteratorConn[c].BId();
467  eId = conn.iteratorConn[c].EId();
468  dD = GRAVITY_FACTOR * (bk.depth[bId] - bk.depth[eId]);
469 
470  valup = 0;
471  rhsup = 0;
472  valdown = 0;
473  rhsdown = 0;
474 
475  for (USI j = 0; j < np; j++) {
476  uId_np_j = conn.upblock[c * np + j] * np + j;
477  if (!bk.phaseExist[uId_np_j]) continue;
478 
479  valupi = 0;
480  valdowni = 0;
481 
482  for (USI i = 0; i < nc; i++) {
483  valupi += bk.vfi[bId * nc + i] * bk.xij[uId_np_j * nc + i];
484  valdowni += bk.vfi[eId * nc + i] * bk.xij[uId_np_j * nc + i];
485  }
486 
487  tmp = bk.xi[uId_np_j] * conn.upblock_Trans[c * np + j] * dt;
488  valup += tmp * valupi;
489  valdown += tmp * valdowni;
490  tmp *= conn.upblock_Rho[c * np + j] * dD -
491  (bk.Pc[bId * np + j] - bk.Pc[eId * np + j]);
492  rhsup += tmp * valupi;
493  rhsdown -= tmp * valdowni;
494  }
495  ls.AddDiag(bId, valup);
496  ls.AddDiag(eId, valdown);
497  ls.NewOffDiag(bId, eId, -valup);
498  ls.NewOffDiag(eId, bId, -valdown);
499  ls.AddRhs(bId, rhsup);
500  ls.AddRhs(eId, rhsdown);
501  }
502 }
503 
504 void IsoT_IMPEC::AssembleMatWells(LinearSystem& ls,
505  const Reservoir& rs,
506  const OCP_DBL& dt) const
507 {
508  for (auto& wl : rs.allWells.wells) {
509  if (wl.IsOpen()) {
510  switch (wl.WellType()) {
511  case INJ:
512  AssembleMatInjWells(ls, rs.bulk, wl, dt);
513  break;
514  case PROD:
515  AssembleMatProdWells(ls, rs.bulk, wl, dt);
516  break;
517  default:
518  OCP_ABORT("Wrong well type");
519  }
520  }
521  }
522 
523  // for Reinjection
524  // for (auto& wG : wellGroup) {
525  // if (wG.reInj) {
526  // for (auto& prod : wellGroup[wG.prodGroup].wIdPROD) {
527  // if (wells[prod].IsOpen()) {
528  // wells[prod].AssembleMatReinjection_IMPEC(myBulk, myLS, dt, wells,
529  // wG.wIdINJ);
530  // }
531  // }
532  // }
533  //}
534 }
535 
536 void IsoT_IMPEC::AssembleMatInjWells(LinearSystem& ls,
537  const Bulk& bk,
538  const Well& wl,
539  const OCP_DBL& dt) const
540 {
541  const OCP_USI wId = ls.AddDim(1) - 1;
542  ls.NewDiag(wId, 0.0);
543 
544  OCP_DBL Vfi_zi, valb, valw, bb, bw;
545 
546  const USI nc = bk.numCom;
547 
548  for (USI p = 0; p < wl.PerfNum(); p++) {
549  const OCP_USI n = wl.PerfLocation(p);
550 
551  Vfi_zi = 0;
552  for (USI i = 0; i < nc; i++) {
553  Vfi_zi += bk.vfi[n * nc + i] * wl.InjZi(i);
554  }
555 
556  valw = dt * wl.PerfXi(p) * wl.PerfTransInj(p);
557  bw = valw * wl.DG(p);
558  valb = valw * Vfi_zi;
559  bb = valb * wl.DG(p);
560 
561  // Bulk to Well
562  ls.AddDiag(n, valb);
563  ls.NewOffDiag(n, wId, -valb);
564  ls.AddRhs(n, bb);
565 
566  // Well to Bulk
567  switch (wl.OptMode()) {
568  case RATE_MODE:
569  case ORATE_MODE:
570  case GRATE_MODE:
571  case WRATE_MODE:
572  case LRATE_MODE:
573  ls.AddDiag(wId, valw);
574  ls.NewOffDiag(wId, n, -valw);
575  ls.AddRhs(wId, -bw);
576  break;
577  case BHP_MODE:
578  ls.NewOffDiag(wId, n, 0);
579  break;
580  default:
581  OCP_ABORT("Wrong well option mode!");
582  }
583  }
584 
585  // Well Self
586  switch (wl.OptMode()) {
587  case RATE_MODE:
588  case ORATE_MODE:
589  case GRATE_MODE:
590  case WRATE_MODE:
591  case LRATE_MODE:
592  ls.AddRhs(wId, dt * wl.MaxRate());
593  break;
594  case BHP_MODE:
595  ls.AddDiag(wId, dt);
596  ls.AddRhs(wId, dt * wl.MaxBHP());
597  ls.AssignGuess(wId, wl.MaxBHP());
598  break;
599  default:
600  OCP_ABORT("Wrong well option mode!");
601  }
602 }
603 
604 void IsoT_IMPEC::AssembleMatProdWells(LinearSystem& ls,
605  const Bulk& bk,
606  const Well& wl,
607  const OCP_DBL& dt) const
608 {
609 
610  const OCP_USI wId = ls.AddDim(1) - 1;
611  ls.NewDiag(wId, 0.0);
612 
613  // Set Prod Weight
614  if (wl.OptMode() != BHP_MODE) wl.CalProdWeight(bk);
615 
616  const USI np = bk.numPhase;
617  const USI nc = bk.numCom;
618 
619  for (USI p = 0; p < wl.PerfNum(); p++) {
620  const OCP_USI n = wl.PerfLocation(p);
621 
622  OCP_DBL valb = 0;
623  OCP_DBL bb = 0;
624  OCP_DBL valw = 0;
625  OCP_DBL bw = 0;
626 
627  for (USI j = 0; j < np; j++) {
628  if (!bk.phaseExist[n * np + j]) continue;
629 
630  OCP_DBL tempb = 0;
631  OCP_DBL tempw = 0;
632 
633  for (USI i = 0; i < nc; i++) {
634  tempb += bk.vfi[n * nc + i] * bk.xij[n * np * nc + j * nc + i];
635  tempw += wl.ProdWeight(i) * bk.xij[n * np * nc + j * nc + i];
636  }
637  OCP_DBL trans = dt * wl.PerfTransj(p, j) * bk.xi[n * np + j];
638  valb += tempb * trans;
639  valw += tempw * trans;
640 
641  OCP_DBL dP = wl.DG(p) - bk.Pc[n * np + j];
642  bb += tempb * trans * dP;
643  bw += tempw * trans * dP;
644  }
645 
646  // Bulk to Well
647  ls.AddDiag(n, valb);
648  ls.NewOffDiag(n, wId, -valb);
649  ls.AddRhs(n, bb);
650 
651  // Well to Bulk
652  switch (wl.OptMode()) {
653  case RATE_MODE:
654  case ORATE_MODE:
655  case GRATE_MODE:
656  case WRATE_MODE:
657  case LRATE_MODE:
658  ls.AddDiag(wId, -valw);
659  ls.NewOffDiag(wId, n, valw);
660  ls.AddRhs(wId, bw);
661  break;
662  case BHP_MODE:
663  ls.NewOffDiag(wId, n, 0.0);
664  break;
665  default:
666  OCP_ABORT("Wrong well option mode!");
667  }
668  }
669 
670  // Well Self
671  switch (wl.OptMode()) {
672  case RATE_MODE:
673  case ORATE_MODE:
674  case GRATE_MODE:
675  case WRATE_MODE:
676  case LRATE_MODE:
677  ls.AddRhs(wId, dt * wl.MaxRate());
678  break;
679  case BHP_MODE:
680  ls.AddDiag(wId, dt);
681  ls.AddRhs(wId, dt * wl.MinBHP());
682  ls.AssignGuess(wId, wl.MinBHP());
683  break;
684  default:
685  OCP_ABORT("Wrong well option mode!");
686  }
687 }
688 
689 void IsoT_IMPEC::GetSolution(Reservoir& rs, const vector<OCP_DBL>& u)
690 {
691  Bulk& bk = rs.bulk;
692  const OCP_USI nb = bk.numBulk;
693  const USI np = bk.numPhase;
694 
695  // Bulk
696  for (OCP_USI n = 0; n < nb; n++) {
697  bk.P[n] = u[n];
698  for (USI j = 0; j < np; j++) {
699  bk.Pj[n * np + j] = bk.P[n] + bk.Pc[n * np + j];
700  }
701  }
702 
703  // Well
704  USI wId = nb;
705  for (auto& wl : rs.allWells.wells) {
706  if (wl.IsOpen()) {
707  wl.SetBHP(u[wId]);
708  wl.CalPerfP();
709  wId++;
710  }
711  }
712 }
713 
714 void IsoT_IMPEC::ResetToLastTimeStep01(Reservoir& rs, OCPControl& ctrl)
715 {
716  // Bulk
717  rs.bulk.Pj = rs.bulk.lPj;
718  // Bulk Conn
719  rs.conn.upblock = rs.conn.lupblock;
720  rs.conn.upblock_Rho = rs.conn.lupblock_Rho;
721  rs.conn.upblock_Trans = rs.conn.lupblock_Trans;
722  rs.conn.upblock_Velocity = rs.conn.lupblock_Velocity;
723 
724  // Iters
725  ctrl.ResetIterNRLS();
726 }
727 
728 void IsoT_IMPEC::ResetToLastTimeStep02(Reservoir& rs, OCPControl& ctrl)
729 {
730  // Bulk
731  rs.bulk.Ni = rs.bulk.lNi;
732  rs.bulk.Pj = rs.bulk.lPj;
733  // Bulk Conn
734  rs.conn.upblock = rs.conn.lupblock;
735  rs.conn.upblock_Rho = rs.conn.lupblock_Rho;
736  rs.conn.upblock_Trans = rs.conn.lupblock_Trans;
737  rs.conn.upblock_Velocity = rs.conn.lupblock_Velocity;
738 
739  // Iters
740  ctrl.ResetIterNRLS();
741 }
742 
743 void IsoT_IMPEC::ResetToLastTimeStep03(Reservoir& rs, OCPControl& ctrl)
744 {
745  Bulk& bk = rs.bulk;
746  // Rock
747  bk.rockVp = bk.lrockVp;
748  bk.poro = bk.lporo;
749  bk.poroP = bk.lporoP;
750 
751  // Fluid
752  bk.phaseNum = bk.lphaseNum;
753  bk.Nt = bk.lNt;
754  bk.Ni = bk.lNi;
755  bk.vf = bk.lvf;
756  bk.Pj = bk.lPj;
757  bk.phaseExist = bk.lphaseExist;
758  bk.S = bk.lS;
759  bk.vj = bk.lvj;
760  bk.xij = bk.lxij;
761  bk.rho = bk.lrho;
762  bk.xi = bk.lxi;
763  bk.mu = bk.lmu;
764 
765  // derivatives
766  bk.vfP = bk.lvfP;
767  bk.vfi = bk.lvfi;
768 
769  // Bulk Conn
770  rs.conn.upblock = rs.conn.lupblock;
771  rs.conn.upblock_Rho = rs.conn.lupblock_Rho;
772  rs.conn.upblock_Trans = rs.conn.lupblock_Trans;
773  rs.conn.upblock_Velocity = rs.conn.lupblock_Velocity;
774 
775  // Optional Features
776  rs.optFeatures.ResetToLastTimeStep();
777 
778  // Iters
779  ctrl.ResetIterNRLS();
780 }
781 
782 void IsoT_IMPEC::UpdateLastTimeStep(Reservoir& rs) const
783 {
784 
785  Bulk& bk = rs.bulk;
786 
787  // Rock
788  bk.lporo = bk.poro;
789  bk.lporoP = bk.poroP;
790  bk.lrockVp = bk.rockVp;
791 
792  // Fluid
793  bk.lphaseNum = bk.phaseNum;
794  bk.lNt = bk.Nt;
795  bk.lNi = bk.Ni;
796  bk.lvf = bk.vf;
797  bk.lP = bk.P;
798  bk.lPj = bk.Pj;
799  bk.lPc = bk.Pc;
800  bk.lphaseExist = bk.phaseExist;
801  bk.lS = bk.S;
802  bk.lvj = bk.vj;
803  bk.lxij = bk.xij;
804  bk.lrho = bk.rho;
805  bk.lxi = bk.xi;
806  bk.lmu = bk.mu;
807  bk.lkr = bk.kr;
808 
809  // derivatives
810  bk.lvfP = bk.vfP;
811  bk.lvfi = bk.vfi;
812 
813  BulkConn& conn = rs.conn;
814 
815  conn.lupblock = conn.upblock;
816  conn.lupblock_Rho = conn.upblock_Rho;
817  conn.lupblock_Trans = conn.upblock_Trans;
818  conn.lupblock_Velocity = conn.upblock_Velocity;
819 
820  rs.allWells.UpdateLastTimeStepBHP();
821  rs.optFeatures.UpdateLastTimeStep();
822 }
823 
825 // IsoT_FIM
827 
828 void IsoT_FIM::Setup(Reservoir& rs, LinearSystem& ls, const OCPControl& ctrl)
829 {
830  // Allocate memory for reservoir
831  AllocateReservoir(rs);
832  // Allocate memory for linear system
833  AllocateLinearSystem(ls, rs, ctrl);
834 }
835 
836 void IsoT_FIM::InitReservoir(Reservoir& rs) const
837 {
838  // Calculate initial bulk pressure and temperature and water saturation
839  rs.bulk.InitPTSw(50);
840  // Initialize rock property
841  InitRock(rs.bulk);
842  CalRock(rs.bulk);
843  // Initialize fluid properties
844  InitFlash(rs.bulk);
845  CalKrPc(rs.bulk);
846  // Initialize well pressure
847  rs.allWells.InitBHP(rs.bulk);
848  // Update variables at last time step
849  UpdateLastTimeStep(rs);
850 }
851 
852 void IsoT_FIM::Prepare(Reservoir& rs, const OCP_DBL& dt)
853 {
854  // Calculate well property at the beginning of next time step
855  rs.allWells.PrepareWell(rs.bulk);
856  // Calculate initial residual
857  CalRes(rs, dt, OCP_TRUE);
858 }
859 
860 void IsoT_FIM::AssembleMat(LinearSystem& ls,
861  const Reservoir& rs,
862  const OCP_DBL& dt) const
863 {
864  // Assemble matrix
865 #ifdef OCP_OLD_FIM
866  AssembleMatBulks(ls, rs, dt);
867  AssembleMatWells(ls, rs, dt);
868  // rs.allWells.AssemblaMatFIM(ls, rs.bulk, dt);
869 #else
870  AssembleMatBulksNew(ls, rs, dt);
871  AssembleMatWellsNew(ls, rs, dt);
872 #endif // OCP_OLD_FIM
873  // Assemble rhs -- from residual
874  ls.AssembleRhsCopy(rs.bulk.res.resAbs);
875 }
876 
877 void IsoT_FIM::SolveLinearSystem(LinearSystem& ls,
878  Reservoir& rs,
879  OCPControl& ctrl) const
880 {
881 #ifdef DEBUG
882  // Check if inf or nan occurs in A and b
883  ls.CheckEquation();
884 #endif // DEBUG
885 
886  // Assemble external linear solver with internal A and b
887  ls.AssembleMatLinearSolver();
888  // Solve linear system
889  GetWallTime Timer;
890  Timer.Start();
891  int status = ls.Solve();
892  if (status < 0) {
893  status = ls.GetNumIters();
894  }
895  // Record time, iterations
896  ctrl.RecordTimeLS(Timer.Stop() / 1000);
897  ctrl.UpdateIterLS(status);
898  ctrl.UpdateIterNR();
899 
900 #ifdef DEBUG
901  // Output A, b, x
902  // ls.OutputLinearSystem("testA_FIM.out", "testb_FIM.out");
903  // ls.OutputSolution("testx_FIM.out");
904  // Check if inf or nan occurs in solution
905  ls.CheckSolution();
906 #endif // DEBUG
907 
908  // Get solution from linear system to Reservoir
909  GetSolution(rs, ls.GetSolution(), ctrl);
910  // rs.PrintSolFIM(ctrl.workDir + "testPNi.out");
911  ls.ClearData();
912 }
913 
914 OCP_BOOL IsoT_FIM::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
915 {
916  OCP_DBL& dt = ctrl.current_dt;
917 
918  if (!ctrl.Check(rs, {"BulkNi", "BulkP"})) {
919  ResetToLastTimeStep(rs, ctrl);
920  cout << "Cut time step size and repeat! current dt = " << fixed
921  << setprecision(3) << dt << " days\n";
922  return OCP_FALSE;
923  }
924 
925  // Update fluid property
926  CalFlash(rs.bulk);
927  CalKrPc(rs.bulk);
928  // Update rock property
929  CalRock(rs.bulk);
930  // Update well property
931  rs.allWells.CalTrans(rs.bulk);
932  rs.allWells.CalFlux(rs.bulk);
933  // Update residual
934  CalRes(rs, dt, OCP_FALSE);
935 
936  return OCP_TRUE;
937 }
938 
939 OCP_BOOL IsoT_FIM::FinishNR(Reservoir& rs, OCPControl& ctrl)
940 {
941  OCP_USI dSn;
942 
943  const OCP_DBL NRdSmax = rs.GetNRdSmax(dSn);
944  const OCP_DBL NRdPmax = rs.GetNRdPmax();
945  // const OCP_DBL NRdNmax = rs.GetNRdNmax();
946 
947  if (((rs.bulk.res.maxRelRes_V <= rs.bulk.res.maxRelRes0_V * ctrl.ctrlNR.NRtol ||
948  rs.bulk.res.maxRelRes_V <= ctrl.ctrlNR.NRtol ||
949  rs.bulk.res.maxRelRes_N <= ctrl.ctrlNR.NRtol) &&
950  rs.bulk.res.maxWellRelRes_mol <= ctrl.ctrlNR.NRtol) ||
951  (fabs(NRdPmax) <= ctrl.ctrlNR.NRdPmin &&
952  fabs(NRdSmax) <= ctrl.ctrlNR.NRdSmin)) {
953 
954  if (!ctrl.Check(rs, {"WellP"})) {
955  ResetToLastTimeStep(rs, ctrl);
956  return OCP_FALSE;
957  } else {
958  return OCP_TRUE;
959  }
960  } else if (ctrl.iterNR >= ctrl.ctrlNR.maxNRiter) {
961  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
962  ResetToLastTimeStep(rs, ctrl);
963  cout << "### WARNING: NR not fully converged! Cut time step size and repeat! "
964  "current dt = "
965  << fixed << setprecision(3) << ctrl.current_dt << " days\n";
966  return OCP_FALSE;
967  } else {
968  return OCP_FALSE;
969  }
970 }
971 
972 void IsoT_FIM::FinishStep(Reservoir& rs, OCPControl& ctrl)
973 {
974  rs.CalIPRT(ctrl.GetCurDt());
975  rs.CalMaxChange();
976  UpdateLastTimeStep(rs);
977  ctrl.CalNextTimeStep(rs, {"dP", "dS", "iter"});
978 }
979 
980 void IsoT_FIM::AllocateReservoir(Reservoir& rs)
981 {
982  Bulk& bk = rs.bulk;
983  const OCP_USI nb = bk.numBulk;
984  const USI np = bk.numPhase;
985  const USI nc = bk.numCom;
986 
987  // Rock
988  bk.poro.resize(nb);
989  bk.rockVp.resize(nb);
990 
991  bk.lporo.resize(nb);
992  bk.lrockVp.resize(nb);
993 
994  // derivatives
995  bk.poroP.resize(nb);
996  bk.lporoP.resize(nb);
997 
998  // Fluid
999  bk.phaseNum.resize(nb);
1000  bk.Nt.resize(nb);
1001  bk.Ni.resize(nb * nc);
1002  bk.vf.resize(nb);
1003  bk.T.resize(nb);
1004  bk.P.resize(nb);
1005  bk.Pb.resize(nb);
1006  bk.Pj.resize(nb * np);
1007  bk.Pc.resize(nb * np);
1008  bk.phaseExist.resize(nb * np);
1009  bk.S.resize(nb * np);
1010  bk.xij.resize(nb * np * nc);
1011  bk.rho.resize(nb * np);
1012  bk.xi.resize(nb * np);
1013  bk.mu.resize(nb * np);
1014  bk.kr.resize(nb * np);
1015 
1016  bk.lphaseNum.resize(nb);
1017  bk.lNt.resize(nb);
1018  bk.lNi.resize(nb * nc);
1019  bk.lvf.resize(nb);
1020  bk.lT.resize(nb);
1021  bk.lP.resize(nb);
1022  bk.lPj.resize(nb * np);
1023  bk.lPc.resize(nb * np);
1024  bk.lphaseExist.resize(nb * np);
1025  bk.lS.resize(nb * np);
1026  bk.lxij.resize(nb * np * nc);
1027  bk.lrho.resize(nb * np);
1028  bk.lxi.resize(nb * np);
1029  bk.lmu.resize(nb * np);
1030  bk.lkr.resize(nb * np);
1031 
1032  // derivatives
1033  bk.vfP.resize(nb);
1034  bk.vfi.resize(nb * nc);
1035  bk.rhoP.resize(nb * np);
1036  bk.rhox.resize(nb * nc * np);
1037  bk.xiP.resize(nb * np);
1038  bk.xix.resize(nb * nc * np);
1039  bk.muP.resize(nb * np);
1040  bk.mux.resize(nb * nc * np);
1041  bk.dPcj_dS.resize(nb * np * np);
1042  bk.dKr_dS.resize(nb * np * np);
1043 
1044  bk.lvfP.resize(nb);
1045  bk.lvfi.resize(nb * nc);
1046  bk.lrhoP.resize(nb * np);
1047  bk.lrhox.resize(nb * nc * np);
1048  bk.lxiP.resize(nb * np);
1049  bk.lxix.resize(nb * nc * np);
1050  bk.lmuP.resize(nb * np);
1051  bk.lmux.resize(nb * nc * np);
1052  bk.ldPcj_dS.resize(nb * np * np);
1053  bk.ldKr_dS.resize(nb * np * np);
1054 
1055  // FIM-Specified
1056  bk.maxLendSdP = (nc + 1) * (nc + 1) * np;
1057  bk.dSec_dPri.resize(nb * bk.maxLendSdP);
1058  bk.bRowSizedSdP.resize(nb);
1059  bk.pSderExist.resize(nb * np);
1060  bk.pVnumCom.resize(nb * np);
1061 
1062  bk.ldSec_dPri.resize(nb * bk.maxLendSdP);
1063  bk.lbRowSizedSdP.resize(nb);
1064  bk.lpSderExist.resize(nb * np);
1065  bk.lpVnumCom.resize(nb * np);
1066 
1067  // Allocate Residual
1068  bk.res.Setup_IsoT(nb, rs.allWells.numWell, nc);
1069 
1070  // NR
1071  bk.NRstep.resize(nb);
1072  bk.NRphaseNum.resize(nb);
1073  bk.dSNR.resize(nb * np);
1074  bk.dSNRP.resize(nb * np);
1075  bk.dNNR.resize(nb * nc);
1076  bk.dPNR.resize(nb);
1077 
1078  // BulkConn
1079  BulkConn& conn = rs.conn;
1080 
1081  conn.upblock.resize(conn.numConn * np);
1082  conn.upblock_Rho.resize(conn.numConn * np);
1083  conn.upblock_Velocity.resize(conn.numConn * np);
1084 }
1085 
1086 void IsoT_FIM::AllocateLinearSystem(LinearSystem& ls,
1087  const Reservoir& rs,
1088  const OCPControl& ctrl)
1089 {
1090  ls.AllocateRowMem(rs.GetBulkNum() + rs.GetWellNum(), rs.GetComNum() + 1);
1091  ls.AllocateColMem(rs.conn.GetNeighborNum(), rs.allWells.GetWell2Bulk());
1092  ls.SetupLinearSolver(VECTORFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
1093 }
1094 
1095 void IsoT_FIM::InitFlash(Bulk& bk) const
1096 {
1097  for (OCP_USI n = 0; n < bk.numBulk; n++) {
1098  bk.flashCal[bk.PVTNUM[n]]->InitFlashFIM(bk.P[n], bk.Pb[n], bk.T[n],
1099  &bk.S[n * bk.numPhase], bk.rockVp[n],
1100  bk.Ni.data() + n * bk.numCom, n);
1101  for (USI i = 0; i < bk.numCom; i++) {
1102  bk.Ni[n * bk.numCom + i] = bk.flashCal[bk.PVTNUM[n]]->GetNi(i);
1103  }
1104  PassFlashValue(bk, n);
1105  }
1106 }
1107 
1108 void IsoT_FIM::CalFlash(Bulk& bk)
1109 {
1110  bk.maxNRdSSP = 0;
1111  bk.index_maxNRdSSP = 0;
1112 
1113  for (OCP_USI n = 0; n < bk.numBulk; n++) {
1114 
1115  bk.flashCal[bk.PVTNUM[n]]->FlashFIM(bk.P[n], bk.T[n], &bk.Ni[n * bk.numCom],
1116  &bk.S[n * bk.numPhase], bk.phaseNum[n],
1117  &bk.xij[n * bk.numPhase * bk.numCom], n);
1118  PassFlashValue(bk, n);
1119  }
1120 }
1121 
1122 void IsoT_FIM::PassFlashValue(Bulk& bk, const OCP_USI& n) const
1123 {
1124  const USI np = bk.numPhase;
1125  const USI nc = bk.numCom;
1126  const OCP_USI bIdp = n * np;
1127  const USI pvtnum = bk.PVTNUM[n];
1128  USI len = 0;
1129 
1130  bk.phaseNum[n] = 0;
1131  bk.Nt[n] = bk.flashCal[pvtnum]->GetNt();
1132  bk.vf[n] = bk.flashCal[pvtnum]->GetVf();
1133 
1134  for (USI j = 0; j < np; j++) {
1135  // Important! Saturation must be passed no matter if the phase exists. This is
1136  // because it will be used to calculate relative permeability and capillary
1137  // pressure at each time step. Make sure that all saturations are updated at
1138  // each step!
1139  bk.S[bIdp + j] = bk.flashCal[pvtnum]->GetS(j);
1140  bk.dSNR[bIdp + j] = bk.S[bIdp + j] - bk.dSNR[bIdp + j];
1141  if (bk.phaseExist[bIdp + j]) {
1142  if (fabs(bk.maxNRdSSP) < fabs(bk.dSNR[bIdp + j] - bk.dSNRP[bIdp + j])) {
1143  bk.maxNRdSSP = bk.dSNR[bIdp + j] - bk.dSNRP[bIdp + j];
1144  bk.index_maxNRdSSP = n;
1145  }
1146  }
1147 
1148  bk.phaseExist[bIdp + j] = bk.flashCal[pvtnum]->GetPhaseExist(j);
1149  if (bk.phaseExist[bIdp + j]) {
1150  bk.phaseNum[n]++;
1151  bk.rho[bIdp + j] = bk.flashCal[pvtnum]->GetRho(j);
1152  bk.xi[bIdp + j] = bk.flashCal[pvtnum]->GetXi(j);
1153  bk.mu[bIdp + j] = bk.flashCal[pvtnum]->GetMu(j);
1154 
1155  // Derivatives
1156  bk.rhoP[bIdp + j] = bk.flashCal[pvtnum]->GetRhoP(j);
1157  bk.xiP[bIdp + j] = bk.flashCal[pvtnum]->GetXiP(j);
1158  bk.muP[bIdp + j] = bk.flashCal[pvtnum]->GetMuP(j);
1159 
1160  for (USI i = 0; i < nc; i++) {
1161  bk.xij[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetXij(j, i);
1162  bk.rhox[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetRhoX(j, i);
1163  bk.xix[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetXiX(j, i);
1164  bk.mux[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetMuX(j, i);
1165  }
1166  }
1167 
1168  bk.pSderExist[bIdp + j] = bk.flashCal[pvtnum]->GetPSderExist(j);
1169  bk.pVnumCom[bIdp + j] = bk.flashCal[pvtnum]->GetPVnumCom(j);
1170  if (bk.pSderExist[bIdp + j]) len++;
1171  len += bk.pVnumCom[bIdp + j];
1172  }
1173 
1174  bk.vfP[n] = bk.flashCal[pvtnum]->GetVfP();
1175  for (USI i = 0; i < nc; i++) {
1176  bk.vfi[n * nc + i] = bk.flashCal[pvtnum]->GetVfi(i);
1177  }
1178 
1179 #ifdef OCP_OLD_FIM
1180  Dcopy(bk.maxLendSdP, &bk.dSec_dPri[n * bk.maxLendSdP],
1181  &bk.flashCal[pvtnum]->GetDXsDXp()[0]);
1182 #else
1183  bk.bRowSizedSdP[n] = len;
1184  len *= (nc + 1);
1185  Dcopy(len, &bk.dSec_dPri[n * bk.maxLendSdP], &bk.flashCal[pvtnum]->GetDXsDXp()[0]);
1186 #endif // OCP_OLD_FIM
1187 }
1188 
1189 void IsoT_FIM::CalKrPc(Bulk& bk) const
1190 {
1191  const USI& np = bk.numPhase;
1192  for (OCP_USI n = 0; n < bk.numBulk; n++) {
1193  const OCP_USI bId = n * np;
1194  bk.flow[bk.SATNUM[n]]->CalKrPcDeriv(&bk.S[bId], &bk.kr[bId], &bk.Pc[bId],
1195  &bk.dKr_dS[bId * np], &bk.dPcj_dS[bId * np],
1196  n);
1197  for (USI j = 0; j < np; j++) bk.Pj[bId + j] = bk.P[n] + bk.Pc[bId + j];
1198  }
1199 }
1200 
1201 void IsoT_FIM::CalRes(Reservoir& rs, const OCP_DBL& dt, const OCP_BOOL& resetRes0) const
1202 {
1203  const Bulk& bk = rs.bulk;
1204  const USI nb = bk.numBulk;
1205  const USI np = bk.numPhase;
1206  const USI nc = bk.numCom;
1207  const USI len = nc + 1;
1208  OCPRes& Res = bk.res;
1209  BulkConn& conn = rs.conn;
1210 
1211  Res.SetZero();
1212 
1213  // Bulk to Bulk
1214 
1215  OCP_USI bId, eId, uId, bIdb;
1216  // Accumalation Term
1217  for (OCP_USI n = 0; n < nb; n++) {
1218  bId = n * len;
1219  bIdb = n * nc;
1220  Res.resAbs[bId] = bk.rockVp[n] - bk.vf[n];
1221  for (USI i = 0; i < nc; i++) {
1222  Res.resAbs[bId + 1 + i] = bk.Ni[bIdb + i] - bk.lNi[bIdb + i];
1223  }
1224  }
1225 
1226  OCP_USI bId_np_j, eId_np_j, uId_np_j;
1227  OCP_DBL rho, dP, dNi, Akd;
1228 
1229  // Flux Term
1230  // Calculate the upblock at the same time.
1231  for (OCP_USI c = 0; c < conn.numConn; c++) {
1232  bId = conn.iteratorConn[c].BId();
1233  eId = conn.iteratorConn[c].EId();
1234  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
1235 
1236  for (USI j = 0; j < np; j++) {
1237  bId_np_j = bId * np + j;
1238  eId_np_j = eId * np + j;
1239 
1240  OCP_BOOL exbegin = bk.phaseExist[bId_np_j];
1241  OCP_BOOL exend = bk.phaseExist[eId_np_j];
1242 
1243  if ((exbegin) && (exend)) {
1244  rho = (bk.rho[bId_np_j] + bk.rho[eId_np_j]) / 2;
1245  } else if (exbegin && (!exend)) {
1246  rho = bk.rho[bId_np_j];
1247  } else if ((!exbegin) && (exend)) {
1248  rho = bk.rho[eId_np_j];
1249  } else {
1250  conn.upblock[c * np + j] = bId;
1251  conn.upblock_Rho[c * np + j] = 0;
1252  continue;
1253  }
1254 
1255  uId = bId;
1256  dP = (bk.Pj[bId_np_j] - GRAVITY_FACTOR * rho * bk.depth[bId]) -
1257  (bk.Pj[eId_np_j] - GRAVITY_FACTOR * rho * bk.depth[eId]);
1258  if (dP < 0) {
1259  uId = eId;
1260  }
1261  conn.upblock_Rho[c * np + j] = rho;
1262  conn.upblock[c * np + j] = uId;
1263  uId_np_j = uId * np + j;
1264 
1265  if (bk.phaseExist[uId_np_j]) {
1266  conn.upblock_Velocity[c * np + j] =
1267  Akd * bk.kr[uId_np_j] / bk.mu[uId_np_j] * dP;
1268  } else {
1269  conn.upblock_Velocity[c * np + j] = 0;
1270  continue;
1271  }
1272 
1273  OCP_DBL tmp =
1274  dt * Akd * bk.xi[uId_np_j] * bk.kr[uId_np_j] / bk.mu[uId_np_j] * dP;
1275 
1276  for (USI i = 0; i < nc; i++) {
1277  dNi = tmp * bk.xij[uId_np_j * nc + i];
1278  Res.resAbs[bId * len + 1 + i] += dNi;
1279  Res.resAbs[eId * len + 1 + i] -= dNi;
1280  }
1281  }
1282  }
1283 
1284  // Well to Bulk
1285  USI wId = nb * len;
1286  for (const auto& wl : rs.allWells.wells) {
1287  if (wl.IsOpen()) {
1288  // Well to Bulk
1289  for (USI p = 0; p < wl.PerfNum(); p++) {
1290  const OCP_USI k = wl.PerfLocation(p);
1291  for (USI i = 0; i < nc; i++) {
1292  Res.resAbs[k * len + 1 + i] += wl.PerfQi_lbmol(p, i) * dt;
1293  }
1294  }
1295  // Well Self
1296  if (wl.WellType() == INJ) {
1297  // Injection
1298  switch (wl.OptMode()) {
1299  case BHP_MODE:
1300  // bhp = opt.maxBHP;
1301  // Res.resAbs[bId] = bhp - opt.maxBHP;
1302  Res.resAbs[wId] = wl.BHP() - wl.MaxBHP();
1303  break;
1304  case RATE_MODE:
1305  case ORATE_MODE:
1306  case GRATE_MODE:
1307  case WRATE_MODE:
1308  case LRATE_MODE:
1309  Res.resAbs[wId] = wl.MaxRate();
1310  for (USI i = 0; i < nc; i++) {
1311  Res.resAbs[wId] += wl.Qi_lbmol(i);
1312  }
1313  // if (opt.reInj) {
1314  // for (auto& w : opt.connWell) {
1315  // OCP_DBL tmp = 0;
1316  // for (USI i = 0; i < numCom; i++) {
1317  // tmp += allWell[w].qi_lbmol[i];
1318  // }
1319  // tmp *= opt.reInjFactor;
1320  // Res.resAbs[bId] += tmp;
1321  // }
1322  // }
1323  Res.maxWellRelRes_mol =
1324  max(Res.maxWellRelRes_mol,
1325  fabs(Res.resAbs[wId] / wl.MaxRate()));
1326  break;
1327  default:
1328  OCP_ABORT("Wrong well opt mode!");
1329  break;
1330  }
1331  } else {
1332  // Production
1333  switch (wl.OptMode()) {
1334  case BHP_MODE:
1335  // bhp = opt.minBHP;
1336  // Res.resAbs[bId] = bhp - opt.minBHP;
1337  Res.resAbs[wId] = wl.BHP() - wl.MinBHP();
1338  break;
1339  case RATE_MODE:
1340  case ORATE_MODE:
1341  case GRATE_MODE:
1342  case WRATE_MODE:
1343  case LRATE_MODE:
1344  wl.CalProdWeight(bk);
1345  Res.resAbs[wId] = -wl.MaxRate();
1346  for (USI i = 0; i < nc; i++) {
1347  Res.resAbs[wId] += wl.Qi_lbmol(i) * wl.ProdWeight(i);
1348  }
1349  Res.maxWellRelRes_mol =
1350  max(Res.maxWellRelRes_mol,
1351  fabs(Res.resAbs[wId] / wl.MaxRate()));
1352  break;
1353  default:
1354  OCP_ABORT("Wrong well opt mode!");
1355  break;
1356  }
1357  }
1358  wId += len;
1359  }
1360  }
1361 
1362  // Calculate RelRes
1363  OCP_DBL tmp;
1364  for (OCP_USI n = 0; n < nb; n++) {
1365 
1366  for (USI i = 0; i < len; i++) {
1367  tmp = fabs(Res.resAbs[n * len + i] / bk.rockVp[n]);
1368  if (Res.maxRelRes_V < tmp) {
1369  Res.maxRelRes_V = tmp;
1370  Res.maxId_V = n;
1371  }
1372  Res.resRelV[n] += tmp * tmp;
1373  }
1374  Res.resRelV[n] = sqrt(Res.resRelV[n]);
1375 
1376  for (USI i = 1; i < len; i++) {
1377  tmp = fabs(Res.resAbs[n * len + i] / bk.Nt[n]);
1378  if (Res.maxRelRes_N < tmp) {
1379  Res.maxRelRes_N = tmp;
1380  Res.maxId_N = n;
1381  }
1382  Res.resRelN[n] += tmp * tmp;
1383  }
1384  Res.resRelN[n] = sqrt(Res.resRelN[n]);
1385  }
1386 
1387  Dscalar(Res.resAbs.size(), -1.0, Res.resAbs.data());
1388  if (resetRes0) Res.SetInitRes();
1389 }
1390 
1391 void IsoT_FIM::AssembleMatBulks(LinearSystem& ls,
1392  const Reservoir& rs,
1393  const OCP_DBL& dt) const
1394 {
1395  const Bulk& bk = rs.bulk;
1396  const BulkConn& conn = rs.conn;
1397  const OCP_USI nb = bk.numBulk;
1398  const USI np = bk.numPhase;
1399  const USI nc = bk.numCom;
1400  const USI ncol = nc + 1;
1401  const USI ncol2 = np * nc + np;
1402  const USI bsize = ncol * ncol;
1403  const USI bsize2 = ncol * ncol2;
1404 
1405  ls.AddDim(nb);
1406 
1407  vector<OCP_DBL> bmat(bsize, 0);
1408  // Accumulation term
1409  for (USI i = 1; i < ncol; i++) {
1410  bmat[i * ncol + i] = 1;
1411  }
1412  for (OCP_USI n = 0; n < nb; n++) {
1413  bmat[0] = bk.v[n] * bk.poroP[n] - bk.vfP[n];
1414  for (USI i = 0; i < nc; i++) {
1415  bmat[i + 1] = -bk.vfi[n * nc + i];
1416  }
1417  ls.NewDiag(n, bmat);
1418  }
1419 
1420  // flux term
1421  OCP_DBL Akd;
1422  OCP_DBL transJ, transIJ;
1423  vector<OCP_DBL> dFdXpB(bsize, 0); // begin bulk: dF / dXp
1424  vector<OCP_DBL> dFdXpE(bsize, 0); // end bulk: dF / dXp
1425  vector<OCP_DBL> dFdXsB(bsize2, 0); // begin bulk: dF / dXs
1426  vector<OCP_DBL> dFdXsE(bsize2, 0); // end bulk: dF / dXs
1427  OCP_DBL* dFdXpU; // up bulk: dF / dXp
1428  OCP_DBL* dFdXpD; // down bulk: dF / dXp
1429  OCP_DBL* dFdXsU; // up bulk: dF / dXs
1430  OCP_DBL* dFdXsD; // down bulk: dF / dXs
1431 
1432  OCP_USI bId, eId, uId;
1433  OCP_USI bId_np_j, eId_np_j, uId_np_j, dId_np_j;
1434  OCP_BOOL phaseExistBj, phaseExistEj, phaseExistDj;
1435  OCP_DBL kr, mu, xi, xij, xiP, muP, rhox, xix, mux;
1436  OCP_DBL dP, dGamma;
1437  OCP_DBL rhoWghtU, rhoWghtD;
1438  OCP_DBL tmp;
1439 
1440  for (OCP_USI c = 0; c < conn.numConn; c++) {
1441 
1442  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
1443  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
1444  fill(dFdXsB.begin(), dFdXsB.end(), 0.0);
1445  fill(dFdXsE.begin(), dFdXsE.end(), 0.0);
1446 
1447  bId = conn.iteratorConn[c].BId();
1448  eId = conn.iteratorConn[c].EId();
1449  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
1450  dGamma = GRAVITY_FACTOR * (bk.depth[bId] - bk.depth[eId]);
1451 
1452  for (USI j = 0; j < np; j++) {
1453  uId = conn.upblock[c * np + j];
1454  uId_np_j = uId * np + j;
1455  if (!bk.phaseExist[uId_np_j]) continue;
1456  bId_np_j = bId * np + j;
1457  eId_np_j = eId * np + j;
1458  phaseExistBj = bk.phaseExist[bId_np_j];
1459  phaseExistEj = bk.phaseExist[eId_np_j];
1460 
1461  if (bId == uId) {
1462  dFdXpU = &dFdXpB[0];
1463  dFdXpD = &dFdXpE[0];
1464  dFdXsU = &dFdXsB[0];
1465  dFdXsD = &dFdXsE[0];
1466  phaseExistDj = phaseExistEj;
1467  dId_np_j = eId_np_j;
1468  } else {
1469  dFdXpU = &dFdXpE[0];
1470  dFdXpD = &dFdXpB[0];
1471  dFdXsU = &dFdXsE[0];
1472  dFdXsD = &dFdXsB[0];
1473  phaseExistDj = phaseExistBj;
1474  dId_np_j = bId_np_j;
1475  }
1476  if (phaseExistDj) {
1477  rhoWghtU = 0.5;
1478  rhoWghtD = 0.5;
1479  } else {
1480  rhoWghtU = 1;
1481  rhoWghtD = 0;
1482  }
1483 
1484  dP = bk.Pj[bId_np_j] - bk.Pj[eId_np_j] -
1485  conn.upblock_Rho[c * np + j] * dGamma;
1486  xi = bk.xi[uId_np_j];
1487  kr = bk.kr[uId_np_j];
1488  mu = bk.mu[uId_np_j];
1489  muP = bk.muP[uId_np_j];
1490  xiP = bk.xiP[uId_np_j];
1491  transJ = Akd * kr / mu;
1492 
1493  for (USI i = 0; i < nc; i++) {
1494  xij = bk.xij[uId_np_j * nc + i];
1495  transIJ = xij * xi * transJ;
1496 
1497  // dP
1498  dFdXpB[(i + 1) * ncol] += transIJ;
1499  dFdXpE[(i + 1) * ncol] -= transIJ;
1500 
1501  tmp = transJ * xiP * xij * dP;
1502  tmp += -transIJ * muP / mu * dP;
1503  dFdXpU[(i + 1) * ncol] +=
1504  (tmp - transIJ * rhoWghtU * bk.rhoP[uId_np_j] * dGamma);
1505  dFdXpD[(i + 1) * ncol] +=
1506  -transIJ * rhoWghtD * bk.rhoP[dId_np_j] * dGamma;
1507 
1508  // dS
1509  for (USI k = 0; k < np; k++) {
1510  dFdXsB[(i + 1) * ncol2 + k] +=
1511  transIJ * bk.dPcj_dS[bId_np_j * np + k];
1512  dFdXsE[(i + 1) * ncol2 + k] -=
1513  transIJ * bk.dPcj_dS[eId_np_j * np + k];
1514  dFdXsU[(i + 1) * ncol2 + k] +=
1515  Akd * bk.dKr_dS[uId_np_j * np + k] / mu * xi * xij * dP;
1516  }
1517  // dxij
1518  for (USI k = 0; k < nc; k++) {
1519  rhox = bk.rhox[uId_np_j * nc + k];
1520  xix = bk.xix[uId_np_j * nc + k];
1521  mux = bk.mux[uId_np_j * nc + k];
1522  tmp = -transIJ * rhoWghtU * rhox * dGamma;
1523  tmp += transJ * xix * xij * dP;
1524  tmp += -transIJ * mux / mu * dP;
1525  dFdXsU[(i + 1) * ncol2 + np + j * nc + k] += tmp;
1526  dFdXsD[(i + 1) * ncol2 + np + j * nc + k] +=
1527  -transIJ * rhoWghtD * bk.rhox[dId_np_j * nc + k] * dGamma;
1528  }
1529  dFdXsU[(i + 1) * ncol2 + np + j * nc + i] += transJ * xi * dP;
1530  }
1531  }
1532 
1533  // Assemble
1534  bmat = dFdXpB;
1535  DaABpbC(ncol, ncol, ncol2, 1, dFdXsB.data(), &bk.dSec_dPri[bId * bsize2], 1,
1536  bmat.data());
1537  Dscalar(bsize, dt, bmat.data());
1538  // Begin - Begin -- add
1539  ls.AddDiag(bId, bmat);
1540  // End - Begin -- insert
1541  Dscalar(bsize, -1, bmat.data());
1542  ls.NewOffDiag(eId, bId, bmat);
1543 
1544 #ifdef OCP_NANCHECK
1545  if (!CheckNan(bmat.size(), &bmat[0])) {
1546  OCP_ABORT("INF or INF in bmat !");
1547  }
1548 #endif
1549 
1550  // End
1551  bmat = dFdXpE;
1552  DaABpbC(ncol, ncol, ncol2, 1, dFdXsE.data(), &bk.dSec_dPri[eId * bsize2], 1,
1553  bmat.data());
1554  Dscalar(bsize, dt, bmat.data());
1555  // Begin - End -- insert
1556  ls.NewOffDiag(bId, eId, bmat);
1557  // End - End -- add
1558  Dscalar(bsize, -1, bmat.data());
1559  ls.AddDiag(eId, bmat);
1560 
1561 #ifdef OCP_NANCHECK
1562  if (!CheckNan(bmat.size(), &bmat[0])) {
1563  OCP_ABORT("INF or INF in bmat !");
1564  }
1565 #endif
1566  }
1567 }
1568 
1569 void IsoT_FIM::AssembleMatBulksNew(LinearSystem& ls,
1570  const Reservoir& rs,
1571  const OCP_DBL& dt) const
1572 {
1573  const Bulk& bk = rs.bulk;
1574  const BulkConn& conn = rs.conn;
1575  const OCP_USI nb = bk.numBulk;
1576  const USI np = bk.numPhase;
1577  const USI nc = bk.numCom;
1578  const USI ncol = nc + 1;
1579  const USI ncol2 = np * nc + np;
1580  const USI bsize = ncol * ncol;
1581  const USI bsize2 = ncol * ncol2;
1582  const USI lendSdP = bk.maxLendSdP;
1583 
1584  ls.AddDim(nb);
1585 
1586  vector<OCP_DBL> bmat(bsize, 0);
1587  // Accumulation term
1588  for (USI i = 1; i < ncol; i++) {
1589  bmat[i * ncol + i] = 1;
1590  }
1591  for (OCP_USI n = 0; n < nb; n++) {
1592  bmat[0] = bk.v[n] * bk.poroP[n] - bk.vfP[n];
1593  for (USI i = 0; i < nc; i++) {
1594  bmat[i + 1] = -bk.vfi[n * nc + i];
1595  }
1596  ls.NewDiag(n, bmat);
1597  }
1598 
1599  // flux term
1600  OCP_DBL Akd;
1601  OCP_DBL transJ, transIJ;
1602  vector<OCP_DBL> dFdXpB(bsize, 0);
1603  vector<OCP_DBL> dFdXpE(bsize, 0);
1604  vector<OCP_DBL> dFdXsB(bsize2, 0);
1605  vector<OCP_DBL> dFdXsE(bsize2, 0);
1606  vector<OCP_BOOL> phaseExistB(np, OCP_FALSE);
1607  vector<OCP_BOOL> phaseExistE(np, OCP_FALSE);
1608  OCP_BOOL phaseExistU;
1609  vector<OCP_BOOL> phasedS_B(np, OCP_FALSE);
1610  vector<OCP_BOOL> phasedS_E(np, OCP_FALSE);
1611  vector<USI> pVnumComB(np, 0);
1612  vector<USI> pVnumComE(np, 0);
1613  USI ncolB, ncolE;
1614 
1615  OCP_USI bId, eId, uId;
1616  OCP_USI bId_np_j, eId_np_j, uId_np_j;
1617  OCP_DBL kr, mu, xi, xij, rhoP, xiP, muP, rhox, xix, mux;
1618  OCP_DBL dP, dGamma;
1619  OCP_DBL tmp;
1620 
1621  for (OCP_USI c = 0; c < conn.numConn; c++) {
1622  bId = conn.iteratorConn[c].BId();
1623  eId = conn.iteratorConn[c].EId();
1624  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
1625  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
1626  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
1627  fill(dFdXsB.begin(), dFdXsB.end(), 0.0);
1628  fill(dFdXsE.begin(), dFdXsE.end(), 0.0);
1629  dGamma = GRAVITY_FACTOR * (bk.depth[bId] - bk.depth[eId]);
1630 
1631  USI jxB = 0;
1632  USI jxE = 0;
1633  ncolB = 0;
1634  ncolE = 0;
1635 
1636  for (USI j = 0; j < np; j++) {
1637  phaseExistB[j] = bk.phaseExist[bId * np + j];
1638  phaseExistE[j] = bk.phaseExist[eId * np + j];
1639  phasedS_B[j] = bk.pSderExist[bId * np + j];
1640  phasedS_E[j] = bk.pSderExist[eId * np + j];
1641  if (phasedS_B[j]) jxB++;
1642  if (phasedS_E[j]) jxE++;
1643  pVnumComB[j] = bk.pVnumCom[bId * np + j];
1644  pVnumComE[j] = bk.pVnumCom[eId * np + j];
1645  ncolB += pVnumComB[j];
1646  ncolE += pVnumComE[j];
1647  }
1648  ncolB += jxB;
1649  ncolE += jxE;
1650 
1651  for (USI j = 0; j < np; j++) {
1652  uId = conn.upblock[c * np + j];
1653 
1654  phaseExistU = (uId == bId ? phaseExistB[j] : phaseExistE[j]);
1655  if (!phaseExistU) {
1656  jxB += pVnumComB[j];
1657  jxE += pVnumComE[j];
1658  continue;
1659  }
1660 
1661  bId_np_j = bId * np + j;
1662  eId_np_j = eId * np + j;
1663  uId_np_j = uId * np + j;
1664  dP = bk.Pj[bId_np_j] - bk.Pj[eId_np_j] -
1665  conn.upblock_Rho[c * np + j] * dGamma;
1666  xi = bk.xi[uId_np_j];
1667  kr = bk.kr[uId_np_j];
1668  mu = bk.mu[uId_np_j];
1669  muP = bk.muP[uId_np_j];
1670  xiP = bk.xiP[uId_np_j];
1671  rhoP = bk.rhoP[uId_np_j];
1672  transJ = Akd * kr / mu;
1673 
1674  for (USI i = 0; i < nc; i++) {
1675  xij = bk.xij[uId_np_j * nc + i];
1676  transIJ = xij * xi * transJ;
1677 
1678  // Pressure -- Primary var
1679  dFdXpB[(i + 1) * ncol] += transIJ;
1680  dFdXpE[(i + 1) * ncol] -= transIJ;
1681 
1682  tmp = xij * transJ * xiP * dP;
1683  tmp += -transIJ * muP / mu * dP;
1684  if (!phaseExistE[j]) {
1685  tmp += transIJ * (-rhoP * dGamma);
1686  dFdXpB[(i + 1) * ncol] += tmp;
1687  } else if (!phaseExistB[j]) {
1688  tmp += transIJ * (-rhoP * dGamma);
1689  dFdXpE[(i + 1) * ncol] += tmp;
1690  } else {
1691  dFdXpB[(i + 1) * ncol] +=
1692  transIJ * (-bk.rhoP[bId_np_j] * dGamma) / 2;
1693  dFdXpE[(i + 1) * ncol] +=
1694  transIJ * (-bk.rhoP[eId_np_j] * dGamma) / 2;
1695  if (bId == uId) {
1696  dFdXpB[(i + 1) * ncol] += tmp;
1697  } else {
1698  dFdXpE[(i + 1) * ncol] += tmp;
1699  }
1700  }
1701 
1702  // Second var
1703  USI j1SB = 0;
1704  USI j1SE = 0;
1705  if (bId == uId) {
1706  // Saturation
1707  for (USI j1 = 0; j1 < np; j1++) {
1708  if (phasedS_B[j1]) {
1709  dFdXsB[(i + 1) * ncolB + j1SB] +=
1710  transIJ * bk.dPcj_dS[bId_np_j * np + j1];
1711  tmp = Akd * xij * xi / mu * bk.dKr_dS[uId_np_j * np + j1] *
1712  dP;
1713  dFdXsB[(i + 1) * ncolB + j1SB] += tmp;
1714  j1SB++;
1715  }
1716  if (phasedS_E[j1]) {
1717  dFdXsE[(i + 1) * ncolE + j1SE] -=
1718  transIJ * bk.dPcj_dS[eId_np_j * np + j1];
1719  j1SE++;
1720  }
1721  }
1722  // Cij
1723  if (!phaseExistE[j]) {
1724  for (USI k = 0; k < pVnumComB[j]; k++) {
1725  rhox = bk.rhox[uId_np_j * nc + k];
1726  xix = bk.xix[uId_np_j * nc + k];
1727  mux = bk.mux[uId_np_j * nc + k];
1728  tmp = -transIJ * rhox * dGamma;
1729  tmp += xij * transJ * xix * dP;
1730  tmp += -transIJ * mux / mu * dP;
1731  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
1732  }
1733  // WARNING !!!
1734  if (i < pVnumComB[j])
1735  dFdXsB[(i + 1) * ncolB + jxB + i] += xi * transJ * dP;
1736  } else {
1737  for (USI k = 0; k < pVnumComB[j]; k++) {
1738  rhox = bk.rhox[bId_np_j * nc + k] / 2;
1739  xix = bk.xix[uId_np_j * nc + k];
1740  mux = bk.mux[uId_np_j * nc + k];
1741  tmp = -transIJ * rhox * dGamma;
1742  tmp += xij * transJ * xix * dP;
1743  tmp += -transIJ * mux / mu * dP;
1744  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
1745  dFdXsE[(i + 1) * ncolE + jxE + k] +=
1746  -transIJ * bk.rhox[eId_np_j * nc + k] / 2 * dGamma;
1747  }
1748  // WARNING !!!
1749  if (i < pVnumComB[j])
1750  dFdXsB[(i + 1) * ncolB + jxB + i] += xi * transJ * dP;
1751  }
1752  } else {
1753  // Saturation
1754  for (USI j1 = 0; j1 < np; j1++) {
1755  if (phasedS_B[j1]) {
1756  dFdXsB[(i + 1) * ncolB + j1SB] +=
1757  transIJ * bk.dPcj_dS[bId_np_j * np + j1];
1758  j1SB++;
1759  }
1760  if (phasedS_E[j1]) {
1761  dFdXsE[(i + 1) * ncolE + j1SE] -=
1762  transIJ * bk.dPcj_dS[eId_np_j * np + j1];
1763  tmp = Akd * xij * xi / mu * bk.dKr_dS[uId_np_j * np + j1] *
1764  dP;
1765  dFdXsE[(i + 1) * ncolE + j1SE] += tmp;
1766  j1SE++;
1767  }
1768  }
1769  // Cij
1770  if (!phaseExistB[j]) {
1771  for (USI k = 0; k < pVnumComE[j]; k++) {
1772  rhox = bk.rhox[uId_np_j * nc + k];
1773  xix = bk.xix[uId_np_j * nc + k];
1774  mux = bk.mux[uId_np_j * nc + k];
1775  tmp = -transIJ * rhox * dGamma;
1776  tmp += xij * transJ * xix * dP;
1777  tmp += -transIJ * mux / mu * dP;
1778  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
1779  }
1780  // WARNING !!!
1781  if (i < pVnumComE[j])
1782  dFdXsE[(i + 1) * ncolE + jxE + i] += xi * transJ * dP;
1783  } else {
1784  for (USI k = 0; k < pVnumComE[j]; k++) {
1785  rhox = bk.rhox[eId_np_j * nc + k] / 2;
1786  xix = bk.xix[uId_np_j * nc + k];
1787  mux = bk.mux[uId_np_j * nc + k];
1788  tmp = -transIJ * rhox * dGamma;
1789  tmp += xij * transJ * xix * dP;
1790  tmp += -transIJ * mux / mu * dP;
1791  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
1792  dFdXsB[(i + 1) * ncolB + jxB + k] +=
1793  -transIJ * bk.rhox[bId_np_j * nc + k] / 2 * dGamma;
1794  }
1795  // WARNING !!!
1796  if (i < pVnumComE[j])
1797  dFdXsE[(i + 1) * ncolE + jxE + i] += xi * transJ * dP;
1798  }
1799  }
1800  }
1801  jxB += pVnumComB[j];
1802  jxE += pVnumComE[j];
1803  }
1804 
1805  // Assemble
1806  bmat = dFdXpB;
1807  DaABpbC(ncol, ncol, ncolB, 1, dFdXsB.data(), &bk.dSec_dPri[bId * lendSdP], 1,
1808  bmat.data());
1809  Dscalar(bsize, dt, bmat.data());
1810  // Begin - Begin -- add
1811  ls.AddDiag(bId, bmat);
1812  // End - Begin -- insert
1813  Dscalar(bsize, -1, bmat.data());
1814  ls.NewOffDiag(eId, bId, bmat);
1815 
1816 #ifdef OCP_NANCHECK
1817  if (!CheckNan(bmat.size(), &bmat[0])) {
1818  OCP_ABORT("INF or NAN in bmat !");
1819  }
1820 #endif
1821 
1822  bmat = dFdXpE;
1823  DaABpbC(ncol, ncol, ncolE, 1, dFdXsE.data(), &bk.dSec_dPri[eId * lendSdP], 1,
1824  bmat.data());
1825 
1826  Dscalar(bsize, dt, bmat.data());
1827  // Begin - End -- insert
1828  ls.NewOffDiag(bId, eId, bmat);
1829  // End - End -- add
1830  Dscalar(bsize, -1, bmat.data());
1831  ls.AddDiag(eId, bmat);
1832 
1833 #ifdef OCP_NANCHECK
1834  if (!CheckNan(bmat.size(), &bmat[0])) {
1835  OCP_ABORT("INF or INF in bmat !");
1836  }
1837 #endif
1838  }
1839 }
1840 
1841 void IsoT_FIM::AssembleMatBulksNewS(LinearSystem& ls,
1842  const Reservoir& rs,
1843  const OCP_DBL& dt) const
1844 {
1845 
1846  const Bulk& bk = rs.bulk;
1847  const BulkConn& conn = rs.conn;
1848  const OCP_USI nb = bk.numCom;
1849  const USI np = bk.numPhase;
1850  const USI nc = bk.numCom;
1851  const USI ncol = nc + 1;
1852  const USI ncol2 = np * nc + np;
1853  const USI bsize = ncol * ncol;
1854  const USI bsize2 = ncol * ncol2;
1855  const USI lendSdP = bk.maxLendSdP;
1856 
1857  ls.AddDim(nb);
1858 
1859  vector<OCP_DBL> bmat(bsize, 0);
1860  // Accumulation term
1861  for (USI i = 1; i < ncol; i++) {
1862  bmat[i * ncol + i] = 1;
1863  }
1864  for (OCP_USI n = 0; n < nb; n++) {
1865  bmat[0] = bk.v[n] * bk.poroP[n] - bk.vfP[n];
1866  for (USI i = 0; i < nc; i++) {
1867  bmat[i + 1] = -bk.vfi[n * nc + i];
1868  }
1869  ls.NewDiag(n, bmat);
1870  }
1871 
1872  // flux term
1873  OCP_DBL Akd;
1874  OCP_DBL transJ, transIJ;
1875  vector<OCP_DBL> dFdXpB(bsize, 0);
1876  vector<OCP_DBL> dFdXpE(bsize, 0);
1877  vector<OCP_DBL> dFdXsB(bsize2, 0);
1878  vector<OCP_DBL> dFdXsE(bsize2, 0);
1879  vector<OCP_BOOL> phaseExistB(np, OCP_FALSE);
1880  vector<OCP_BOOL> phaseExistE(np, OCP_FALSE);
1881  OCP_BOOL phaseExistU;
1882  vector<USI> pEnumComB(np, 0);
1883  vector<USI> pEnumComE(np, 0);
1884  USI ncolB, ncolE;
1885 
1886  OCP_USI bId, eId, uId;
1887  OCP_USI bId_np_j, eId_np_j, uId_np_j;
1888  OCP_DBL kr, mu, xi, xij, rhoP, xiP, muP, rhox, xix, mux;
1889  OCP_DBL dP, dGamma;
1890  OCP_DBL tmp;
1891  OCP_DBL wghtb, wghte;
1892 
1893  for (OCP_USI c = 0; c < conn.numConn; c++) {
1894  bId = conn.iteratorConn[c].BId();
1895  eId = conn.iteratorConn[c].EId();
1896  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
1897  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
1898  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
1899  fill(dFdXsB.begin(), dFdXsB.end(), 0.0);
1900  fill(dFdXsE.begin(), dFdXsE.end(), 0.0);
1901  dGamma = GRAVITY_FACTOR * (bk.depth[bId] - bk.depth[eId]);
1902 
1903  const USI npB = bk.phaseNum[bId];
1904  ncolB = npB;
1905  const USI npE = bk.phaseNum[eId];
1906  ncolE = npE;
1907 
1908  for (USI j = 0; j < np; j++) {
1909  phaseExistB[j] = bk.phaseExist[bId * np + j];
1910  phaseExistE[j] = bk.phaseExist[eId * np + j];
1911  pEnumComB[j] = bk.pVnumCom[bId * np + j];
1912  pEnumComE[j] = bk.pVnumCom[eId * np + j];
1913  ncolB += pEnumComB[j];
1914  ncolE += pEnumComE[j];
1915  }
1916 
1917  USI jxB = npB;
1918  USI jxE = npE;
1919  for (USI j = 0; j < np; j++) {
1920  uId = conn.upblock[c * np + j];
1921 
1922  phaseExistU = (uId == bId ? phaseExistB[j] : phaseExistE[j]);
1923  if (!phaseExistU) {
1924  jxB += pEnumComB[j];
1925  jxE += pEnumComE[j];
1926  continue;
1927  }
1928 
1929  bId_np_j = bId * np + j;
1930  eId_np_j = eId * np + j;
1931  uId_np_j = uId * np + j;
1932  dP = bk.Pj[bId_np_j] - bk.Pj[eId_np_j] -
1933  conn.upblock_Rho[c * np + j] * dGamma;
1934  xi = bk.xi[uId_np_j];
1935  kr = bk.kr[uId_np_j];
1936  mu = bk.mu[uId_np_j];
1937  muP = bk.muP[uId_np_j];
1938  xiP = bk.xiP[uId_np_j];
1939  rhoP = bk.rhoP[uId_np_j];
1940  transJ = Akd * kr / mu;
1941 
1942  for (USI i = 0; i < nc; i++) {
1943  xij = bk.xij[uId_np_j * nc + i];
1944  transIJ = xij * xi * transJ;
1945 
1946  // Pressure -- Primary var
1947  dFdXpB[(i + 1) * ncol] += transIJ;
1948  dFdXpE[(i + 1) * ncol] -= transIJ;
1949 
1950  tmp = xij * transJ * xiP * dP;
1951  tmp += -transIJ * muP / mu * dP;
1952  if (!phaseExistE[j]) {
1953  tmp += transIJ * (-rhoP * dGamma);
1954  dFdXpB[(i + 1) * ncol] += tmp;
1955  } else if (!phaseExistB[j]) {
1956  tmp += transIJ * (-rhoP * dGamma);
1957  dFdXpE[(i + 1) * ncol] += tmp;
1958  } else {
1959  dFdXpB[(i + 1) * ncol] += transIJ * dGamma *
1960  (-bk.rhoP[bId_np_j] * bk.S[bId_np_j]) /
1961  (bk.S[bId_np_j] + bk.S[eId_np_j]);
1962  dFdXpE[(i + 1) * ncol] += transIJ * dGamma *
1963  (-bk.rhoP[eId_np_j] * bk.S[eId_np_j]) /
1964  (bk.S[bId_np_j] + bk.S[eId_np_j]);
1965  if (bId == uId) {
1966  dFdXpB[(i + 1) * ncol] += tmp;
1967  } else {
1968  dFdXpE[(i + 1) * ncol] += tmp;
1969  }
1970  }
1971 
1972  // Second var
1973  USI j1SB = 0;
1974  USI j1SE = 0;
1975  if (bId == uId) {
1976  // Saturation
1977  for (USI j1 = 0; j1 < np; j1++) {
1978 
1979  wghtb = 0;
1980  wghte = 0;
1981  if (j1 == j && phaseExistE[j]) {
1982  tmp = -dGamma / ((bk.S[bId_np_j] + bk.S[eId_np_j]) *
1983  (bk.S[bId_np_j] + bk.S[eId_np_j]));
1984  wghtb = tmp * bk.rho[bId_np_j] * bk.S[eId_np_j];
1985  wghte = tmp * bk.rho[eId_np_j] * bk.S[bId_np_j];
1986  }
1987 
1988  if (phaseExistB[j1]) {
1989  dFdXsB[(i + 1) * ncolB + j1SB] +=
1990  transIJ * (bk.dPcj_dS[bId_np_j * np + j1] + wghtb);
1991  tmp = Akd * xij * xi / mu * bk.dKr_dS[uId_np_j * np + j1] *
1992  dP;
1993  dFdXsB[(i + 1) * ncolB + j1SB] += tmp;
1994  j1SB++;
1995  }
1996  if (phaseExistE[j1]) {
1997  dFdXsE[(i + 1) * ncolE + j1SE] -=
1998  transIJ * (bk.dPcj_dS[eId_np_j * np + j1] + wghte);
1999  j1SE++;
2000  }
2001  }
2002  // Cij
2003  if (!phaseExistE[j]) {
2004  for (USI k = 0; k < pEnumComB[j]; k++) {
2005  rhox = bk.rhox[uId_np_j * nc + k];
2006  xix = bk.xix[uId_np_j * nc + k];
2007  mux = bk.mux[uId_np_j * nc + k];
2008  tmp = -transIJ * rhox * dGamma;
2009  tmp += xij * transJ * xix * dP;
2010  tmp += -transIJ * mux / mu * dP;
2011  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
2012  }
2013  // WARNING !!!
2014  if (i < pEnumComB[j])
2015  dFdXsB[(i + 1) * ncolB + jxB + i] += xi * transJ * dP;
2016  } else {
2017  wghtb = bk.S[bId_np_j] / (bk.S[bId_np_j] + bk.S[eId_np_j]);
2018  wghte = bk.S[eId_np_j] / (bk.S[bId_np_j] + bk.S[eId_np_j]);
2019  for (USI k = 0; k < pEnumComB[j]; k++) {
2020  rhox = bk.rhox[bId_np_j * nc + k] * wghtb;
2021  xix = bk.xix[uId_np_j * nc + k];
2022  mux = bk.mux[uId_np_j * nc + k];
2023  tmp = -transIJ * rhox * dGamma;
2024  tmp += xij * transJ * xix * dP;
2025  tmp += -transIJ * mux / mu * dP;
2026  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
2027  dFdXsE[(i + 1) * ncolE + jxE + k] +=
2028  -transIJ * dGamma * bk.rhox[eId_np_j * nc + k] * wghte;
2029  }
2030  // WARNING !!!
2031  if (i < pEnumComB[j])
2032  dFdXsB[(i + 1) * ncolB + jxB + i] += xi * transJ * dP;
2033  }
2034  } else {
2035  // Saturation
2036  for (USI j1 = 0; j1 < np; j1++) {
2037 
2038  wghtb = 0;
2039  wghte = 0;
2040  if (j1 == j && phaseExistB[j]) {
2041  tmp = -dGamma / ((bk.S[bId_np_j] + bk.S[eId_np_j]) *
2042  (bk.S[bId_np_j] + bk.S[eId_np_j]));
2043  wghtb = tmp * bk.rho[bId_np_j] * bk.S[eId_np_j];
2044  wghte = tmp * bk.rho[eId_np_j] * bk.S[bId_np_j];
2045  }
2046 
2047  if (phaseExistB[j1]) {
2048  dFdXsB[(i + 1) * ncolB + j1SB] +=
2049  transIJ * (bk.dPcj_dS[bId_np_j * np + j1] + wghtb);
2050  j1SB++;
2051  }
2052  if (phaseExistE[j1]) {
2053  dFdXsE[(i + 1) * ncolE + j1SE] -=
2054  transIJ * (bk.dPcj_dS[eId_np_j * np + j1] + wghte);
2055  tmp = Akd * xij * xi / mu * bk.dKr_dS[uId_np_j * np + j1] *
2056  dP;
2057  dFdXsE[(i + 1) * ncolE + j1SE] += tmp;
2058  j1SE++;
2059  }
2060  }
2061  // Cij
2062  if (!phaseExistB[j]) {
2063  for (USI k = 0; k < pEnumComE[j]; k++) {
2064  rhox = bk.rhox[uId_np_j * nc + k];
2065  xix = bk.xix[uId_np_j * nc + k];
2066  mux = bk.mux[uId_np_j * nc + k];
2067  tmp = -transIJ * rhox * dGamma;
2068  tmp += xij * transJ * xix * dP;
2069  tmp += -transIJ * mux / mu * dP;
2070  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
2071  }
2072  // WARNING !!!
2073  if (i < pEnumComE[j])
2074  dFdXsE[(i + 1) * ncolE + jxE + i] += xi * transJ * dP;
2075  } else {
2076  wghtb = bk.S[bId_np_j] / (bk.S[bId_np_j] + bk.S[eId_np_j]);
2077  wghte = bk.S[eId_np_j] / (bk.S[bId_np_j] + bk.S[eId_np_j]);
2078  for (USI k = 0; k < pEnumComE[j]; k++) {
2079  rhox = bk.rhox[eId_np_j * nc + k] * wghte;
2080  xix = bk.xix[uId_np_j * nc + k];
2081  mux = bk.mux[uId_np_j * nc + k];
2082  tmp = -transIJ * rhox * dGamma;
2083  tmp += xij * transJ * xix * dP;
2084  tmp += -transIJ * mux / mu * dP;
2085  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
2086  dFdXsB[(i + 1) * ncolB + jxB + k] +=
2087  -transIJ * dGamma * bk.rhox[bId_np_j * nc + k] * wghtb;
2088  }
2089  // WARNING !!!
2090  if (i < pEnumComE[j])
2091  dFdXsE[(i + 1) * ncolE + jxE + i] += xi * transJ * dP;
2092  }
2093  }
2094  }
2095  jxB += pEnumComB[j];
2096  jxE += pEnumComE[j];
2097  }
2098 
2099  // Assemble
2100  bmat = dFdXpB;
2101  DaABpbC(ncol, ncol, ncolB, 1, dFdXsB.data(), &bk.dSec_dPri[bId * lendSdP], 1,
2102  bmat.data());
2103  Dscalar(bsize, dt, bmat.data());
2104  // Begin - Begin -- add
2105  ls.AddDiag(bId, bmat);
2106  // End - Begin -- insert
2107  Dscalar(bsize, -1, bmat.data());
2108  ls.NewOffDiag(eId, bId, bmat);
2109 
2110 #ifdef OCP_NANCHECK
2111  if (!CheckNan(bmat.size(), &bmat[0])) {
2112  OCP_ABORT("INF or NAN in bmat !");
2113  }
2114 #endif
2115 
2116  // End
2117  bmat = dFdXpE;
2118  DaABpbC(ncol, ncol, ncolE, 1, dFdXsE.data(), &bk.dSec_dPri[eId * lendSdP], 1,
2119  bmat.data());
2120  Dscalar(bsize, dt, bmat.data());
2121  // Begin - End -- insert
2122  ls.NewOffDiag(bId, eId, bmat);
2123  // End - End -- add
2124  Dscalar(bsize, -1, bmat.data());
2125  ls.AddDiag(eId, bmat);
2126 
2127 #ifdef OCP_NANCHECK
2128  if (!CheckNan(bmat.size(), &bmat[0])) {
2129  OCP_ABORT("INF or INF in bmat !");
2130  }
2131 #endif
2132  }
2133 }
2134 
2135 void IsoT_FIM::AssembleMatWells(LinearSystem& ls,
2136  const Reservoir& rs,
2137  const OCP_DBL& dt) const
2138 {
2139  for (auto& wl : rs.allWells.wells) {
2140  if (wl.IsOpen()) {
2141 
2142  switch (wl.WellType()) {
2143  case INJ:
2144  AssembleMatInjWells(ls, rs.bulk, wl, dt);
2145  break;
2146  case PROD:
2147  AssembleMatProdWells(ls, rs.bulk, wl, dt);
2148  break;
2149  default:
2150  OCP_ABORT("Wrong well type");
2151  }
2152  }
2153  }
2154 
2156  // for (auto& wG : rs.allWells.wellGroup) {
2157  // if (wG.IfReInj()) {
2158  // for (auto& prod : rs.allWells.wellGroup[wG.prodGroup].wIdPROD) {
2159  // if (rs.allWells.wells[prod].IsOpen()) {
2160  // rs.allWells.wells[prod].AssembleMatReinjection_FIM(rs.bulk, ls,
2161  // dt, rs.allWells.wells,
2162  // wG.wIdINJ);
2163  // }
2164  // }
2165  // }
2166  // }
2167 }
2168 
2169 void IsoT_FIM::AssembleMatInjWells(LinearSystem& ls,
2170  const Bulk& bk,
2171  const Well& wl,
2172  const OCP_DBL& dt) const
2173 {
2174  const USI nc = bk.numCom;
2175  const USI np = bk.numPhase;
2176  const USI ncol = nc + 1;
2177  const USI ncol2 = np * nc + np;
2178  const USI bsize = ncol * ncol;
2179  const USI bsize2 = ncol * ncol2;
2180  OCP_USI n_np_j;
2181 
2182  vector<OCP_DBL> bmat(bsize, 0);
2183  vector<OCP_DBL> bmat2(bsize, 0);
2184  vector<OCP_DBL> dQdXpB(bsize, 0);
2185  vector<OCP_DBL> dQdXpW(bsize, 0);
2186  vector<OCP_DBL> dQdXsB(bsize2, 0);
2187 
2188  OCP_DBL mu, muP, dP, transIJ;
2189 
2190  const OCP_USI wId = ls.AddDim(1) - 1;
2191  ls.NewDiag(wId, bmat);
2192 
2193  for (USI p = 0; p < wl.PerfNum(); p++) {
2194  const OCP_USI n = wl.PerfLocation(p);
2195  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2196  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2197  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2198 
2199  dP = bk.P[n] - wl.BHP() - wl.DG(p);
2200 
2201  for (USI j = 0; j < np; j++) {
2202  n_np_j = n * np + j;
2203  if (!bk.phaseExist[n_np_j]) continue;
2204 
2205  mu = bk.mu[n_np_j];
2206  muP = bk.muP[n_np_j];
2207 
2208  for (USI i = 0; i < nc; i++) {
2209  // dQ / dP
2210  transIJ = wl.PerfTransj(p, j) * wl.PerfXi(p) * wl.InjZi(i);
2211  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
2212  dQdXpW[(i + 1) * ncol] += -transIJ;
2213 
2214  // dQ / dS
2215  for (USI k = 0; k < np; k++) {
2216  dQdXsB[(i + 1) * ncol2 + k] +=
2217  CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * wl.PerfXi(p) *
2218  wl.InjZi(i) * bk.dKr_dS[n_np_j * np + k] * dP / mu;
2219  }
2220  // dQ / dxij
2221  for (USI k = 0; k < nc; k++) {
2222  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] +=
2223  -transIJ * dP / mu * bk.mux[n_np_j * nc + k];
2224  }
2225  }
2226  }
2227 
2228  // Bulk to Well
2229  bmat = dQdXpB;
2230  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2], 1,
2231  bmat.data());
2232  Dscalar(bsize, dt, bmat.data());
2233  // Bulk - Bulk -- add
2234  ls.AddDiag(n, bmat);
2235 
2236  // Bulk - Well -- insert
2237  bmat = dQdXpW;
2238  Dscalar(bsize, dt, bmat.data());
2239  ls.NewOffDiag(n, wId, bmat);
2240 
2241  // Well
2242  switch (wl.OptMode()) {
2243  case RATE_MODE:
2244  case ORATE_MODE:
2245  case GRATE_MODE:
2246  case WRATE_MODE:
2247  case LRATE_MODE:
2248  // Well - Well -- add
2249  fill(bmat.begin(), bmat.end(), 0.0);
2250  for (USI i = 0; i < nc; i++) {
2251  bmat[0] += dQdXpW[(i + 1) * ncol];
2252  bmat[(i + 1) * ncol + i + 1] = 1;
2253  }
2254  ls.AddDiag(wId, bmat);
2255 
2256  // Well - Bulk -- insert
2257  bmat = dQdXpB;
2258  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2],
2259  1, bmat.data());
2260  fill(bmat2.begin(), bmat2.end(), 0.0);
2261  for (USI i = 0; i < nc; i++) {
2262  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
2263  }
2264  ls.NewOffDiag(wId, n, bmat2);
2265  break;
2266 
2267  case BHP_MODE:
2268  // Well - Well -- add
2269  fill(bmat.begin(), bmat.end(), 0.0);
2270  for (USI i = 0; i < ncol; i++) {
2271  bmat[i * ncol + i] = 1;
2272  }
2273  ls.AddDiag(wId, bmat);
2274 
2275  // Well - Bulk -- insert
2276  fill(bmat.begin(), bmat.end(), 0.0);
2277  ls.NewOffDiag(wId, n, bmat);
2278  break;
2279 
2280  default:
2281  OCP_ABORT("Wrong Well Opt mode!");
2282  break;
2283  }
2284  }
2285 }
2286 
2287 void IsoT_FIM::AssembleMatProdWells(LinearSystem& ls,
2288  const Bulk& bk,
2289  const Well& wl,
2290  const OCP_DBL& dt) const
2291 {
2292  const USI nc = bk.numCom;
2293  const USI np = bk.numPhase;
2294  const USI ncol = nc + 1;
2295  const USI ncol2 = np * nc + np;
2296  const USI bsize = ncol * ncol;
2297  const USI bsize2 = ncol * ncol2;
2298  OCP_USI n_np_j;
2299 
2300  vector<OCP_DBL> bmat(bsize, 0);
2301  vector<OCP_DBL> bmat2(bsize, 0);
2302  vector<OCP_DBL> dQdXpB(bsize, 0);
2303  vector<OCP_DBL> dQdXpW(bsize, 0);
2304  vector<OCP_DBL> dQdXsB(bsize2, 0);
2305 
2306  OCP_DBL xij, xi, mu, muP, xiP, dP, transIJ, tmp;
2307 
2308  const OCP_USI wId = ls.AddDim(1) - 1;
2309  ls.NewDiag(wId, bmat);
2310 
2311  // Set Prod Weight
2312  if (wl.OptMode() != BHP_MODE) wl.CalProdWeight(bk);
2313 
2314  for (USI p = 0; p < wl.PerfNum(); p++) {
2315  const OCP_USI n = wl.PerfLocation(p);
2316  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2317  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2318  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2319 
2320  for (USI j = 0; j < np; j++) {
2321  n_np_j = n * np + j;
2322  if (!bk.phaseExist[n_np_j]) continue;
2323 
2324  dP = bk.Pj[n_np_j] - wl.BHP() - wl.DG(p);
2325  xi = bk.xi[n_np_j];
2326  mu = bk.mu[n_np_j];
2327  muP = bk.muP[n_np_j];
2328  xiP = bk.xiP[n_np_j];
2329 
2330  for (USI i = 0; i < nc; i++) {
2331  xij = bk.xij[n_np_j * nc + i];
2332  // dQ / dP
2333  transIJ = wl.PerfTransj(p, j) * xi * xij;
2334  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu) +
2335  dP * wl.PerfTransj(p, j) * xij * xiP;
2336  dQdXpW[(i + 1) * ncol] += -transIJ;
2337 
2338  // dQ / dS
2339  for (USI k = 0; k < np; k++) {
2340  tmp = CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * dP / mu * xi *
2341  xij * bk.dKr_dS[n_np_j * np + k];
2342  // capillary pressure
2343  tmp += transIJ * bk.dPcj_dS[n_np_j * np + k];
2344  dQdXsB[(i + 1) * ncol2 + k] += tmp;
2345  }
2346  // dQ / dCij
2347  for (USI k = 0; k < nc; k++) {
2348  tmp = dP * wl.PerfTransj(p, j) * xij *
2349  (bk.xix[n_np_j * nc + k] - xi / mu * bk.mux[n_np_j * nc + k]);
2350  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] += tmp;
2351  }
2352  dQdXsB[(i + 1) * ncol2 + np + j * nc + i] +=
2353  wl.PerfTransj(p, j) * xi * dP;
2354  }
2355  }
2356 
2357  // Bulk - Bulk -- add
2358  bmat = dQdXpB;
2359  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2], 1,
2360  bmat.data());
2361 
2362  Dscalar(bsize, dt, bmat.data());
2363  ls.AddDiag(n, bmat);
2364 
2365  // Bulk - Well -- insert
2366  bmat = dQdXpW;
2367  Dscalar(bsize, dt, bmat.data());
2368  ls.NewOffDiag(n, wId, bmat);
2369 
2370  // Well
2371  switch (wl.OptMode()) {
2372  case RATE_MODE:
2373  case ORATE_MODE:
2374  case GRATE_MODE:
2375  case WRATE_MODE:
2376  case LRATE_MODE:
2377  // Well - Well -- add
2378  fill(bmat.begin(), bmat.end(), 0.0);
2379  for (USI i = 0; i < nc; i++) {
2380  bmat[0] += dQdXpW[(i + 1) * ncol] * wl.ProdWeight(i);
2381  bmat[(i + 1) * ncol + i + 1] = 1;
2382  }
2383  ls.AddDiag(wId, bmat);
2384 
2385  // Well - Bulk -- insert
2386  bmat = dQdXpB;
2387  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2],
2388  1, bmat.data());
2389  fill(bmat2.begin(), bmat2.end(), 0.0);
2390  for (USI i = 0; i < nc; i++) {
2391  Daxpy(ncol, wl.ProdWeight(i), bmat.data() + (i + 1) * ncol,
2392  bmat2.data());
2393  }
2394  ls.NewOffDiag(wId, n, bmat2);
2395  break;
2396 
2397  case BHP_MODE:
2398  // Well - Well -- add
2399  fill(bmat.begin(), bmat.end(), 0.0);
2400  for (USI i = 0; i < ncol; i++) {
2401  bmat[i * ncol + i] = 1;
2402  }
2403  ls.AddDiag(wId, bmat);
2404 
2405  // Well - Bulk -- insert
2406  fill(bmat.begin(), bmat.end(), 0.0);
2407  ls.NewOffDiag(wId, n, bmat);
2408  break;
2409 
2410  default:
2411  OCP_ABORT("Wrong Well Opt mode!");
2412  break;
2413  }
2414  }
2415 }
2416 
2417 void IsoT_FIM::AssembleMatWellsNew(LinearSystem& ls,
2418  const Reservoir& rs,
2419  const OCP_DBL& dt) const
2420 {
2421  for (auto& wl : rs.allWells.wells) {
2422  if (wl.IsOpen()) {
2423 
2424  switch (wl.WellType()) {
2425  case INJ:
2426  AssembleMatInjWellsNew(ls, rs.bulk, wl, dt);
2427  break;
2428  case PROD:
2429  AssembleMatProdWellsNew(ls, rs.bulk, wl, dt);
2430  break;
2431  default:
2432  OCP_ABORT("Wrong well type");
2433  }
2434  }
2435  }
2436 }
2437 
2438 void IsoT_FIM::AssembleMatInjWellsNew(LinearSystem& ls,
2439  const Bulk& bk,
2440  const Well& wl,
2441  const OCP_DBL& dt) const
2442 {
2443  const USI nc = bk.numCom;
2444  const USI np = bk.numPhase;
2445  const USI lendSdP = bk.maxLendSdP;
2446  const USI ncol = nc + 1;
2447  const USI ncol2 = np * nc + np;
2448  const USI bsize = ncol * ncol;
2449  const USI bsize2 = ncol * ncol2;
2450 
2451  vector<OCP_DBL> bmat(bsize, 0);
2452  vector<OCP_DBL> bmat2(bsize, 0);
2453  vector<OCP_DBL> dQdXpB(bsize, 0);
2454  vector<OCP_DBL> dQdXpW(bsize, 0);
2455  vector<OCP_DBL> dQdXsB(bsize2, 0);
2456  vector<OCP_BOOL> phaseExistB(np, OCP_FALSE);
2457  vector<OCP_BOOL> phasedS_B(np, OCP_FALSE);
2458  vector<USI> pVnumComB(np, 0);
2459  USI ncolB;
2460  OCP_USI n_np_j;
2461 
2462  OCP_DBL mu, muP, dP, transIJ;
2463 
2464  const OCP_USI wId = ls.AddDim(1) - 1;
2465  ls.NewDiag(wId, bmat);
2466 
2467  for (USI p = 0; p < wl.PerfNum(); p++) {
2468  const OCP_USI n = wl.PerfLocation(p);
2469  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2470  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2471  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2472 
2473  dP = bk.P[n] - wl.BHP() - wl.DG(p);
2474 
2475  USI jxB = 0;
2476  ncolB = 0;
2477 
2478  for (USI j = 0; j < np; j++) {
2479  phaseExistB[j] = bk.phaseExist[n * np + j];
2480  phasedS_B[j] = bk.pSderExist[n * np + j];
2481  if (phasedS_B[j]) jxB++;
2482  pVnumComB[j] = bk.pVnumCom[n * np + j];
2483  ncolB += pVnumComB[j];
2484  }
2485  ncolB += jxB;
2486 
2487  for (USI j = 0; j < np; j++) {
2488 
2489  if (!phaseExistB[j]) {
2490  jxB += pVnumComB[j];
2491  continue;
2492  }
2493 
2494  n_np_j = n * np + j;
2495  mu = bk.mu[n_np_j];
2496  muP = bk.muP[n_np_j];
2497 
2498  for (USI i = 0; i < nc; i++) {
2499  // dQ / dP
2500  transIJ = wl.PerfTransj(p, j) * wl.PerfXi(p) * wl.InjZi(i);
2501  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
2502  dQdXpW[(i + 1) * ncol] += -transIJ;
2503 
2504  // dQ / dS
2505  USI j1B = 0;
2506  for (USI j1 = 0; j1 < np; j1++) {
2507  if (phasedS_B[j1]) {
2508  dQdXsB[(i + 1) * ncolB + j1B] +=
2509  CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * wl.PerfXi(p) *
2510  wl.InjZi(i) * bk.dKr_dS[n_np_j * np + j1] * dP / mu;
2511  j1B++;
2512  }
2513  }
2514 
2515  // dQ / dxij
2516  for (USI k = 0; k < pVnumComB[j]; k++) {
2517  dQdXsB[(i + 1) * ncolB + jxB + k] +=
2518  -transIJ * dP / mu * bk.mux[n_np_j * nc + k];
2519  }
2520  }
2521  jxB += pVnumComB[j];
2522  }
2523 
2524  bmat = dQdXpB;
2525  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP], 1,
2526  bmat.data());
2527  Dscalar(bsize, dt, bmat.data());
2528  // Bulk - Bulk -- add
2529  ls.AddDiag(n, bmat);
2530 
2531  // Bulk - Well -- insert
2532  bmat = dQdXpW;
2533  Dscalar(bsize, dt, bmat.data());
2534  ls.NewOffDiag(n, wId, bmat);
2535 
2536  // Well
2537  switch (wl.OptMode()) {
2538  case RATE_MODE:
2539  case ORATE_MODE:
2540  case GRATE_MODE:
2541  case WRATE_MODE:
2542  case LRATE_MODE:
2543  // Well - Well -- add
2544  fill(bmat.begin(), bmat.end(), 0.0);
2545  for (USI i = 0; i < nc; i++) {
2546  bmat[0] += dQdXpW[(i + 1) * ncol];
2547  bmat[(i + 1) * ncol + i + 1] = 1;
2548  }
2549  ls.AddDiag(wId, bmat);
2550 
2551  // Well - Bulk -- insert
2552  bmat = dQdXpB;
2553  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP],
2554  1, bmat.data());
2555  fill(bmat2.begin(), bmat2.end(), 0.0);
2556  for (USI i = 0; i < nc; i++) {
2557  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
2558  }
2559  ls.NewOffDiag(wId, n, bmat2);
2560  break;
2561 
2562  case BHP_MODE:
2563  // Well - Well -- add
2564  fill(bmat.begin(), bmat.end(), 0.0);
2565  for (USI i = 0; i < ncol; i++) {
2566  bmat[i * ncol + i] = 1;
2567  }
2568  ls.AddDiag(wId, bmat);
2569 
2570  // Well - Bulk -- insert
2571  fill(bmat.begin(), bmat.end(), 0.0);
2572  ls.NewOffDiag(wId, n, bmat);
2573  break;
2574 
2575  default:
2576  OCP_ABORT("Wrong Well Opt mode!");
2577  break;
2578  }
2579  }
2580 }
2581 
2582 void IsoT_FIM::AssembleMatProdWellsNew(LinearSystem& ls,
2583  const Bulk& bk,
2584  const Well& wl,
2585  const OCP_DBL& dt) const
2586 {
2587  const USI nc = bk.numCom;
2588  const USI np = bk.numPhase;
2589  const USI lendSdP = bk.maxLendSdP;
2590  const USI ncol = nc + 1;
2591  const USI ncol2 = np * nc + np;
2592  const USI bsize = ncol * ncol;
2593  const USI bsize2 = ncol * ncol2;
2594 
2595  vector<OCP_DBL> bmat(bsize, 0);
2596  vector<OCP_DBL> bmat2(bsize, 0);
2597  vector<OCP_DBL> dQdXpB(bsize, 0);
2598  vector<OCP_DBL> dQdXpW(bsize, 0);
2599  vector<OCP_DBL> dQdXsB(bsize2, 0);
2600  vector<OCP_BOOL> phaseExistB(np, OCP_FALSE);
2601  vector<OCP_BOOL> phasedS_B(np, OCP_FALSE);
2602  vector<USI> pVnumComB(np, 0);
2603  USI ncolB;
2604  OCP_USI n_np_j;
2605 
2606  OCP_DBL xij, xi, mu, muP, xiP, dP, transIJ, tmp;
2607 
2608  const OCP_USI wId = ls.AddDim(1) - 1;
2609  ls.NewDiag(wId, bmat);
2610 
2611  // Set Prod Weight
2612  if (wl.OptMode() != BHP_MODE) wl.CalProdWeight(bk);
2613 
2614  for (USI p = 0; p < wl.PerfNum(); p++) {
2615  OCP_USI n = wl.PerfLocation(p);
2616  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
2617  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
2618  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
2619 
2620  USI jxB = 0;
2621  ncolB = 0;
2622  for (USI j = 0; j < np; j++) {
2623  phaseExistB[j] = bk.phaseExist[n * np + j];
2624  phasedS_B[j] = bk.pSderExist[n * np + j];
2625  if (phasedS_B[j]) jxB++;
2626  pVnumComB[j] = bk.pVnumCom[n * np + j];
2627  ncolB += pVnumComB[j];
2628  }
2629  ncolB += jxB;
2630 
2631  for (USI j = 0; j < np; j++) {
2632 
2633  if (!phaseExistB[j]) {
2634  jxB += pVnumComB[j];
2635  continue;
2636  }
2637 
2638  n_np_j = n * np + j;
2639  dP = bk.Pj[n_np_j] - wl.BHP() - wl.DG(p);
2640  xi = bk.xi[n_np_j];
2641  mu = bk.mu[n_np_j];
2642  muP = bk.muP[n_np_j];
2643  xiP = bk.xiP[n_np_j];
2644 
2645  for (USI i = 0; i < nc; i++) {
2646  xij = bk.xij[n_np_j * nc + i];
2647  // dQ / dP
2648  transIJ = wl.PerfTransj(p, j) * xi * xij;
2649  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu) +
2650  dP * wl.PerfTransj(p, j) * xij * xiP;
2651  dQdXpW[(i + 1) * ncol] += -transIJ;
2652 
2653  // dQ / dS
2654  USI j1B = 0;
2655  for (USI j1 = 0; j1 < np; j1++) {
2656  if (phasedS_B[j1]) {
2657  tmp = CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * dP / mu *
2658  xi * xij * bk.dKr_dS[n_np_j * np + j1];
2659  // capillary pressure
2660  tmp += transIJ * bk.dPcj_dS[n_np_j * np + j1];
2661  dQdXsB[(i + 1) * ncolB + j1B] += tmp;
2662  j1B++;
2663  }
2664  }
2665 
2666  for (USI k = 0; k < pVnumComB[j]; k++) {
2667  tmp = dP * wl.PerfTransj(p, j) * xij *
2668  (bk.xix[n_np_j * nc + k] - xi / mu * bk.mux[n_np_j * nc + k]);
2669  dQdXsB[(i + 1) * ncolB + jxB + k] += tmp;
2670  }
2671  // WARNING !!!
2672  if (i < pVnumComB[j])
2673  dQdXsB[(i + 1) * ncolB + jxB + i] += wl.PerfTransj(p, j) * xi * dP;
2674  }
2675  jxB += pVnumComB[j];
2676  }
2677  // Bulk - Bulk -- add
2678  bmat = dQdXpB;
2679  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP], 1,
2680  bmat.data());
2681  Dscalar(bsize, dt, bmat.data());
2682  ls.AddDiag(n, bmat);
2683  // Bulk - Well -- insert
2684  bmat = dQdXpW;
2685  Dscalar(bsize, dt, bmat.data());
2686  ls.NewOffDiag(n, wId, bmat);
2687 
2688  // Well
2689  switch (wl.OptMode()) {
2690  case RATE_MODE:
2691  case ORATE_MODE:
2692  case GRATE_MODE:
2693  case WRATE_MODE:
2694  case LRATE_MODE:
2695  // Well - Well -- add
2696  fill(bmat.begin(), bmat.end(), 0.0);
2697  for (USI i = 0; i < nc; i++) {
2698  bmat[0] += dQdXpW[(i + 1) * ncol] * wl.ProdWeight(i);
2699  bmat[(i + 1) * ncol + i + 1] = 1;
2700  }
2701  ls.AddDiag(wId, bmat);
2702 
2703  // Well - Bulk -- insert
2704  bmat = dQdXpB;
2705  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP],
2706  1, bmat.data());
2707  fill(bmat2.begin(), bmat2.end(), 0.0);
2708  for (USI i = 0; i < nc; i++) {
2709  Daxpy(ncol, wl.ProdWeight(i), bmat.data() + (i + 1) * ncol,
2710  bmat2.data());
2711  }
2712  ls.NewOffDiag(wId, n, bmat2);
2713  break;
2714 
2715  case BHP_MODE:
2716  // Well - Well -- add
2717  fill(bmat.begin(), bmat.end(), 0.0);
2718  for (USI i = 0; i < ncol; i++) {
2719  bmat[i * ncol + i] = 1;
2720  }
2721  ls.AddDiag(wId, bmat);
2722 
2723  // Well - Bulk -- insert
2724  fill(bmat.begin(), bmat.end(), 0.0);
2725  ls.NewOffDiag(wId, n, bmat);
2726  break;
2727 
2728  default:
2729  OCP_ABORT("Wrong Well Opt mode!");
2730  break;
2731  }
2732  }
2733 }
2734 
2735 void IsoT_FIM::GetSolution(Reservoir& rs,
2736  const vector<OCP_DBL>& u,
2737  const OCPControl& ctrl) const
2738 {
2739  // Bulk
2740  const OCP_DBL dSmaxlim = ctrl.ctrlNR.NRdSmax;
2741  // const OCP_DBL dPmaxlim = ctrl.ctrlNR.NRdPmax;
2742 
2743  Bulk& bk = rs.bulk;
2744  const OCP_USI nb = bk.numBulk;
2745  const USI np = bk.numPhase;
2746  const USI nc = bk.numCom;
2747  const USI row = np * (nc + 1);
2748  const USI col = nc + 1;
2749  vector<OCP_DBL> dtmp(row, 0);
2750  OCP_DBL chopmin = 1;
2751  OCP_DBL choptmp = 0;
2752 
2753  bk.dSNR = bk.S;
2754  bk.NRphaseNum = bk.phaseNum;
2755  bk.NRdPmax = 0;
2756  bk.NRdNmax = 0;
2757 
2758  for (OCP_USI n = 0; n < nb; n++) {
2759  // const vector<OCP_DBL>& scm = satcm[SATNUM[n]];
2760 
2761  chopmin = 1;
2762  // compute the chop
2763  fill(dtmp.begin(), dtmp.end(), 0.0);
2764 
2765 #ifdef OCP_OLD_FIM
2766  DaAxpby(row, col, 1, &bk.dSec_dPri[n * bk.maxLendSdP], u.data() + n * col, 1,
2767  dtmp.data());
2768  const OCP_BOOL newFIM = OCP_FALSE;
2769 #else
2770  DaAxpby(bk.bRowSizedSdP[n], col, 1, &bk.dSec_dPri[n * bk.maxLendSdP],
2771  u.data() + n * col, 1, dtmp.data());
2772  const OCP_BOOL newFIM = OCP_TRUE;
2773 #endif // OCP_OLD_FIM
2774 
2775  USI js = 0;
2776  for (USI j = 0; j < np; j++) {
2777  if (!bk.pSderExist[n * np + j] && newFIM) {
2778  continue;
2779  }
2780 
2781  choptmp = 1;
2782  if (fabs(dtmp[js]) > dSmaxlim) {
2783  choptmp = dSmaxlim / fabs(dtmp[js]);
2784  } else if (bk.S[n * np + j] + dtmp[js] < 0.0) {
2785  choptmp = 0.9 * bk.S[n * np + j] / fabs(dtmp[js]);
2786  }
2787 
2788  // if (fabs(S[n_np_j] - scm[j]) > TINY &&
2789  // (S[n_np_j] - scm[j]) / (choptmp * dtmp[js]) < 0)
2790  // choptmp *= min(1.0, -((S[n_np_j] - scm[j]) / (choptmp * dtmp[js])));
2791 
2792  chopmin = min(chopmin, choptmp);
2793  js++;
2794  }
2795 
2796  // dS
2797  js = 0;
2798  for (USI j = 0; j < np; j++) {
2799  if (!bk.pSderExist[n * np + j] && newFIM) {
2800  bk.dSNRP[n * np + j] = 0;
2801  continue;
2802  }
2803  bk.dSNRP[n * np + j] = chopmin * dtmp[js];
2804  bk.S[n * np + j] += bk.dSNRP[n * np + j];
2805  js++;
2806  }
2807 
2808  // dxij ---- Compositional model only
2809  if (bk.IfUseEoS()) {
2810  if (bk.phaseNum[n] >= 3) {
2811  // num of Hydroncarbon phase >= 2
2812  OCP_USI bId = 0;
2813  for (USI j = 0; j < 2; j++) {
2814  bId = n * np * nc + j * nc;
2815  for (USI i = 0; i < bk.numComH; i++) {
2816  bk.xij[bId + i] += chopmin * dtmp[js];
2817  js++;
2818  }
2819 #ifdef OCP_OLD_FIM
2820  js++;
2821 #endif // OCP_OLD_FIM
2822  }
2823  }
2824  }
2825 
2826  // dP
2827  OCP_DBL dP = u[n * col];
2828  // choptmp = dPmaxlim / fabs(dP);
2829  // chopmin = min(chopmin, choptmp);
2830  if (fabs(bk.NRdPmax) < fabs(dP)) bk.NRdPmax = dP;
2831  bk.P[n] += dP; // seems better
2832  bk.dPNR[n] = dP;
2833 
2834  // dNi
2835  bk.NRstep[n] = chopmin;
2836  for (USI i = 0; i < nc; i++) {
2837  bk.dNNR[n * nc + i] = u[n * col + 1 + i] * chopmin;
2838  if (fabs(bk.NRdNmax) < fabs(bk.dNNR[n * nc + i]) / bk.Nt[n])
2839  bk.NRdNmax = bk.dNNR[n * nc + i] / bk.Nt[n];
2840 
2841  bk.Ni[n * nc + i] += bk.dNNR[n * nc + i];
2842 
2843  // if (bk.Ni[n * nc + i] < 0 && bk.Ni[n * nc + i] > -1E-3) {
2844  // bk.Ni[n * nc + i] = 1E-20;
2845  // }
2846  }
2847  }
2848 
2849  // Well
2850  OCP_USI wId = nb * col;
2851  for (auto& wl : rs.allWells.wells) {
2852  if (wl.IsOpen()) {
2853  wl.SetBHP(wl.BHP() + u[wId]);
2854  wId += col;
2855  }
2856  }
2857 }
2858 
2859 void IsoT_FIM::ResetToLastTimeStep(Reservoir& rs, OCPControl& ctrl)
2860 {
2861  Bulk& bk = rs.bulk;
2862 
2863  // Rock
2864  bk.poro = bk.lporo;
2865  bk.poroP = bk.lporoP;
2866  bk.rockVp = bk.lrockVp;
2867  // Fluid
2868  bk.phaseNum = bk.lphaseNum;
2869  bk.Nt = bk.lNt;
2870  bk.Ni = bk.lNi;
2871  bk.vf = bk.lvf;
2872  bk.P = bk.lP;
2873  bk.Pj = bk.lPj;
2874  bk.Pc = bk.lPc;
2875  bk.phaseExist = bk.lphaseExist;
2876  bk.S = bk.lS;
2877  bk.xij = bk.lxij;
2878  bk.rho = bk.lrho;
2879  bk.xi = bk.lxi;
2880  bk.mu = bk.lmu;
2881  bk.kr = bk.lkr;
2882  // derivatives
2883  bk.vfP = bk.lvfP;
2884  bk.vfi = bk.lvfi;
2885  bk.rhoP = bk.lrhoP;
2886  bk.rhox = bk.lrhox;
2887  bk.xiP = bk.lxiP;
2888  bk.xix = bk.lxix;
2889  bk.muP = bk.lmuP;
2890  bk.mux = bk.lmux;
2891  bk.dPcj_dS = bk.ldPcj_dS;
2892  bk.dKr_dS = bk.ldKr_dS;
2893  // FIM-Specified
2894  bk.bRowSizedSdP = bk.lbRowSizedSdP;
2895  bk.dSec_dPri = bk.ldSec_dPri;
2896  bk.pSderExist = bk.lpSderExist;
2897  bk.pVnumCom = bk.lpVnumCom;
2898 
2899  // Wells
2900  rs.allWells.ResetBHP();
2901  rs.allWells.CalTrans(bk);
2902  rs.allWells.CaldG(bk);
2903  rs.allWells.CalFlux(bk);
2904 
2905  // Optional Features
2906  rs.optFeatures.ResetToLastTimeStep();
2907 
2908  // Iters
2909  ctrl.ResetIterNRLS();
2910 
2911  // Residual
2912  CalRes(rs, ctrl.GetCurDt(), OCP_TRUE);
2913 }
2914 
2915 void IsoT_FIM::UpdateLastTimeStep(Reservoir& rs) const
2916 {
2917  Bulk& bk = rs.bulk;
2918 
2919  // Rock
2920  bk.lporo = bk.poro;
2921  bk.lporoP = bk.poroP;
2922  bk.lrockVp = bk.rockVp;
2923 
2924  // Fluid
2925  bk.lphaseNum = bk.phaseNum;
2926  bk.lNt = bk.Nt;
2927  bk.lNi = bk.Ni;
2928  bk.lvf = bk.vf;
2929  bk.lP = bk.P;
2930  bk.lPj = bk.Pj;
2931  bk.lPc = bk.Pc;
2932  bk.lphaseExist = bk.phaseExist;
2933  bk.lS = bk.S;
2934  bk.lxij = bk.xij;
2935  bk.lrho = bk.rho;
2936  bk.lxi = bk.xi;
2937  bk.lmu = bk.mu;
2938  bk.lkr = bk.kr;
2939 
2940  // derivatives
2941  bk.lvfP = bk.vfP;
2942  bk.lvfi = bk.vfi;
2943  bk.lrhoP = bk.rhoP;
2944  bk.lrhox = bk.rhox;
2945  bk.lxiP = bk.xiP;
2946  bk.lxix = bk.xix;
2947  bk.lmuP = bk.muP;
2948  bk.lmux = bk.mux;
2949  bk.ldPcj_dS = bk.dPcj_dS;
2950  bk.ldKr_dS = bk.dKr_dS;
2951 
2952  // FIM-Specified
2953  bk.lbRowSizedSdP = bk.bRowSizedSdP;
2954  bk.ldSec_dPri = bk.dSec_dPri;
2955  bk.lpSderExist = bk.pSderExist;
2956  bk.lpVnumCom = bk.pVnumCom;
2957 
2958  rs.allWells.UpdateLastTimeStepBHP();
2959  rs.optFeatures.UpdateLastTimeStep();
2960 }
2961 
2963 // OCP_FIMn
2965 
2967 void IsoT_FIMn::Setup(Reservoir& rs, LinearSystem& ls, const OCPControl& ctrl)
2968 {
2969  // Allocate Bulk and BulkConn Memory
2970  AllocateReservoir(rs);
2971  // Allocate memory for internal matrix structure
2972  IsoT_FIM::AllocateLinearSystem(ls, rs, ctrl);
2973 }
2974 
2975 void IsoT_FIMn::InitReservoir(Reservoir& rs) const
2976 {
2977  rs.bulk.InitPTSw(50);
2978 
2979  InitRock(rs.bulk);
2980  CalRock(rs.bulk);
2981 
2982  InitFlash(rs.bulk);
2983  IsoT_FIM::CalKrPc(rs.bulk);
2984 
2985  rs.allWells.InitBHP(rs.bulk);
2986 
2987  UpdateLastTimeStep(rs);
2988 }
2989 
2991 void IsoT_FIMn::Prepare(Reservoir& rs, const OCP_DBL& dt)
2992 {
2993  rs.allWells.PrepareWell(rs.bulk);
2994  IsoT_FIM::CalRes(rs, dt, OCP_TRUE);
2995 }
2996 
2998 void IsoT_FIMn::AssembleMat(LinearSystem& ls,
2999  const Reservoir& rs,
3000  const OCP_DBL& dt) const
3001 {
3002  AssembleMatBulksNew(ls, rs, dt);
3003  AssembleMatWellsNew(ls, rs, dt);
3004  ls.AssembleRhsAccumulate(rs.bulk.res.resAbs);
3005 }
3006 
3008 void IsoT_FIMn::SolveLinearSystem(LinearSystem& ls,
3009  Reservoir& rs,
3010  OCPControl& ctrl) const
3011 {
3012 #ifdef DEBUG
3013  ls.CheckEquation();
3014 #endif // DEBUG
3015 
3016  ls.AssembleMatLinearSolver();
3017 
3018  GetWallTime Timer;
3019  Timer.Start();
3020  int status = ls.Solve();
3021  if (status < 0) {
3022  status = ls.GetNumIters();
3023  }
3024  // cout << "LS step = " << status << endl;
3025 
3026 #ifdef DEBUG
3027  // ls.OutputLinearSystem("testA_FIMn.out", "testb_FIMn.out");
3028  // ls.OutputSolution("testx_FIMn.out");
3029  // ls.CheckSolution();
3030 #endif // DEBUG
3031 
3032  ctrl.RecordTimeLS(Timer.Stop() / 1000);
3033  ctrl.UpdateIterLS(status);
3034  ctrl.UpdateIterNR();
3035 
3036  GetSolution(rs, ls.GetSolution(), ctrl);
3037  // rs.PrintSolFIM(ctrl.workDir + "testPNi.out");
3038  ls.ClearData();
3039 }
3040 
3042 OCP_BOOL IsoT_FIMn::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
3043 {
3044  OCP_DBL& dt = ctrl.current_dt;
3045 
3046  // First check: Ni check and bulk Pressure check
3047  if (!ctrl.Check(rs, {"BulkNi", "BulkP"})) {
3048  ResetToLastTimeStep(rs, ctrl);
3049  cout << "Cut time step size and repeat! current dt = " << fixed
3050  << setprecision(3) << dt << " days\n";
3051  return OCP_FALSE;
3052  }
3053 
3054  // Update reservoir properties
3055  CalFlash(rs.bulk);
3056  IsoT_FIM::CalKrPc(rs.bulk);
3057  CalRock(rs.bulk);
3058 
3059  rs.allWells.CalTrans(rs.bulk);
3060  rs.allWells.CalFlux(rs.bulk);
3061 
3062  CalRes(rs, dt, OCP_FALSE);
3063 
3064  return OCP_TRUE;
3065 }
3066 
3068 OCP_BOOL IsoT_FIMn::FinishNR(Reservoir& rs, OCPControl& ctrl)
3069 {
3070  OCP_USI dSn;
3071 
3072  const OCP_DBL NRdSmax = rs.GetNRdSmax(dSn);
3073  const OCP_DBL NRdPmax = rs.GetNRdPmax();
3074  // const OCP_DBL NRdNmax = rs.GetNRdNmax();
3075 
3076  if (((rs.bulk.res.maxRelRes_V <= rs.bulk.res.maxRelRes0_V * ctrl.ctrlNR.NRtol ||
3077  rs.bulk.res.maxRelRes_V <= ctrl.ctrlNR.NRtol ||
3078  rs.bulk.res.maxRelRes_N <= ctrl.ctrlNR.NRtol) &&
3079  rs.bulk.res.maxWellRelRes_mol <= ctrl.ctrlNR.NRtol) ||
3080  (fabs(NRdPmax) <= ctrl.ctrlNR.NRdPmin &&
3081  fabs(NRdSmax) <= ctrl.ctrlNR.NRdSmin)) {
3082 
3083  if (!ctrl.Check(rs, {"WellP"})) {
3084  ResetToLastTimeStep(rs, ctrl);
3085  return OCP_FALSE;
3086  } else {
3087  return OCP_TRUE;
3088  }
3089 
3090  } else if (ctrl.iterNR > ctrl.ctrlNR.maxNRiter) {
3091  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
3092  ResetToLastTimeStep(rs, ctrl);
3093  cout << "### WARNING: NR not fully converged! Cut time step size and repeat! "
3094  "current dt = "
3095  << fixed << setprecision(3) << ctrl.current_dt << " days\n";
3096  return OCP_FALSE;
3097  } else {
3098  return OCP_FALSE;
3099  }
3100 }
3101 
3103 void IsoT_FIMn::FinishStep(Reservoir& rs, OCPControl& ctrl) const
3104 {
3105  rs.CalIPRT(ctrl.GetCurDt());
3106  rs.CalMaxChange();
3107  UpdateLastTimeStep(rs);
3108  ctrl.CalNextTimeStep(rs, {"dP", "dS", "iter"});
3109 }
3110 
3111 void IsoT_FIMn::AllocateReservoir(Reservoir& rs)
3112 {
3113  IsoT_FIM::AllocateReservoir(rs);
3114 
3115  Bulk& bk = rs.bulk;
3116  const OCP_USI nb = bk.numBulk;
3117  const OCP_USI np = bk.numPhase;
3118  const OCP_USI nc = bk.numCom;
3119 
3120  bk.nj.resize(nb * np);
3121  bk.res_n.resize(nb * (np + np * nc));
3122  bk.resPc.resize(nb);
3123 
3124  bk.lnj.resize(nb * np);
3125  bk.lres_n.resize(nb * (np + np * nc));
3126  bk.lresPc.resize(nb);
3127 }
3128 
3129 void IsoT_FIMn::InitFlash(Bulk& bk) const
3130 {
3131  OCP_FUNCNAME;
3132 
3133  const OCP_USI nb = bk.numBulk;
3134  const OCP_USI np = bk.numPhase;
3135  const OCP_USI nc = bk.numCom;
3136 
3137  if (bk.ifComps) {
3138  for (OCP_USI n = 0; n < nb; n++) {
3139  bk.flashCal[bk.PVTNUM[n]]->InitFlashFIMn(bk.P[n], bk.Pb[n], bk.T[n],
3140  &bk.S[n * np], bk.rockVp[n],
3141  bk.Ni.data() + n * nc, n);
3142  for (USI i = 0; i < nc; i++) {
3143  bk.Ni[n * nc + i] = bk.flashCal[bk.PVTNUM[n]]->GetNi(i);
3144  }
3145  PassFlashValue(bk, n);
3146  }
3147  } else {
3148  OCP_ABORT("Not Completed in BLKOIL MODEL!");
3149  }
3150 }
3151 
3152 void IsoT_FIMn::CalFlash(Bulk& bk)
3153 {
3154 
3155  const OCP_USI nb = bk.numBulk;
3156  const OCP_USI np = bk.numPhase;
3157  const OCP_USI nc = bk.numCom;
3158 
3159  if (bk.ifComps) {
3160  vector<USI> flagB(np, 0);
3161  bk.maxNRdSSP = 0;
3162  bk.index_maxNRdSSP = 0;
3163 
3164  for (OCP_USI n = 0; n < nb; n++) {
3165 
3166  for (USI j = 0; j < np; j++) flagB[j] = bk.phaseExist[n * np + j];
3167 
3168  bk.flashCal[bk.PVTNUM[n]]->FlashFIMn(
3169  bk.P[n], bk.T[n], &bk.Ni[n * nc], &bk.S[n * np], &bk.xij[n * np * nc],
3170  &bk.nj[n * np], &flagB[0], bk.phaseNum[n], n);
3171 
3172  PassFlashValue(bk, n);
3173  }
3174  } else {
3175  OCP_ABORT("Not completed!");
3176  }
3177 }
3178 
3179 void IsoT_FIMn::PassFlashValue(Bulk& bk, const OCP_USI& n) const
3180 {
3181  IsoT_FIM::PassFlashValue(bk, n);
3182 
3183  const USI np = bk.numPhase;
3184  const USI nc = bk.numCom;
3185  const USI pvtnum = bk.PVTNUM[n];
3186 
3187  for (USI j = 0; j < bk.numPhase; j++) {
3188  if (bk.phaseExist[n * np + j]) {
3189  bk.nj[n * np + j] = bk.flashCal[pvtnum]->GetNj(j);
3190  }
3191  }
3192 
3193  Dcopy(bk.bRowSizedSdP[n], &bk.res_n[0] + n * (np * (nc + 1)),
3194  &bk.flashCal[pvtnum]->GetRes()[0]);
3195  bk.resPc[n] = bk.flashCal[pvtnum]->GetResPc();
3196 }
3197 
3198 void IsoT_FIMn::AssembleMatBulksNew(LinearSystem& ls,
3199  const Reservoir& rs,
3200  const OCP_DBL& dt) const
3201 {
3202 
3203  const Bulk& bk = rs.bulk;
3204  const BulkConn& conn = rs.conn;
3205  const OCP_USI nb = bk.numBulk;
3206  const USI np = bk.numPhase;
3207  const USI nc = bk.numCom;
3208  const USI ncol = nc + 1;
3209  const USI ncol2 = np * nc + np;
3210  const USI bsize = ncol * ncol;
3211  const USI bsize2 = ncol * ncol2;
3212  const USI lendSdP = bk.maxLendSdP;
3213 
3214  ls.AddDim(nb);
3215 
3216  vector<OCP_DBL> bmat(bsize, 0);
3217  vector<OCP_DBL> tmpb(ncol, 0);
3218 
3219  // Accumulation term
3220  for (USI i = 1; i < ncol; i++) {
3221  bmat[i * ncol + i] = 1;
3222  }
3223  for (OCP_USI n = 0; n < nb; n++) {
3224  bmat[0] = bk.v[n] * bk.poroP[n] - bk.vfP[n];
3225  for (USI i = 0; i < nc; i++) {
3226  bmat[i + 1] = -bk.vfi[n * nc + i];
3227  }
3228 
3229  ls.NewDiag(n, bmat);
3230  ls.AddRhs(n * ncol, -bk.resPc[n]);
3231  }
3232 
3233  // flux term
3234  OCP_DBL Akd;
3235  OCP_DBL transJ, transIJ;
3236  vector<OCP_DBL> dFdXpB(bsize, 0);
3237  vector<OCP_DBL> dFdXpE(bsize, 0);
3238  vector<OCP_DBL> dFdXsB(bsize2, 0);
3239  vector<OCP_DBL> dFdXsE(bsize2, 0);
3240  vector<OCP_BOOL> phaseExistB(np, OCP_FALSE);
3241  vector<OCP_BOOL> phaseExistE(np, OCP_FALSE);
3242  OCP_BOOL phaseExistU;
3243  vector<OCP_BOOL> phasedS_B(np, OCP_FALSE);
3244  vector<OCP_BOOL> phasedS_E(np, OCP_FALSE);
3245  vector<USI> pVnumComB(np, 0);
3246  vector<USI> pVnumComE(np, 0);
3247  USI ncolB, ncolE;
3248 
3249  OCP_USI bId, eId, uId;
3250  OCP_USI bId_np_j, eId_np_j, uId_np_j;
3251  OCP_DBL kr, mu, xi, xij, rhoP, xiP, muP, rhox, xix, mux;
3252  OCP_DBL dP, dGamma;
3253  OCP_DBL tmp;
3254 
3255  for (OCP_USI c = 0; c < conn.numConn; c++) {
3256  bId = conn.iteratorConn[c].BId();
3257  eId = conn.iteratorConn[c].EId();
3258  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
3259  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
3260  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
3261  fill(dFdXsB.begin(), dFdXsB.end(), 0.0);
3262  fill(dFdXsE.begin(), dFdXsE.end(), 0.0);
3263  dGamma = GRAVITY_FACTOR * (bk.depth[bId] - bk.depth[eId]);
3264 
3265  const USI npB = bk.phaseNum[bId];
3266  const USI npE = bk.phaseNum[eId];
3267 
3268  USI jxB = 0;
3269  USI jxE = 0;
3270  ncolB = 0;
3271  ncolE = 0;
3272  for (USI j = 0; j < np; j++) {
3273  phaseExistB[j] = bk.phaseExist[bId * np + j];
3274  phaseExistE[j] = bk.phaseExist[eId * np + j];
3275  phasedS_B[j] = bk.pSderExist[bId * np + j];
3276  phasedS_E[j] = bk.pSderExist[eId * np + j];
3277  if (phasedS_B[j]) jxB++;
3278  if (phasedS_E[j]) jxE++;
3279  pVnumComB[j] = bk.pVnumCom[bId * np + j];
3280  pVnumComE[j] = bk.pVnumCom[eId * np + j];
3281  ncolB += pVnumComB[j];
3282  ncolE += pVnumComE[j];
3283  }
3284  ncolB += jxB;
3285  ncolE += jxE;
3286 
3287  for (USI j = 0; j < np; j++) {
3288  uId = conn.upblock[c * np + j];
3289 
3290  phaseExistU = (uId == bId ? phaseExistB[j] : phaseExistE[j]);
3291  if (!phaseExistU) {
3292  jxB += pVnumComB[j];
3293  jxE += pVnumComE[j];
3294  continue;
3295  }
3296 
3297  bId_np_j = bId * np + j;
3298  eId_np_j = eId * np + j;
3299  uId_np_j = uId * np + j;
3300  dP = bk.Pj[bId_np_j] - bk.Pj[eId_np_j] -
3301  conn.upblock_Rho[c * np + j] * dGamma;
3302  xi = bk.xi[uId_np_j];
3303  kr = bk.kr[uId_np_j];
3304  mu = bk.mu[uId_np_j];
3305  muP = bk.muP[uId_np_j];
3306  xiP = bk.xiP[uId_np_j];
3307  rhoP = bk.rhoP[uId_np_j];
3308  transJ = Akd * kr / mu;
3309 
3310  for (USI i = 0; i < nc; i++) {
3311  xij = bk.xij[uId_np_j * nc + i];
3312  transIJ = xij * xi * transJ;
3313 
3314  // Pressure -- Primary var
3315  dFdXpB[(i + 1) * ncol] += transIJ;
3316  dFdXpE[(i + 1) * ncol] -= transIJ;
3317 
3318  tmp = xij * transJ * xiP * dP;
3319  tmp += -transIJ * muP / mu * dP;
3320  if (!phaseExistE[j]) {
3321  tmp += transIJ * (-rhoP * dGamma);
3322  dFdXpB[(i + 1) * ncol] += tmp;
3323  } else if (!phaseExistB[j]) {
3324  tmp += transIJ * (-rhoP * dGamma);
3325  dFdXpE[(i + 1) * ncol] += tmp;
3326  } else {
3327  dFdXpB[(i + 1) * ncol] +=
3328  transIJ * (-bk.rhoP[bId_np_j] * dGamma) / 2;
3329  dFdXpE[(i + 1) * ncol] +=
3330  transIJ * (-bk.rhoP[eId_np_j] * dGamma) / 2;
3331  if (bId == uId) {
3332  dFdXpB[(i + 1) * ncol] += tmp;
3333  } else {
3334  dFdXpE[(i + 1) * ncol] += tmp;
3335  }
3336  }
3337 
3338  // Second var
3339  USI j1SB = 0;
3340  USI j1SE = 0;
3341  if (bId == uId) {
3342  // Saturation
3343  for (USI j1 = 0; j1 < np; j1++) {
3344  if (phasedS_B[j1]) {
3345  dFdXsB[(i + 1) * ncolB + j1SB] +=
3346  transIJ * bk.dPcj_dS[bId_np_j * np + j1];
3347  tmp = Akd * xij * xi / mu * bk.dKr_dS[uId_np_j * np + j1] *
3348  dP;
3349  dFdXsB[(i + 1) * ncolB + j1SB] += tmp;
3350  j1SB++;
3351  }
3352  if (phasedS_E[j1]) {
3353  dFdXsE[(i + 1) * ncolE + j1SE] -=
3354  transIJ * bk.dPcj_dS[eId_np_j * np + j1];
3355  j1SE++;
3356  }
3357  }
3358  // Cij
3359  if (!phaseExistE[j]) {
3360  for (USI k = 0; k < pVnumComB[j]; k++) {
3361  rhox = bk.rhox[uId_np_j * nc + k];
3362  xix = bk.xix[uId_np_j * nc + k];
3363  mux = bk.mux[uId_np_j * nc + k];
3364  tmp = -transIJ * rhox * dGamma;
3365  tmp += xij * transJ * xix * dP;
3366  tmp += -transIJ * mux / mu * dP;
3367  tmp -= transIJ * dP / bk.nj[uId_np_j];
3368  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
3369  }
3370  // WARNING !!!
3371  if (i < pVnumComB[j])
3372  dFdXsB[(i + 1) * ncolB + jxB + i] +=
3373  xi * transJ * dP / bk.nj[uId_np_j];
3374  } else {
3375  for (USI k = 0; k < pVnumComB[j]; k++) {
3376  rhox = bk.rhox[bId_np_j * nc + k] / 2;
3377  xix = bk.xix[uId_np_j * nc + k];
3378  mux = bk.mux[uId_np_j * nc + k];
3379  tmp = -transIJ * rhox * dGamma;
3380  tmp += xij * transJ * xix * dP;
3381  tmp += -transIJ * mux / mu * dP;
3382  tmp -= transIJ * dP / bk.nj[uId_np_j];
3383  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
3384  dFdXsE[(i + 1) * ncolE + jxE + k] +=
3385  -transIJ * bk.rhox[eId_np_j * nc + k] / 2 * dGamma;
3386  }
3387  // WARNING !!!
3388  if (i < pVnumComB[j])
3389  dFdXsB[(i + 1) * ncolB + jxB + i] +=
3390  xi * transJ * dP / bk.nj[uId_np_j];
3391  }
3392  } else {
3393  // Saturation
3394  for (USI j1 = 0; j1 < np; j1++) {
3395  if (phasedS_B[j1]) {
3396  dFdXsB[(i + 1) * ncolB + j1SB] +=
3397  transIJ * bk.dPcj_dS[bId_np_j * np + j1];
3398  j1SB++;
3399  }
3400  if (phasedS_E[j1]) {
3401  dFdXsE[(i + 1) * ncolE + j1SE] -=
3402  transIJ * bk.dPcj_dS[eId_np_j * np + j1];
3403  tmp = Akd * xij * xi / mu * bk.dKr_dS[uId_np_j * np + j1] *
3404  dP;
3405  dFdXsE[(i + 1) * ncolE + j1SE] += tmp;
3406  j1SE++;
3407  }
3408  }
3409  // Cij
3410  if (!phaseExistB[j]) {
3411  for (USI k = 0; k < pVnumComE[j]; k++) {
3412  rhox = bk.rhox[uId_np_j * nc + k];
3413  xix = bk.xix[uId_np_j * nc + k];
3414  mux = bk.mux[uId_np_j * nc + k];
3415  tmp = -transIJ * rhox * dGamma;
3416  tmp += xij * transJ * xix * dP;
3417  tmp += -transIJ * mux / mu * dP;
3418  tmp -= transIJ * dP / bk.nj[uId_np_j];
3419  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
3420  }
3421  // WARNING !!!
3422  if (i < pVnumComE[j])
3423  dFdXsE[(i + 1) * ncolE + jxE + i] +=
3424  xi * transJ * dP / bk.nj[uId_np_j];
3425  } else {
3426  for (USI k = 0; k < pVnumComE[j]; k++) {
3427  rhox = bk.rhox[eId_np_j * nc + k] / 2;
3428  xix = bk.xix[uId_np_j * nc + k];
3429  mux = bk.mux[uId_np_j * nc + k];
3430  tmp = -transIJ * rhox * dGamma;
3431  tmp += xij * transJ * xix * dP;
3432  tmp += -transIJ * mux / mu * dP;
3433  tmp -= transIJ * dP / bk.nj[uId_np_j];
3434  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
3435  dFdXsB[(i + 1) * ncolB + jxB + k] +=
3436  -transIJ * bk.rhox[bId_np_j * nc + k] / 2 * dGamma;
3437  }
3438  // WARNING !!!
3439  if (i < pVnumComE[j])
3440  dFdXsE[(i + 1) * ncolE + jxE + i] +=
3441  xi * transJ * dP / bk.nj[uId_np_j];
3442  }
3443  }
3444  }
3445  jxB += pVnumComB[j];
3446  jxE += pVnumComE[j];
3447  }
3448 
3449  // Assemble rhs
3450  // Begin
3451  if (npB > 2) {
3452  fill(tmpb.begin(), tmpb.end(), 0.0);
3453  DaAxpby(ncol, ncolB, -1.0, dFdXsB.data(), &bk.res_n[bId * ncol2], 1.0,
3454  &tmpb[0]);
3455  ls.AddRhs(bId, tmpb);
3456  }
3457 
3458  // Assemble mat
3459  bmat = dFdXpB;
3460  DaABpbC(ncol, ncol, ncolB, 1, dFdXsB.data(), &bk.dSec_dPri[bId * lendSdP], 1,
3461  bmat.data());
3462  Dscalar(bsize, dt, bmat.data());
3463  // Begin - Begin -- add
3464  ls.AddDiag(bId, bmat);
3465  // End - Begin -- insert
3466  Dscalar(bsize, -1, bmat.data());
3467  ls.NewOffDiag(eId, bId, bmat);
3468 
3469 #ifdef OCP_NANCHECK
3470  if (!CheckNan(bmat.size(), &bmat[0])) {
3471  OCP_ABORT("INF or NAN in bmat !");
3472  }
3473 #endif
3474 
3475  // Assemble rhs
3476  // End
3477  if (npE > 2) {
3478  fill(tmpb.begin(), tmpb.end(), 0.0);
3479  DaAxpby(ncol, ncolE, -1.0, dFdXsE.data(), &bk.res_n[eId * ncol2], 1.0,
3480  &tmpb[0]);
3481  ls.AddRhs(eId, tmpb);
3482  }
3483 
3484  // End
3485  bmat = dFdXpE;
3486  DaABpbC(ncol, ncol, ncolE, 1, dFdXsE.data(), &bk.dSec_dPri[eId * lendSdP], 1,
3487  bmat.data());
3488  Dscalar(bsize, dt, bmat.data());
3489  // Begin - End -- insert
3490  ls.NewOffDiag(bId, eId, bmat);
3491  // End - End -- add
3492  Dscalar(bsize, -1, bmat.data());
3493  ls.AddDiag(eId, bmat);
3494 
3495 #ifdef OCP_NANCHECK
3496  if (!CheckNan(bmat.size(), &bmat[0])) {
3497  OCP_ABORT("INF or INF in bmat !");
3498  }
3499 #endif
3500  }
3501 }
3502 
3503 void IsoT_FIMn::AssembleMatWellsNew(LinearSystem& ls,
3504  const Reservoir& rs,
3505  const OCP_DBL& dt) const
3506 {
3507  for (auto& wl : rs.allWells.wells) {
3508  if (wl.IsOpen()) {
3509 
3510  switch (wl.WellType()) {
3511  case INJ:
3512  AssembleMatInjWellsNew(ls, rs.bulk, wl, dt);
3513  break;
3514  case PROD:
3515  AssembleMatProdWellsNew(ls, rs.bulk, wl, dt);
3516  break;
3517  default:
3518  OCP_ABORT("Wrong well type");
3519  }
3520  }
3521  }
3522 }
3523 
3524 void IsoT_FIMn::AssembleMatInjWellsNew(LinearSystem& ls,
3525  const Bulk& bk,
3526  const Well& wl,
3527  const OCP_DBL& dt) const
3528 {
3529 
3530  const USI nc = bk.numCom;
3531  const USI np = bk.numPhase;
3532  const USI lendSdP = bk.maxLendSdP;
3533  const USI ncol = nc + 1;
3534  const USI ncol2 = np * nc + np;
3535  const USI bsize = ncol * ncol;
3536  const USI bsize2 = ncol * ncol2;
3537  OCP_USI n_np_j;
3538 
3539  vector<OCP_DBL> tmpb(ncol, 0);
3540  vector<OCP_DBL> bmat(bsize, 0);
3541  vector<OCP_DBL> bmat2(bsize, 0);
3542  vector<OCP_DBL> dQdXpB(bsize, 0);
3543  vector<OCP_DBL> dQdXpW(bsize, 0);
3544  vector<OCP_DBL> dQdXsB(bsize2, 0);
3545  vector<char> phaseExistB(np, OCP_FALSE);
3546  vector<USI> pVnumComB(np, 0);
3547  vector<char> phasedS_B(np, OCP_FALSE);
3548  USI ncolB;
3549 
3550  OCP_DBL mu, muP, dP, transIJ;
3551 
3552  const OCP_USI wId = ls.AddDim(1) - 1;
3553  ls.NewDiag(wId, bmat);
3554 
3555  for (USI p = 0; p < wl.PerfNum(); p++) {
3556  const OCP_USI n = wl.PerfLocation(p);
3557  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
3558  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
3559  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
3560 
3561  dP = bk.P[n] - wl.BHP() - wl.DG(p);
3562 
3563  const USI npB = bk.phaseNum[n];
3564  USI jxB = 0;
3565  ncolB = 0;
3566  for (USI j = 0; j < np; j++) {
3567  phaseExistB[j] = bk.phaseExist[n * np + j];
3568  phasedS_B[j] = bk.pSderExist[n * np + j];
3569  if (phasedS_B[j]) jxB++;
3570  pVnumComB[j] = bk.pVnumCom[n * np + j];
3571  ncolB += pVnumComB[j];
3572  }
3573  ncolB += jxB;
3574 
3575  for (USI j = 0; j < np; j++) {
3576 
3577  if (!phaseExistB[j]) {
3578  jxB += pVnumComB[j];
3579  continue;
3580  }
3581 
3582  n_np_j = n * np + j;
3583  mu = bk.mu[n_np_j];
3584  muP = bk.muP[n_np_j];
3585 
3586  for (USI i = 0; i < nc; i++) {
3587  // dQ / dP
3588  transIJ = wl.PerfTransj(p, j) * wl.PerfXi(p) * wl.InjZi(i);
3589  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
3590  dQdXpW[(i + 1) * ncol] += -transIJ;
3591 
3592  // dQ / dS
3593  USI j1B = 0;
3594  for (USI j1 = 0; j1 < np; j1++) {
3595  if (phasedS_B[j1]) {
3596  dQdXsB[(i + 1) * ncolB + j1B] +=
3597  CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * wl.PerfXi(p) *
3598  wl.InjZi(i) * bk.dKr_dS[n_np_j * np + j1] * dP / mu;
3599  j1B++;
3600  }
3601  }
3602 
3603  // dQ / dxij
3604  for (USI k = 0; k < pVnumComB[j]; k++) {
3605  dQdXsB[(i + 1) * ncolB + jxB + k] +=
3606  -transIJ * dP / mu * bk.mux[n_np_j * nc + k];
3607  }
3608  }
3609  jxB += pVnumComB[j];
3610  }
3611 
3612  // Assemble rhs
3613  // Well
3614  if (npB > 2) {
3615  fill(tmpb.begin(), tmpb.end(), 0.0);
3616  DaAxpby(ncol, ncolB, -1.0, dQdXsB.data(), &bk.res_n[n * ncol2], 1.0,
3617  &tmpb[0]);
3618  ls.AddRhs(n, tmpb);
3619  }
3620 
3621  // Bulk to Well
3622  bmat = dQdXpB;
3623  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP], 1,
3624  bmat.data());
3625  Dscalar(bsize, dt, bmat.data());
3626  // Bulk - Bulk -- add
3627  ls.AddDiag(n, bmat);
3628  // Bulk - Well -- insert
3629  bmat = dQdXpW;
3630  Dscalar(bsize, dt, bmat.data());
3631  ls.NewOffDiag(n, wId, bmat);
3632 
3633  // Well
3634  switch (wl.OptMode()) {
3635  case RATE_MODE:
3636  case ORATE_MODE:
3637  case GRATE_MODE:
3638  case WRATE_MODE:
3639  case LRATE_MODE:
3640  // Well - Well -- add
3641  fill(bmat.begin(), bmat.end(), 0.0);
3642  for (USI i = 0; i < nc; i++) {
3643  bmat[0] += dQdXpW[(i + 1) * ncol];
3644  bmat[(i + 1) * ncol + i + 1] = 1;
3645  }
3646  ls.AddDiag(wId, bmat);
3647 
3648  // Well - Bulk -- insert
3649  bmat = dQdXpB;
3650  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP],
3651  1, bmat.data());
3652  fill(bmat2.begin(), bmat2.end(), 0.0);
3653  for (USI i = 0; i < nc; i++) {
3654  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
3655  }
3656  ls.NewOffDiag(wId, n, bmat2);
3657  break;
3658 
3659  case BHP_MODE:
3660  // Well - Well -- add
3661  fill(bmat.begin(), bmat.end(), 0.0);
3662  for (USI i = 0; i < ncol; i++) {
3663  bmat[i * ncol + i] = 1;
3664  }
3665  ls.AddDiag(wId, bmat);
3666 
3667  // Well - Bulk -- insert
3668  fill(bmat.begin(), bmat.end(), 0.0);
3669  ls.NewOffDiag(wId, n, bmat);
3670  break;
3671 
3672  default:
3673  OCP_ABORT("Wrong Well Opt mode!");
3674  break;
3675  }
3676  }
3677 }
3678 
3679 void IsoT_FIMn::AssembleMatProdWellsNew(LinearSystem& ls,
3680  const Bulk& bk,
3681  const Well& wl,
3682  const OCP_DBL& dt) const
3683 {
3684  const USI nc = bk.numCom;
3685  const USI np = bk.numPhase;
3686  const USI lendSdP = bk.maxLendSdP;
3687  const USI ncol = nc + 1;
3688  const USI ncol2 = np * nc + np;
3689  const USI bsize = ncol * ncol;
3690  const USI bsize2 = ncol * ncol2;
3691  OCP_USI n_np_j;
3692 
3693  vector<OCP_DBL> tmpb(ncol, 0);
3694  vector<OCP_DBL> bmat(bsize, 0);
3695  vector<OCP_DBL> bmat2(bsize, 0);
3696  vector<OCP_DBL> dQdXpB(bsize, 0);
3697  vector<OCP_DBL> dQdXpW(bsize, 0);
3698  vector<OCP_DBL> dQdXsB(bsize2, 0);
3699  vector<char> phaseExistB(np, OCP_FALSE);
3700  vector<char> phasedS_B(np, OCP_FALSE);
3701  vector<USI> pVnumComB(np, 0);
3702  USI ncolB;
3703 
3704  OCP_DBL xij, xi, mu, muP, xiP, dP, transIJ, tmp;
3705 
3706  const OCP_USI wId = ls.AddDim(1) - 1;
3707  ls.NewDiag(wId, bmat);
3708 
3709  // Set Prod Weight
3710  if (wl.OptMode() != BHP_MODE) wl.CalProdWeight(bk);
3711 
3712  for (USI p = 0; p < wl.PerfNum(); p++) {
3713  OCP_USI n = wl.PerfLocation(p);
3714  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
3715  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
3716  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
3717 
3718  const USI npB = bk.phaseNum[n];
3719  USI jxB = 0;
3720  ncolB = 0;
3721  for (USI j = 0; j < np; j++) {
3722  phaseExistB[j] = bk.phaseExist[n * np + j];
3723  phasedS_B[j] = bk.pSderExist[n * np + j];
3724  if (phasedS_B[j]) jxB++;
3725  pVnumComB[j] = bk.pVnumCom[n * np + j];
3726  ncolB += pVnumComB[j];
3727  }
3728  ncolB += jxB;
3729 
3730  for (USI j = 0; j < np; j++) {
3731 
3732  if (!phaseExistB[j]) {
3733  jxB += pVnumComB[j];
3734  continue;
3735  }
3736 
3737  n_np_j = n * np + j;
3738  dP = bk.Pj[n_np_j] - wl.BHP() - wl.DG(p);
3739  xi = bk.xi[n_np_j];
3740  mu = bk.mu[n_np_j];
3741  muP = bk.muP[n_np_j];
3742  xiP = bk.xiP[n_np_j];
3743 
3744  for (USI i = 0; i < nc; i++) {
3745  xij = bk.xij[n_np_j * nc + i];
3746  // dQ / dP
3747  transIJ = wl.PerfTransj(p, j) * xi * xij;
3748  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu) +
3749  dP * wl.PerfTransj(p, j) * xij * xiP;
3750  dQdXpW[(i + 1) * ncol] += -transIJ;
3751 
3752  // dQ / dS
3753  USI j1B = 0;
3754  for (USI j1 = 0; j1 < np; j1++) {
3755  if (phasedS_B[j1]) {
3756  tmp = CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * dP / mu *
3757  xi * xij * bk.dKr_dS[n_np_j * np + j1];
3758  // capillary pressure
3759  tmp += transIJ * bk.dPcj_dS[n_np_j * np + j1];
3760  dQdXsB[(i + 1) * ncolB + j1B] += tmp;
3761  j1B++;
3762  }
3763  }
3764 
3765  for (USI k = 0; k < pVnumComB[j]; k++) {
3766  tmp = dP * wl.PerfTransj(p, j) * xij *
3767  (bk.xix[n_np_j * nc + k] - xi / mu * bk.mux[n_np_j * nc + k]);
3768  tmp -= transIJ * dP / bk.nj[n_np_j];
3769  dQdXsB[(i + 1) * ncolB + jxB + k] += tmp;
3770  }
3771  // WARNING !!!
3772  if (i < pVnumComB[j])
3773  dQdXsB[(i + 1) * ncolB + jxB + i] +=
3774  wl.PerfTransj(p, j) * xi * dP / bk.nj[n_np_j];
3775  }
3776  jxB += pVnumComB[j];
3777  }
3778 
3779  // Assemble rhs
3780  // Well
3781  if (npB > 2) {
3782  fill(tmpb.begin(), tmpb.end(), 0.0);
3783  DaAxpby(ncol, ncolB, -1.0, dQdXsB.data(), &bk.res_n[n * ncol2], 1.0,
3784  &tmpb[0]);
3785  ls.AddRhs(n, tmpb);
3786  }
3787 
3788  // Bulk - Bulk -- add
3789  bmat = dQdXpB;
3790  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP], 1,
3791  bmat.data());
3792  Dscalar(bsize, dt, bmat.data());
3793  ls.AddDiag(n, bmat);
3794  // Bulk - Well -- insert
3795  bmat = dQdXpW;
3796  Dscalar(bsize, dt, bmat.data());
3797  ls.NewOffDiag(n, wId, bmat);
3798 
3799  // Well
3800  switch (wl.OptMode()) {
3801  case RATE_MODE:
3802  case ORATE_MODE:
3803  case GRATE_MODE:
3804  case WRATE_MODE:
3805  case LRATE_MODE:
3806  // Well - Well -- add
3807  fill(bmat.begin(), bmat.end(), 0.0);
3808  for (USI i = 0; i < nc; i++) {
3809  bmat[0] += dQdXpW[(i + 1) * ncol] * wl.ProdWeight(i);
3810  bmat[(i + 1) * ncol + i + 1] = 1;
3811  }
3812  ls.AddDiag(wId, bmat);
3813 
3814  // Well - Bulk -- insert
3815  bmat = dQdXpB;
3816  DaABpbC(ncol, ncol, ncolB, 1, dQdXsB.data(), &bk.dSec_dPri[n * lendSdP],
3817  1, bmat.data());
3818  fill(bmat2.begin(), bmat2.end(), 0.0);
3819  for (USI i = 0; i < nc; i++) {
3820  Daxpy(ncol, wl.ProdWeight(i), bmat.data() + (i + 1) * ncol,
3821  bmat2.data());
3822  }
3823  ls.NewOffDiag(wId, n, bmat2);
3824  break;
3825 
3826  case BHP_MODE:
3827  // Well - Well -- add
3828  fill(bmat.begin(), bmat.end(), 0.0);
3829  for (USI i = 0; i < ncol; i++) {
3830  bmat[i * ncol + i] = 1;
3831  }
3832  ls.AddDiag(wId, bmat);
3833 
3834  // Well - Bulk -- insert
3835  fill(bmat.begin(), bmat.end(), 0.0);
3836  ls.NewOffDiag(wId, n, bmat);
3837  break;
3838 
3839  default:
3840  OCP_ABORT("Wrong Well Opt mode!");
3841  break;
3842  }
3843  }
3844 }
3845 
3847 void IsoT_FIMn::GetSolution(Reservoir& rs,
3848  const vector<OCP_DBL>& u,
3849  const OCPControl& ctrl) const
3850 {
3851  // For saturations changes:
3852  // 1. maximum changes must be less than dSmaxlim,
3853  // 2. if phases become mobile/immobile, then set it to crtical point,
3854  // 3. After saturations are determined, then scale the nij to conserve Volume
3855  // equations.
3856 
3857  // Bulk
3858  const OCP_DBL dSmaxlim = ctrl.ctrlNR.NRdSmax;
3859  // const OCP_DBL dPmaxlim = ctrl.ctrlNR.NRdPmax;
3860 
3861  Bulk& bk = rs.bulk;
3862  const OCP_USI nb = bk.numBulk;
3863  const USI np = bk.numPhase;
3864  const USI nc = bk.numCom;
3865  const USI row = np * (nc + 1);
3866  const USI col = nc + 1;
3867  vector<OCP_DBL> dtmp(row, 0);
3868  vector<OCP_DBL> tmpNij(np * nc, 0);
3869  OCP_USI n_np_j;
3870 
3871  bk.dSNR = bk.S;
3872  bk.NRphaseNum = bk.phaseNum;
3873  bk.NRdPmax = 0;
3874  bk.NRdNmax = 0;
3875 
3876  OCP_DBL dSmax;
3877  OCP_DBL dP;
3878  OCP_DBL chop = 1;
3879 
3880  for (OCP_USI n = 0; n < nb; n++) {
3881  const vector<OCP_DBL>& scm = bk.satcm[bk.SATNUM[n]];
3882 
3883  dP = u[n * col];
3884  bk.dPNR[n] = dP;
3885  bk.P[n] += dP; // seems better
3886  if (fabs(bk.NRdPmax) < fabs(dP)) bk.NRdPmax = dP;
3887 
3888  // rockVp[n] = rockVntg[n] * (1 + rockC1 * (P[n] - rockPref));
3889  // cout << scientific << setprecision(6) << dP << " " << n << endl;
3890 
3891  dSmax = 0;
3892  chop = 1;
3893 
3894  const USI cNp = bk.phaseNum[n];
3895  const USI len = bk.bRowSizedSdP[n];
3896  Dcopy(len, &dtmp[0], &bk.res_n[n * row]);
3897  DaAxpby(len, col, 1.0, &bk.dSec_dPri[n * bk.maxLendSdP], u.data() + n * col,
3898  1.0, dtmp.data());
3899 
3900  // Calculate dSmax
3901  USI js = 0;
3902  for (USI j = 0; j < np; j++) {
3903  n_np_j = n * np + j;
3904  if (bk.pSderExist[n_np_j]) {
3905  dSmax = max(fabs(dtmp[js]), dSmax);
3906  js++;
3907  }
3908  }
3909 
3910  // Calculate chop
3911  if (dSmax > dSmaxlim) {
3912  chop = min(chop, dSmaxlim / dSmax);
3913  dSmax = dSmaxlim;
3914  }
3915  js = 0;
3916  for (USI j = 0; j < np; j++) {
3917  n_np_j = n * np + j;
3918  if (bk.pSderExist[n_np_j]) {
3919  if (fabs(bk.S[n_np_j] - scm[j]) > TINY &&
3920  (bk.S[n_np_j] - scm[j]) / (chop * dtmp[js]) < 0)
3921  chop *= min(1.0, -((bk.S[n_np_j] - scm[j]) / (chop * dtmp[js])));
3922  if (bk.S[n_np_j] + chop * dtmp[js] < 0)
3923  chop *= min(1.0, -bk.S[n_np_j] / (chop * dtmp[js]));
3924  js++;
3925  }
3926  }
3927 
3928  // Calculate S, phaseExist, xij, nj
3929  fill(tmpNij.begin(), tmpNij.end(), 0.0);
3930  // fill(Ni.begin() + n * nc, Ni.begin() + n * nc + nc, 0.0);
3931 
3932  js = 0;
3933  USI jx = cNp;
3934  for (USI j = 0; j < np; j++) {
3935  n_np_j = n * np + j;
3936  if (bk.pSderExist[n_np_j]) {
3937  bk.dSNRP[n_np_j] = chop * dtmp[js];
3938  js++;
3939 
3940  // S[n_np_j] += chop * dtmp[js];
3941  // if (S[n_np_j] < TINY) {
3942  // pSderExist[n_np_j] = OCP_FALSE;
3943  // }
3944  // js++;
3945  Daxpy(nc, bk.nj[n_np_j], &bk.xij[n_np_j * nc], &tmpNij[j * nc]);
3946  for (USI i = 0; i < bk.pVnumCom[j]; i++) {
3947  tmpNij[j * nc + i] += chop * dtmp[jx + i];
3948  }
3949  jx += bk.pVnumCom[j];
3950  bk.nj[n_np_j] = Dnorm1(nc, &tmpNij[j * nc]);
3951  for (USI i = 0; i < nc; i++) {
3952  bk.xij[n_np_j * nc + i] = tmpNij[j * nc + i] / bk.nj[n_np_j];
3953  // Ni[n * nc + i] += tmpNij[j * nc + i];
3954  }
3955  }
3956  }
3957 
3958  // Vf correction
3959  /* OCP_DBL dVmin = 100 * rockVp[n];
3960  USI index = 0;
3961  for (USI j = 0; j < numPhase; j++) {
3962  if (phaseExist[n * numPhase + j]) {
3963  tmpxi[j] = flashCal[PVTNUM[n]]->XiPhase(P[n], T, &xij[(n *
3964  numPhase + j) * nc]); tmpVf[j] = nj[n * numPhase + j] / (S[n * numPhase +
3965  j] * tmpxi[j]); if (fabs(tmpVf[j] - rockVp[n]) < dVmin) { dVmin =
3966  fabs(tmpVf[j] - rockVp[n]); index = j;
3967  }
3968  }
3969  } */
3970  // for (USI j = 0; j < numPhase; j++) {
3971  // if (phaseExist[n * numPhase + j]) {
3972  // nj[n * numPhase + j] *= (tmpVf[index] / tmpVf[j]);
3973  // for (USI i = 0; i < nc; i++) {
3974  // Ni[n * nc + i] += nj[n * numPhase + j] * xij[(n * numPhase +
3975  // j) * nc + i];
3976  // }
3977  // }
3978  // }
3979 
3980  bk.NRstep[n] = chop;
3981  for (USI i = 0; i < nc; i++) {
3982  bk.dNNR[n * nc + i] = u[n * col + 1 + i] * chop;
3983  if (fabs(bk.NRdNmax) < fabs(bk.dNNR[n * nc + i]) / bk.Nt[n])
3984  bk.NRdNmax = bk.dNNR[n * nc + i] / bk.Nt[n];
3985  bk.Ni[n * nc + i] += bk.dNNR[n * nc + i];
3986  }
3987  }
3988 
3989  // Well
3990  OCP_USI wId = nb * col;
3991  for (auto& wl : rs.allWells.wells) {
3992  if (wl.IsOpen()) {
3993  wl.SetBHP(wl.BHP() + u[wId]);
3994  wId += col;
3995  }
3996  }
3997 }
3998 
3999 void IsoT_FIMn::ResetToLastTimeStep(Reservoir& rs, OCPControl& ctrl)
4000 {
4001  rs.bulk.nj = rs.bulk.lnj;
4002  rs.bulk.res_n = rs.bulk.lres_n;
4003  rs.bulk.resPc = rs.bulk.lresPc;
4004 
4005  IsoT_FIM::ResetToLastTimeStep(rs, ctrl);
4006 }
4007 
4008 void IsoT_FIMn::UpdateLastTimeStep(Reservoir& rs) const
4009 {
4010  rs.bulk.lnj = rs.bulk.nj;
4011  rs.bulk.lres_n = rs.bulk.res_n;
4012  rs.bulk.lresPc = rs.bulk.resPc;
4013  IsoT_FIM::UpdateLastTimeStep(rs);
4014 }
4015 
4017 // IsoT_AIMc
4019 
4020 void IsoT_AIMc::Setup(Reservoir& rs, LinearSystem& ls, const OCPControl& ctrl)
4021 {
4022  // Allocate Bulk and BulkConn Memory
4023  AllocateReservoir(rs);
4024  // Allocate memory for internal matrix structure
4025  IsoT_FIM::AllocateLinearSystem(ls, rs, ctrl);
4026 }
4027 
4028 void IsoT_AIMc::InitReservoir(Reservoir& rs) const
4029 {
4030  rs.bulk.InitPTSw(50);
4031 
4032  InitRock(rs.bulk);
4033  CalRock(rs.bulk);
4034 
4035  IsoT_IMPEC::InitFlash(rs.bulk);
4036  IsoT_IMPEC::CalKrPc(rs.bulk);
4037 
4038  rs.allWells.InitBHP(rs.bulk);
4039 
4040  UpdateLastTimeStep(rs);
4041 }
4042 
4043 void IsoT_AIMc::Prepare(Reservoir& rs, const OCP_DBL& dt)
4044 {
4045  rs.allWells.PrepareWell(rs.bulk);
4046  CalRes(rs, dt, OCP_TRUE);
4047 
4048  // Set FIM Bulk
4049  rs.CalCFL(dt);
4050  rs.allWells.SetupWellBulk(rs.bulk);
4051  SetFIMBulk(rs);
4052  // Calculate FIM Bulk properties
4053  CalFlashI(rs.bulk);
4054  CalKrPcI(rs.bulk);
4055 
4056  UpdateLastTimeStep(rs);
4057  // rs.bulk.ShowFIMBulk(OCP_FALSE);
4058 }
4059 
4060 void IsoT_AIMc::AssembleMat(LinearSystem& ls,
4061  const Reservoir& rs,
4062  const OCP_DBL& dt) const
4063 {
4064  AssembleMatBulks(ls, rs, dt);
4065  IsoT_FIM::AssembleMatWellsNew(ls, rs, dt);
4066  ls.AssembleRhsCopy(rs.bulk.res.resAbs);
4067 }
4068 
4069 void IsoT_AIMc::SolveLinearSystem(LinearSystem& ls, Reservoir& rs, OCPControl& ctrl)
4070 {
4071 #ifdef DEBUG
4072  ls.CheckEquation();
4073 #endif // DEBUG
4074 
4075  ls.AssembleMatLinearSolver();
4076 
4077  GetWallTime Timer;
4078  Timer.Start();
4079  int status = ls.Solve();
4080  if (status < 0) {
4081  status = ls.GetNumIters();
4082  }
4083 
4084 #ifdef DEBUG
4085  // ls.OutputLinearSystem("testA_AIMc.out", "testb_AIMc.out");
4086  // ls.OutputSolution("testx_AIMc.out");
4087  // ls.CheckSolution();
4088 #endif // DEBUG
4089 
4090  ctrl.RecordTimeLS(Timer.Stop() / 1000);
4091  ctrl.UpdateIterLS(status);
4092  ctrl.UpdateIterNR();
4093 
4094  GetSolution(rs, ls.GetSolution(), ctrl);
4095  ls.ClearData();
4096 }
4097 
4098 OCP_BOOL IsoT_AIMc::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
4099 {
4100  OCP_DBL& dt = ctrl.current_dt;
4101 
4102  // First check: Ni check and bulk Pressure check
4103  if (!ctrl.Check(rs, {"BulkNi", "BulkP"})) {
4104  ResetToLastTimeStep(rs, ctrl);
4105  cout << "Cut time step size and repeat! current dt = " << fixed
4106  << setprecision(3) << dt << " days\n";
4107  return OCP_FALSE;
4108  }
4109 
4110  CalFlashI(rs.bulk);
4111  CalFlashEp(rs.bulk);
4112  CalKrPcI(rs.bulk);
4113 
4114  CalRock(rs.bulk);
4115 
4116  rs.allWells.CalTrans(rs.bulk);
4117  rs.allWells.CalFlux(rs.bulk);
4118 
4119  CalRes(rs, dt, OCP_FALSE);
4120  return OCP_TRUE;
4121 }
4122 
4123 OCP_BOOL IsoT_AIMc::FinishNR(Reservoir& rs, OCPControl& ctrl)
4124 {
4125  OCP_USI dSn;
4126  const OCP_DBL NRdSmax = rs.GetNRdSmax(dSn);
4127  const OCP_DBL NRdPmax = rs.GetNRdPmax();
4128  // const OCP_DBL NRdNmax = rs.GetNRdNmax();
4129 
4130 #ifdef DEBUG
4131  // cout << "### DEBUG: Residuals = " << setprecision(3) << scientific
4132  // << resAIMc.maxRelRes0_V << " " << resAIMc.maxRelRes_V << " "
4133  // << resAIMc.maxRelRes_N << " " << NRdPmax << " " << NRdSmax << endl;
4134 #endif
4135 
4136  if (((rs.bulk.res.maxRelRes_V <= rs.bulk.res.maxRelRes0_V * ctrl.ctrlNR.NRtol ||
4137  rs.bulk.res.maxRelRes_V <= ctrl.ctrlNR.NRtol ||
4138  rs.bulk.res.maxRelRes_N <= ctrl.ctrlNR.NRtol) &&
4139  rs.bulk.res.maxWellRelRes_mol <= ctrl.ctrlNR.NRtol) ||
4140  (fabs(NRdPmax) <= ctrl.ctrlNR.NRdPmin &&
4141  fabs(NRdSmax) <= ctrl.ctrlNR.NRdSmin)) {
4142 
4143  if (!ctrl.Check(rs, {"WellP"})) {
4144  ResetToLastTimeStep(rs, ctrl);
4145  return OCP_FALSE;
4146  } else {
4147  CalFlashEa(rs.bulk);
4148  CalKrPcE(rs.bulk);
4149  return OCP_TRUE;
4150  }
4151 
4152  } else if (ctrl.iterNR > ctrl.ctrlNR.maxNRiter) {
4153  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
4154  ResetToLastTimeStep(rs, ctrl);
4155  cout << "### WARNING: NR not fully converged! Cut time step size and repeat! "
4156  "current dt = "
4157  << fixed << setprecision(3) << ctrl.current_dt << " days\n";
4158  return OCP_FALSE;
4159  } else {
4160  return OCP_FALSE;
4161  }
4162 }
4163 
4165 void IsoT_AIMc::FinishStep(Reservoir& rs, OCPControl& ctrl) const
4166 {
4167  rs.CalIPRT(ctrl.GetCurDt());
4168  rs.CalMaxChange();
4169  UpdateLastTimeStep(rs);
4170  ctrl.CalNextTimeStep(rs, {"dP", "dS", "iter"});
4171 }
4172 
4174 void IsoT_AIMc::AllocateReservoir(Reservoir& rs)
4175 {
4176  IsoT_FIM::AllocateReservoir(rs);
4177 
4178  Bulk& bk = rs.bulk;
4179  const OCP_USI nb = bk.numBulk;
4180  const USI np = bk.numPhase;
4181  const USI nc = bk.numCom;
4182 
4183  bk.vj.resize(nb * np);
4184  bk.lvj.resize(nb * np);
4185 
4186  bk.xijNR.resize(nb * np * nc);
4187  bk.cfl.resize(nb * np);
4188  bk.bulkTypeAIM.Setup(nb);
4189 }
4190 
4191 void IsoT_AIMc::SetFIMBulk(Reservoir& rs)
4192 {
4193  Bulk& bk = rs.bulk;
4194  const BulkConn& conn = rs.conn;
4195  const OCP_USI nb = bk.numBulk;
4196  const USI np = bk.numPhase;
4197  const USI nc = bk.numCom;
4198 
4199  bk.bulkTypeAIM.Init();
4200 
4201  OCP_USI bIdp, bIdc;
4202  OCP_BOOL flag;
4203 
4204  for (OCP_USI n = 0; n < nb; n++) {
4205  bIdp = n * np;
4206  bIdc = n * nc;
4207  flag = OCP_FALSE;
4208  // CFL
4209  for (USI j = 0; j < np; j++) {
4210  if (bk.cfl[bIdp + j] > 0.8) {
4211  flag = OCP_TRUE;
4212  break;
4213  }
4214  }
4215  // Volume error
4216  if (!flag) {
4217  if ((fabs(bk.vf[n] - bk.rockVp[n]) / bk.rockVp[n]) > 1E-3) {
4218  flag = OCP_TRUE;
4219  }
4220  }
4221 
4222  // NR Step
4223  if (!flag && OCP_FALSE) {
4224  // dP
4225  if (fabs(bk.dPNR[n] / bk.P[n]) > 1E-3) {
4226  flag = OCP_TRUE;
4227  }
4228  // dNi
4229  if (!flag) {
4230  for (USI i = 0; i < bk.numCom; i++) {
4231  if (fabs(bk.dNNR[bIdc + i] / bk.Ni[bIdc + i]) > 1E-3) {
4232  flag = OCP_TRUE;
4233  break;
4234  }
4235  }
4236  }
4237  }
4238 
4239  if (flag) {
4240  // find it
4241  // bk.map_Bulk2FIM[n] = 1;
4242  for (auto& v : conn.neighbor[n]) {
4243  // n is included also
4244  bk.bulkTypeAIM.SetBulkType(v, 1);
4245  }
4246  }
4247  }
4248 
4249  // add WellBulk's 2-neighbor as Implicit bulk
4250  for (auto& p : bk.wellBulkId) {
4251  for (auto& v : conn.neighbor[p]) {
4252  for (auto& v1 : conn.neighbor[v]) bk.bulkTypeAIM.SetBulkType(v1, 1);
4253  }
4254  }
4255 }
4256 
4257 void IsoT_AIMc::CalFlashEp(Bulk& bk)
4258 {
4259  const OCP_USI nb = bk.numBulk;
4260  const USI np = bk.numPhase;
4261  const USI nc = bk.numCom;
4262 
4263  for (OCP_USI n = 0; n < nb; n++) {
4264  if (bk.bulkTypeAIM.IfIMPECbulk(n)) {
4265  // Explicit bulk
4266 
4267  bk.flashCal[bk.PVTNUM[n]]->FlashIMPEC(bk.P[n], bk.T[n], &bk.Ni[n * nc],
4268  bk.phaseNum[n],
4269  &bk.xijNR[n * np * nc], n);
4270  // bk.PassFlashValueAIMcEp(n);
4271  PassFlashValueEp(bk, n);
4272  }
4273  }
4274 }
4275 
4276 void IsoT_AIMc::CalFlashEa(Bulk& bk)
4277 {
4278  const OCP_USI nb = bk.numBulk;
4279  const USI np = bk.numPhase;
4280  const USI nc = bk.numCom;
4281 
4282  for (OCP_USI n = 0; n < nb; n++) {
4283  if (bk.bulkTypeAIM.IfIMPECbulk(n)) {
4284  // Explicit bulk
4285 
4286  bk.flashCal[bk.PVTNUM[n]]->FlashIMPEC(bk.P[n], bk.T[n], &bk.Ni[n * nc],
4287  bk.phaseNum[n], &bk.xij[n * np * nc],
4288  n);
4289  // bk.PassFlashValueAIMcEa(n);
4290 
4291  IsoT_IMPEC::PassFlashValue(bk, n);
4292  }
4293  }
4294 }
4295 
4296 void IsoT_AIMc::CalFlashI(Bulk& bk)
4297 {
4298  const OCP_USI nb = bk.numBulk;
4299  const USI np = bk.numPhase;
4300  const USI nc = bk.numCom;
4301 
4302  for (OCP_USI n = 0; n < nb; n++) {
4303  if (bk.bulkTypeAIM.IfFIMbulk(n)) {
4304  // Implicit bulk
4305 
4306  bk.flashCal[bk.PVTNUM[n]]->FlashFIM(bk.P[n], bk.T[n], &bk.Ni[n * nc],
4307  &bk.S[n * np], bk.phaseNum[n],
4308  &bk.xij[n * np * nc], n);
4309  IsoT_FIM::PassFlashValue(bk, n);
4310  for (USI j = 0; j < np; j++) {
4311  bk.vj[n * np + j] = bk.vf[n] * bk.S[n * np + j];
4312  }
4313  }
4314  }
4315 }
4316 
4317 void IsoT_AIMc::PassFlashValueEp(Bulk& bk, const OCP_USI& n)
4318 {
4319  // only var about volume needs, some flash var also
4320  OCP_FUNCNAME;
4321 
4322  const USI np = bk.numPhase;
4323  const USI nc = bk.numCom;
4324  const OCP_USI bIdp = n * np;
4325  const USI pvtnum = bk.PVTNUM[n];
4326 
4327  bk.Nt[n] = bk.flashCal[pvtnum]->GetNt();
4328  bk.vf[n] = bk.flashCal[pvtnum]->GetVf();
4329  bk.vfP[n] = bk.flashCal[pvtnum]->GetVfP();
4330  for (USI i = 0; i < nc; i++) {
4331  bk.vfi[n * nc + i] = bk.flashCal[pvtnum]->GetVfi(i);
4332  }
4333 
4334  bk.phaseNum[n] = 0;
4335  for (USI j = 0; j < np; j++) {
4336  if (bk.flashCal[pvtnum]->GetPhaseExist(j)) {
4337  bk.phaseNum[n]++;
4338 
4339  // IMPORTANT -- need for next Flash
4340  // But xij in nonlinear equations has been modified
4341  for (USI i = 0; i < nc; i++) {
4342  bk.xijNR[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetXij(j, i);
4343  }
4344  }
4345  }
4346 }
4347 
4348 void IsoT_AIMc::CalKrPcE(Bulk& bk)
4349 {
4350  const OCP_USI nb = bk.numBulk;
4351  const USI np = bk.numPhase;
4352 
4353  for (OCP_USI n = 0; n < nb; n++) {
4354  if (bk.bulkTypeAIM.IfIMPECbulk(n)) {
4355  // Explicit bulk
4356  const OCP_USI bId = n * np;
4357  bk.flow[bk.SATNUM[n]]->CalKrPc(&bk.S[bId], &bk.kr[bId], &bk.Pc[bId], n);
4358  for (USI j = 0; j < np; j++) bk.Pj[bId + j] = bk.P[n] + bk.Pc[bId + j];
4359  }
4360  }
4361 }
4362 
4363 void IsoT_AIMc::CalKrPcI(Bulk& bk)
4364 {
4365  const OCP_USI nb = bk.numBulk;
4366  const USI np = bk.numPhase;
4367 
4368  for (OCP_USI n = 0; n < nb; n++) {
4369  if (bk.bulkTypeAIM.IfFIMbulk(n)) {
4370  // Implicit bulk
4371  const OCP_USI bId = n * np;
4372  bk.flow[bk.SATNUM[n]]->CalKrPcDeriv(&bk.S[bId], &bk.kr[bId], &bk.Pc[bId],
4373  &bk.dKr_dS[bId * np],
4374  &bk.dPcj_dS[bId * np], n);
4375  for (USI j = 0; j < np; j++) bk.Pj[bId + j] = bk.P[n] + bk.Pc[bId + j];
4376  }
4377  }
4378 }
4379 
4380 void IsoT_AIMc::AssembleMatBulks(LinearSystem& ls,
4381  const Reservoir& rs,
4382  const OCP_DBL& dt) const
4383 {
4384  const Bulk& bk = rs.bulk;
4385  const BulkConn& conn = rs.conn;
4386  const OCP_USI nb = bk.numBulk;
4387  const USI np = bk.numPhase;
4388  const USI nc = bk.numCom;
4389  const USI ncol = nc + 1;
4390  const USI ncol2 = np * nc + np;
4391  const USI bsize = ncol * ncol;
4392  const USI bsize2 = ncol * ncol2;
4393  const USI lendSdP = bk.maxLendSdP;
4394 
4395  ls.AddDim(nb);
4396 
4397  vector<OCP_DBL> bmat(bsize, 0);
4398  // Accumulation term
4399  for (USI i = 1; i < nc + 1; i++) {
4400  bmat[i * ncol + i] = 1;
4401  }
4402  for (OCP_USI n = 0; n < nb; n++) {
4403  bmat[0] = bk.v[n] * bk.poroP[n] - bk.vfP[n];
4404  for (USI i = 0; i < nc; i++) {
4405  bmat[i + 1] = -bk.vfi[n * nc + i];
4406  }
4407  ls.NewDiag(n, bmat);
4408  }
4409 
4410  // flux term
4411  OCP_DBL Akd;
4412  OCP_DBL transJ, transIJ;
4413  vector<OCP_DBL> dFdXpB(bsize, 0);
4414  vector<OCP_DBL> dFdXpE(bsize, 0);
4415  vector<OCP_DBL> dFdXsB(bsize2, 0);
4416  vector<OCP_DBL> dFdXsE(bsize2, 0);
4417  OCP_DBL* dFdXpI;
4418  OCP_DBL* dFdXsI;
4419  vector<OCP_BOOL> phaseExistB(np, OCP_FALSE);
4420  vector<OCP_BOOL> phaseExistE(np, OCP_FALSE);
4421  vector<OCP_BOOL> phaseExistI(np, OCP_FALSE);
4422  OCP_BOOL phaseExistU;
4423  vector<OCP_BOOL> phasedS_B(np, OCP_FALSE);
4424  vector<OCP_BOOL> phasedS_E(np, OCP_FALSE);
4425  vector<OCP_BOOL> phasedS_I(np, OCP_FALSE);
4426  vector<USI> pVnumComB(np, 0);
4427  vector<USI> pVnumComE(np, 0);
4428  vector<USI> pVnumComI(np, 0);
4429  OCP_USI ncolB, ncolE, ncolI;
4430 
4431  OCP_USI bId, eId, uId;
4432  OCP_USI bId_np_j, eId_np_j, uId_np_j, ibId_np_j;
4433  OCP_DBL kr, mu, xi, xij, rhoP, xiP, muP, rhox, xix, mux;
4434  OCP_DBL dP, dGamma;
4435  OCP_DBL tmp;
4436  OCP_BOOL bIdFIM, eIdFIM, uIdFIM;
4437 
4438  for (OCP_USI c = 0; c < conn.numConn; c++) {
4439  bId = conn.iteratorConn[c].BId();
4440  eId = conn.iteratorConn[c].EId();
4441  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
4442  dGamma = GRAVITY_FACTOR * (bk.depth[bId] - bk.depth[eId]);
4443  bIdFIM = eIdFIM = OCP_FALSE;
4444  if (bk.bulkTypeAIM.IfFIMbulk(bId)) bIdFIM = OCP_TRUE;
4445  if (bk.bulkTypeAIM.IfFIMbulk(eId)) eIdFIM = OCP_TRUE;
4446 
4447  if (!bIdFIM && !eIdFIM) {
4448  // both are explicit
4449  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
4450  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
4451  for (USI j = 0; j < np; j++) {
4452  phaseExistB[j] = bk.phaseExist[bId * np + j];
4453  phaseExistE[j] = bk.phaseExist[eId * np + j];
4454  uId = conn.upblock[c * np + j];
4455  phaseExistU = (uId == bId ? phaseExistB[j] : phaseExistE[j]);
4456  if (!phaseExistU) {
4457  continue;
4458  }
4459  bId_np_j = bId * np + j;
4460  eId_np_j = eId * np + j;
4461  uId_np_j = uId * np + j;
4462  xi = bk.xi[uId_np_j];
4463  kr = bk.kr[uId_np_j];
4464  mu = bk.mu[uId_np_j];
4465  transJ = Akd * xi * kr / mu;
4466  for (USI i = 0; i < nc; i++) {
4467  xij = bk.xij[uId_np_j * nc + i];
4468  transIJ = xij * transJ;
4469  // Pressure -- Primary var
4470  dFdXpB[(i + 1) * ncol] += transIJ;
4471  dFdXpE[(i + 1) * ncol] -= transIJ;
4472  }
4473  }
4474  } else if ((bIdFIM && !eIdFIM) || (!bIdFIM && eIdFIM)) {
4475 
4476  // one is explicit, one is implicit
4477  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
4478  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
4479  fill(dFdXsB.begin(), dFdXsB.end(), 0.0);
4480  fill(dFdXsE.begin(), dFdXsE.end(), 0.0);
4481 
4482  ncolI = 0;
4483  USI jxI = 0;
4484  if (bIdFIM) {
4485  for (USI j = 0; j < np; j++) {
4486  phaseExistI[j] = bk.phaseExist[bId * np + j];
4487  phasedS_I[j] = bk.pSderExist[bId * np + j];
4488  if (phasedS_I[j]) jxI++;
4489  pVnumComI[j] = bk.pVnumCom[bId * np + j];
4490  ncolI += pVnumComI[j];
4491  }
4492  dFdXpI = &dFdXpB[0];
4493  dFdXsI = &dFdXsB[0];
4494  } else {
4495  for (USI j = 0; j < np; j++) {
4496  phaseExistI[j] = bk.phaseExist[eId * np + j];
4497  phasedS_I[j] = bk.pSderExist[eId * np + j];
4498  if (phasedS_I[j]) jxI++;
4499  pVnumComI[j] = bk.pVnumCom[eId * np + j];
4500  ncolI += pVnumComI[j];
4501  }
4502  dFdXpI = &dFdXpE[0];
4503  dFdXsI = &dFdXsE[0];
4504  }
4505  ncolI += jxI;
4506 
4507  for (USI j = 0; j < np; j++) {
4508  phaseExistB[j] = bk.phaseExist[bId * np + j];
4509  phaseExistE[j] = bk.phaseExist[eId * np + j];
4510  uId = conn.upblock[c * np + j];
4511  phaseExistU = (uId == bId ? phaseExistB[j] : phaseExistE[j]);
4512  if (!phaseExistU) {
4513  continue;
4514  }
4515 
4516  bId_np_j = bId * np + j;
4517  eId_np_j = eId * np + j;
4518  uId_np_j = uId * np + j;
4519  dP = bk.Pj[bId_np_j] - bk.Pj[eId_np_j] -
4520  conn.upblock_Rho[c * np + j] * dGamma;
4521  xi = bk.xi[uId_np_j];
4522  kr = bk.kr[uId_np_j];
4523  mu = bk.mu[uId_np_j];
4524  transJ = Akd * kr / mu;
4525  uIdFIM = (uId == bId ? bIdFIM : eIdFIM);
4526 
4527  if (bIdFIM)
4528  ibId_np_j = bId_np_j;
4529  else
4530  ibId_np_j = eId_np_j;
4531 
4532  if (uIdFIM) {
4533  // upblock is implicit
4534  OCP_DBL rhoWgt = 1.0;
4535  if (phaseExistB[j] && phaseExistE[j]) rhoWgt = 0.5;
4536 
4537  muP = bk.muP[ibId_np_j];
4538  xiP = bk.xiP[ibId_np_j];
4539  rhoP = bk.rhoP[ibId_np_j];
4540  for (USI i = 0; i < nc; i++) {
4541  xij = bk.xij[ibId_np_j * nc + i];
4542  transIJ = xij * xi * transJ;
4543 
4544  // Pressure -- Primary var
4545  dFdXpB[(i + 1) * ncol] += transIJ;
4546  dFdXpE[(i + 1) * ncol] -= transIJ;
4547 
4548  tmp = transIJ * (-rhoP * dGamma) * rhoWgt;
4549  tmp += xij * transJ * xiP * dP;
4550  tmp += -transIJ * muP / mu * dP;
4551  dFdXpI[(i + 1) * ncol] += tmp;
4552 
4553  USI j1S = 0;
4554  // Saturation -- Second var
4555  for (USI j1 = 0; j1 < np; j1++) {
4556  if (phasedS_I[j1]) {
4557  dFdXsI[(i + 1) * ncolI + j1S] +=
4558  transIJ * bk.dPcj_dS[ibId_np_j * np + j1];
4559  tmp = Akd * xij * xi / mu *
4560  bk.dKr_dS[ibId_np_j * np + j1] * dP;
4561  dFdXsI[(i + 1) * ncolI + j1S] += tmp;
4562  j1S++;
4563  }
4564  }
4565  // Cij
4566  for (USI k = 0; k < pVnumComI[j]; k++) {
4567  rhox = bk.rhox[ibId_np_j * nc + k] * rhoWgt;
4568  xix = bk.xix[ibId_np_j * nc + k];
4569  mux = bk.mux[ibId_np_j * nc + k];
4570  tmp = -transIJ * rhox * dGamma;
4571  tmp += xij * transJ * xix * dP;
4572  tmp += -transIJ * mux / mu * dP;
4573  dFdXsI[(i + 1) * ncolI + jxI + k] += tmp;
4574  }
4575  // WARNING !!!
4576  if (i < pVnumComI[j])
4577  dFdXsI[(i + 1) * ncolI + jxI + i] += xi * transJ * dP;
4578  }
4579  jxI += pVnumComI[j];
4580  } else {
4581  // downblock is implicit
4582  OCP_DBL rhoWgt = 0;
4583  if (phaseExistB[j] && phaseExistE[j]) rhoWgt = 0.5;
4584  rhoP = bk.rhoP[ibId_np_j];
4585 
4586  for (USI i = 0; i < nc; i++) {
4587  xij = bk.xij[uId_np_j * nc + i];
4588  transIJ = xij * xi * transJ;
4589 
4590  // Pressure -- Primary var
4591  dFdXpB[(i + 1) * ncol] += transIJ;
4592  dFdXpE[(i + 1) * ncol] -= transIJ;
4593 
4594  tmp = transIJ * (-rhoP * dGamma) * rhoWgt;
4595  dFdXpI[(i + 1) * ncol] += tmp;
4596 
4597  USI j1S = 0;
4598  // Saturation -- Second var
4599  for (USI j1 = 0; j1 < np; j1++) {
4600  if (phasedS_I[j1]) {
4601  dFdXsI[(i + 1) * ncolI + j1S] +=
4602  transIJ * bk.dPcj_dS[ibId_np_j * np + j1];
4603  j1S++;
4604  }
4605  }
4606  // Cij
4607  for (USI k = 0; k < pVnumComI[j]; k++) {
4608  rhox = bk.rhox[ibId_np_j * nc + k] * rhoWgt;
4609  tmp = -transIJ * rhox * dGamma;
4610  dFdXsI[(i + 1) * ncolI + jxI + k] += tmp;
4611  }
4612  }
4613  jxI += pVnumComI[j];
4614  }
4615  }
4616  } else {
4617  // both are implicit
4618  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
4619  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
4620  fill(dFdXsB.begin(), dFdXsB.end(), 0.0);
4621  fill(dFdXsE.begin(), dFdXsE.end(), 0.0);
4622 
4623  USI jxB = 0;
4624  USI jxE = 0;
4625  ncolB = 0;
4626  ncolE = 0;
4627 
4628  for (USI j = 0; j < np; j++) {
4629  phaseExistB[j] = bk.phaseExist[bId * np + j];
4630  phaseExistE[j] = bk.phaseExist[eId * np + j];
4631  phasedS_B[j] = bk.pSderExist[bId * np + j];
4632  phasedS_E[j] = bk.pSderExist[eId * np + j];
4633  if (phasedS_B[j]) jxB++;
4634  if (phasedS_E[j]) jxE++;
4635  pVnumComB[j] = bk.pVnumCom[bId * np + j];
4636  pVnumComE[j] = bk.pVnumCom[eId * np + j];
4637  ncolB += pVnumComB[j];
4638  ncolE += pVnumComE[j];
4639  }
4640  ncolB += jxB;
4641  ncolE += jxE;
4642 
4643  for (USI j = 0; j < np; j++) {
4644  uId = conn.upblock[c * np + j];
4645 
4646  phaseExistU = (uId == bId ? phaseExistB[j] : phaseExistE[j]);
4647  if (!phaseExistU) {
4648  jxB += pVnumComB[j];
4649  jxE += pVnumComE[j];
4650  continue;
4651  }
4652 
4653  bId_np_j = bId * np + j;
4654  eId_np_j = eId * np + j;
4655  uId_np_j = uId * np + j;
4656  dP = bk.Pj[bId_np_j] - bk.Pj[eId_np_j] -
4657  conn.upblock_Rho[c * np + j] * dGamma;
4658  xi = bk.xi[uId_np_j];
4659  kr = bk.kr[uId_np_j];
4660  mu = bk.mu[uId_np_j];
4661  muP = bk.muP[uId_np_j];
4662  xiP = bk.xiP[uId_np_j];
4663  rhoP = bk.rhoP[uId_np_j];
4664  transJ = Akd * kr / mu;
4665 
4666  for (USI i = 0; i < nc; i++) {
4667  xij = bk.xij[uId_np_j * nc + i];
4668  transIJ = xij * xi * transJ;
4669 
4670  // Pressure -- Primary var
4671  dFdXpB[(i + 1) * ncol] += transIJ;
4672  dFdXpE[(i + 1) * ncol] -= transIJ;
4673 
4674  tmp = xij * transJ * xiP * dP;
4675  tmp += -transIJ * muP / mu * dP;
4676  if (!phaseExistE[j]) {
4677  tmp += transIJ * (-rhoP * dGamma);
4678  dFdXpB[(i + 1) * ncol] += tmp;
4679  } else if (!phaseExistB[j]) {
4680  tmp += transIJ * (-rhoP * dGamma);
4681  dFdXpE[(i + 1) * ncol] += tmp;
4682  } else {
4683  dFdXpB[(i + 1) * ncol] +=
4684  transIJ * (-bk.rhoP[bId_np_j] * dGamma) / 2;
4685  dFdXpE[(i + 1) * ncol] +=
4686  transIJ * (-bk.rhoP[eId_np_j] * dGamma) / 2;
4687  if (bId == uId) {
4688  dFdXpB[(i + 1) * ncol] += tmp;
4689  } else {
4690  dFdXpE[(i + 1) * ncol] += tmp;
4691  }
4692  }
4693 
4694  // Second var
4695  USI j1SB = 0;
4696  USI j1SE = 0;
4697  if (bId == uId) {
4698  // Saturation
4699  for (USI j1 = 0; j1 < np; j1++) {
4700  if (phasedS_B[j1]) {
4701  dFdXsB[(i + 1) * ncolB + j1SB] +=
4702  transIJ * bk.dPcj_dS[bId_np_j * np + j1];
4703  tmp = Akd * xij * xi / mu *
4704  bk.dKr_dS[uId_np_j * np + j1] * dP;
4705  dFdXsB[(i + 1) * ncolB + j1SB] += tmp;
4706  j1SB++;
4707  }
4708  if (phasedS_E[j1]) {
4709  dFdXsE[(i + 1) * ncolE + j1SE] -=
4710  transIJ * bk.dPcj_dS[eId_np_j * np + j1];
4711  j1SE++;
4712  }
4713  }
4714  // Cij
4715  if (!phaseExistE[j]) {
4716  for (USI k = 0; k < pVnumComB[j]; k++) {
4717  rhox = bk.rhox[uId_np_j * nc + k];
4718  xix = bk.xix[uId_np_j * nc + k];
4719  mux = bk.mux[uId_np_j * nc + k];
4720  tmp = -transIJ * rhox * dGamma;
4721  tmp += xij * transJ * xix * dP;
4722  tmp += -transIJ * mux / mu * dP;
4723  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
4724  }
4725  // WARNING !!!
4726  if (i < pVnumComB[j])
4727  dFdXsB[(i + 1) * ncolB + jxB + i] += xi * transJ * dP;
4728  } else {
4729  for (USI k = 0; k < pVnumComB[j]; k++) {
4730  rhox = bk.rhox[bId_np_j * nc + k] / 2;
4731  xix = bk.xix[uId_np_j * nc + k];
4732  mux = bk.mux[uId_np_j * nc + k];
4733  tmp = -transIJ * rhox * dGamma;
4734  tmp += xij * transJ * xix * dP;
4735  tmp += -transIJ * mux / mu * dP;
4736  dFdXsB[(i + 1) * ncolB + jxB + k] += tmp;
4737  dFdXsE[(i + 1) * ncolE + jxE + k] +=
4738  -transIJ * bk.rhox[eId_np_j * nc + k] / 2 * dGamma;
4739  }
4740  // WARNING !!!
4741  if (i < pVnumComB[j])
4742  dFdXsB[(i + 1) * ncolB + jxB + i] += xi * transJ * dP;
4743  }
4744  } else {
4745  // Saturation
4746  for (USI j1 = 0; j1 < np; j1++) {
4747  if (phasedS_B[j1]) {
4748  dFdXsB[(i + 1) * ncolB + j1SB] +=
4749  transIJ * bk.dPcj_dS[bId_np_j * np + j1];
4750  j1SB++;
4751  }
4752  if (phasedS_E[j1]) {
4753  dFdXsE[(i + 1) * ncolE + j1SE] -=
4754  transIJ * bk.dPcj_dS[eId_np_j * np + j1];
4755  tmp = Akd * xij * xi / mu *
4756  bk.dKr_dS[uId_np_j * np + j1] * dP;
4757  dFdXsE[(i + 1) * ncolE + j1SE] += tmp;
4758  j1SE++;
4759  }
4760  }
4761  // Cij
4762  if (!phaseExistB[j]) {
4763  for (USI k = 0; k < pVnumComE[j]; k++) {
4764  rhox = bk.rhox[uId_np_j * nc + k];
4765  xix = bk.xix[uId_np_j * nc + k];
4766  mux = bk.mux[uId_np_j * nc + k];
4767  tmp = -transIJ * rhox * dGamma;
4768  tmp += xij * transJ * xix * dP;
4769  tmp += -transIJ * mux / mu * dP;
4770  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
4771  }
4772  // WARNING !!!
4773  if (i < pVnumComE[j])
4774  dFdXsE[(i + 1) * ncolE + jxE + i] += xi * transJ * dP;
4775  } else {
4776  for (USI k = 0; k < pVnumComE[j]; k++) {
4777  rhox = bk.rhox[eId_np_j * nc + k] / 2;
4778  xix = bk.xix[uId_np_j * nc + k];
4779  mux = bk.mux[uId_np_j * nc + k];
4780  tmp = -transIJ * rhox * dGamma;
4781  tmp += xij * transJ * xix * dP;
4782  tmp += -transIJ * mux / mu * dP;
4783  dFdXsE[(i + 1) * ncolE + jxE + k] += tmp;
4784  dFdXsB[(i + 1) * ncolB + jxB + k] +=
4785  -transIJ * bk.rhox[bId_np_j * nc + k] / 2 * dGamma;
4786  }
4787  // WARNING !!!
4788  if (i < pVnumComE[j])
4789  dFdXsE[(i + 1) * ncolE + jxE + i] += xi * transJ * dP;
4790  }
4791  }
4792  }
4793  jxB += pVnumComB[j];
4794  jxE += pVnumComE[j];
4795  }
4796  }
4797 
4798  if (bIdFIM && !eIdFIM)
4799  ncolB = ncolI;
4800  else if (!bIdFIM && eIdFIM)
4801  ncolE = ncolI;
4802 
4803  // Assemble
4804  bmat = dFdXpB;
4805  if (bIdFIM) {
4806  DaABpbC(ncol, ncol, ncolB, 1, dFdXsB.data(), &bk.dSec_dPri[bId * lendSdP],
4807  1, bmat.data());
4808  }
4809  Dscalar(bsize, dt, bmat.data());
4810  // Begin - Begin -- add
4811  ls.AddDiag(bId, bmat);
4812  // End - Begin -- insert
4813  Dscalar(bsize, -1, bmat.data());
4814  ls.NewOffDiag(eId, bId, bmat);
4815 
4816 #ifdef OCP_NANCHECK
4817  if (!CheckNan(bmat.size(), &bmat[0])) {
4818  OCP_ABORT("INF or NAN in bmat !");
4819  }
4820 #endif
4821 
4822  // End
4823  bmat = dFdXpE;
4824  if (eIdFIM) {
4825  DaABpbC(ncol, ncol, ncolE, 1, dFdXsE.data(), &bk.dSec_dPri[eId * lendSdP],
4826  1, bmat.data());
4827  }
4828  Dscalar(bsize, dt, bmat.data());
4829  // Begin - End -- insert
4830  ls.NewOffDiag(bId, eId, bmat);
4831  // End - End -- add
4832  Dscalar(bsize, -1, bmat.data());
4833  ls.AddDiag(eId, bmat);
4834 
4835 #ifdef OCP_NANCHECK
4836  if (!CheckNan(bmat.size(), &bmat[0])) {
4837  OCP_ABORT("INF or NAN in bmat !");
4838  }
4839 #endif
4840  }
4841 }
4842 
4843 void IsoT_AIMc::GetSolution(Reservoir& rs,
4844  const vector<OCP_DBL>& u,
4845  const OCPControl& ctrl) const
4846 {
4847  // Bulk
4848  const OCP_DBL dSmaxlim = ctrl.ctrlNR.NRdSmax;
4849  // const OCP_DBL dPmaxlim = ctrl.ctrlNR.NRdPmax;
4850 
4851  Bulk& bk = rs.bulk;
4852  const OCP_USI nb = bk.numBulk;
4853  const USI np = bk.numPhase;
4854  const USI nc = bk.numCom;
4855  const USI row = np * (nc + 1);
4856  const USI col = nc + 1;
4857  vector<OCP_DBL> dtmp(row, 0);
4858  OCP_DBL chopmin = 1;
4859  OCP_DBL choptmp = 0;
4860  OCP_USI n_np_j;
4861 
4862  bk.dSNR = bk.S;
4863  bk.NRphaseNum = bk.phaseNum;
4864  bk.NRdPmax = 0;
4865  bk.NRdNmax = 0;
4866 
4867  for (OCP_USI n = 0; n < nb; n++) {
4868  if (bk.bulkTypeAIM.IfIMPECbulk(n)) {
4869  // IMPEC Bulk
4870  // Pressure
4871  OCP_DBL dP = u[n * col];
4872  bk.NRdPmax = max(bk.NRdPmax, fabs(dP));
4873  bk.P[n] += dP; // seems better
4874  bk.dPNR[n] = dP;
4875  bk.NRstep[n] = 1;
4876  // Ni
4877  for (USI i = 0; i < nc; i++) {
4878  bk.dNNR[n * nc + i] = u[n * col + 1 + i];
4879  bk.Ni[n * nc + i] += bk.dNNR[n * nc + i];
4880 
4881  // if (bk.Ni[n * nc + i] < 0 && bk.Ni[n * nc + i] > -1E-3) {
4882  // bk.Ni[n * nc + i] = 1E-20;
4883  // }
4884  }
4885  // Pj
4886  for (USI j = 0; j < np; j++) {
4887  bk.Pj[n * np + j] = bk.P[n] + bk.Pc[n * np + j];
4888  }
4889  continue;
4890  }
4891 
4892  chopmin = 1;
4893  // compute the chop
4894  fill(dtmp.begin(), dtmp.end(), 0.0);
4895  DaAxpby(bk.bRowSizedSdP[n], col, 1, &bk.dSec_dPri[n * bk.maxLendSdP],
4896  u.data() + n * col, 1, dtmp.data());
4897 
4898  USI js = 0;
4899  for (USI j = 0; j < np; j++) {
4900  if (!bk.pSderExist[n * np + j]) {
4901  continue;
4902  }
4903  n_np_j = n * np + j;
4904 
4905  choptmp = 1;
4906  if (fabs(dtmp[js]) > dSmaxlim) {
4907  choptmp = dSmaxlim / fabs(dtmp[js]);
4908  } else if (bk.S[n_np_j] + dtmp[js] < 0.0) {
4909  choptmp = 0.9 * bk.S[n_np_j] / fabs(dtmp[js]);
4910  }
4911 
4912  // if (fabs(S[n_np_j] - scm[j]) > TINY &&
4913  // (S[n_np_j] - scm[j]) / (choptmp * dtmp[js]) < 0)
4914  // choptmp *= min(1.0, -((S[n_np_j] - scm[j]) / (choptmp * dtmp[js])));
4915 
4916  chopmin = min(chopmin, choptmp);
4917  js++;
4918  }
4919 
4920  // dS
4921  js = 0;
4922  for (USI j = 0; j < np; j++) {
4923  if (!bk.pSderExist[n * np + j]) {
4924  bk.dSNRP[n * np + j] = 0;
4925  continue;
4926  }
4927  bk.dSNRP[n * np + j] = chopmin * dtmp[js];
4928  js++;
4929  }
4930 
4931  // dxij ---- Compositional model only
4932  if (bk.IfUseEoS()) {
4933  if (bk.phaseNum[n] >= 3) {
4934  OCP_USI bId = 0;
4935  for (USI j = 0; j < 2; j++) {
4936  bId = n * np * nc + j * nc;
4937  for (USI i = 0; i < bk.numComH; i++) {
4938  bk.xij[bId + i] += chopmin * dtmp[js];
4939  js++;
4940  }
4941  }
4942  }
4943  }
4944 
4945  // dP
4946  OCP_DBL dP = u[n * col];
4947  if (fabs(bk.NRdPmax) < fabs(dP)) bk.NRdPmax = dP;
4948  bk.P[n] += dP; // seems better
4949  bk.dPNR[n] = dP;
4950 
4951  // dNi
4952  bk.NRstep[n] = chopmin;
4953  for (USI i = 0; i < nc; i++) {
4954  bk.dNNR[n * nc + i] = u[n * col + 1 + i] * chopmin;
4955  if (fabs(bk.NRdNmax) < fabs(bk.dNNR[n * nc + i]) / bk.Nt[n])
4956  bk.NRdNmax = bk.dNNR[n * nc + i] / bk.Nt[n];
4957 
4958  bk.Ni[n * nc + i] += bk.dNNR[n * nc + i];
4959 
4960  // if (bk.Ni[n * nc + i] < 0 && bk.Ni[n * nc + i] > -1E-3) {
4961  // bk.Ni[n * nc + i] = 1E-20;
4962  // }
4963  }
4964  }
4965 
4966  // Well
4967  OCP_USI wId = nb * col;
4968  for (auto& wl : rs.allWells.wells) {
4969  if (wl.IsOpen()) {
4970  wl.SetBHP(wl.BHP() + u[wId]);
4971  wId += col;
4972  }
4973  }
4974 }
4975 
4976 void IsoT_AIMc::ResetToLastTimeStep(Reservoir& rs, OCPControl& ctrl)
4977 {
4978  rs.bulk.vj = rs.bulk.lvj;
4979  rs.bulk.xijNR = rs.bulk.lxij;
4980  IsoT_FIM::ResetToLastTimeStep(rs, ctrl);
4981 }
4982 
4983 void IsoT_AIMc::UpdateLastTimeStep(Reservoir& rs) const
4984 {
4985  IsoT_FIM::UpdateLastTimeStep(rs);
4986 
4987  rs.bulk.lvj = rs.bulk.vj;
4988  rs.bulk.xijNR = rs.bulk.xij;
4989 }
4990 
4991 /*----------------------------------------------------------------------------*/
4992 /* Brief Change History of This File */
4993 /*----------------------------------------------------------------------------*/
4994 /* Author Date Actions */
4995 /*----------------------------------------------------------------------------*/
4996 /* Shizhe Li Nov/01/2021 Create file */
4997 /* Chensong Zhang Jan/08/2022 Update output */
4998 /*----------------------------------------------------------------------------*/
DaABpbC
void DaABpbC(const int &m, const int &n, const int &k, const double &alpha, const double *A, const double *B, const double &beta, double *C)
Computes C' = alpha B'A' + beta C', all matrices are column-major.
Definition: DenseMat.cpp:75
Dscalar
void Dscalar(const int &n, const double &alpha, double *x)
Scales a vector by a constant.
Definition: DenseMat.cpp:61
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:203
CheckNan
bool CheckNan(const int &N, const T *x)
check NaN
Definition: DenseMat.hpp:209
Daxpy
void Daxpy(const int &n, const double &alpha, const double *x, double *y)
Constant times a vector plus a vector.
Definition: DenseMat.cpp:68
Dnorm1
double Dnorm1(const int &N, double *x)
Computes the L1-norm of a vector.
Definition: DenseMat.cpp:49
Dcopy
void Dcopy(const int &N, double *dst, const double *src)
Calculate the minimal eigenvalue for symmetric matrix with mkl lapack.
Definition: DenseMat.cpp:36
VECTORFASP
const USI VECTORFASP
Use vector linear solver in Fasp.
Definition: OCPConst.hpp:88
TINY
const OCP_DBL TINY
Small constant.
Definition: OCPConst.hpp:38
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:23
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:27
CONV1
const OCP_DBL CONV1
1 bbl = CONV1 ft3
Definition: OCPConst.hpp:63
GRAVITY_FACTOR
const OCP_DBL GRAVITY_FACTOR
0.00694444 ft2 psi / lb
Definition: OCPConst.hpp:56
LRATE_MODE
const USI LRATE_MODE
Well option = fixed fluid rate???
Definition: OCPConst.hpp:134
WRATE_MODE
const USI WRATE_MODE
Well option = fixed water rate.
Definition: OCPConst.hpp:133
ORATE_MODE
const USI ORATE_MODE
Well option = fixed oil rate.
Definition: OCPConst.hpp:131
RATE_MODE
const USI RATE_MODE
Well option = fixed total rate???
Definition: OCPConst.hpp:130
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:25
PROD
const USI PROD
Well type = producer.
Definition: OCPConst.hpp:123
SCALARFASP
const USI SCALARFASP
Use scalar linear solver in Fasp.
Definition: OCPConst.hpp:87
INJ
const USI INJ
Well type = injector.
Definition: OCPConst.hpp:122
BHP_MODE
const USI BHP_MODE
Well option = fixed bottom-hole-pressure.
Definition: OCPConst.hpp:135
OCP_BOOL
unsigned int OCP_BOOL
OCP_BOOL in OCP.
Definition: OCPConst.hpp:29
GRATE_MODE
const USI GRATE_MODE
Well option = fixed gas rate.
Definition: OCPConst.hpp:132
CONV2
const OCP_DBL CONV2
Darcy constant in field unit.
Definition: OCPConst.hpp:64
OCPFluidMethod.hpp
Declaration of solution methods for fluid part in OpenCAEPoro.
OCP_FUNCNAME
#define OCP_FUNCNAME
Print Function Name.
Definition: UtilError.hpp:73
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47
AllWells::CalFlux
void CalFlux(const Bulk &myBulk)
Calculate volume flow rate and moles flow rate of each perforation.
Definition: AllWells.cpp:227
AllWells::wells
vector< Well > wells
well set.
Definition: AllWells.hpp:174
AllWells::CalTrans
void CalTrans(const Bulk &myBulk)
Calculate Transmissibility of Wells.
Definition: AllWells.cpp:216
AllWells::CaldG
void CaldG(const Bulk &myBulk)
Calculate dG.
Definition: AllWells.cpp:238
AllWells::SetupWellBulk
void SetupWellBulk(Bulk &myBulk) const
Setup bulks which are penetrated by wells at current stage.
Definition: AllWells.cpp:156
AllWells::InitBHP
void InitBHP(const Bulk &myBulk)
Set the initial well pressure.
Definition: AllWells.cpp:192
AllWells::PrepareWell
void PrepareWell(const Bulk &myBulk)
Calculate well properties at the beginning of each time step.
Definition: AllWells.cpp:201
AllWells::numWell
USI numWell
num of wells.
Definition: AllWells.hpp:173
BulkConn
Properties and operations on connections between bulks (active grids).
Definition: BulkConn.hpp:71
BulkConn::numConn
OCP_USI numConn
Number of connections between bulks.
Definition: BulkConn.hpp:108
BulkConn::upblock
vector< OCP_USI > upblock
Index of upwinding bulk of connections : numConn * numPhase.
Definition: BulkConn.hpp:133
BulkConn::upblock_Rho
vector< OCP_DBL > upblock_Rho
Mass density of phase from upblock: numConn * numPhase.
Definition: BulkConn.hpp:135
BulkConn::lupblock_Velocity
vector< OCP_DBL > lupblock_Velocity
last upblock_Velocity
Definition: BulkConn.hpp:146
BulkConn::upblock_Trans
vector< OCP_DBL > upblock_Trans
Transmissibility of phase from upblock: numConn * numPhase.
Definition: BulkConn.hpp:137
BulkConn::lupblock_Rho
vector< OCP_DBL > lupblock_Rho
last upblock_Rho
Definition: BulkConn.hpp:144
BulkConn::upblock_Velocity
vector< OCP_DBL > upblock_Velocity
Definition: BulkConn.hpp:138
BulkConn::neighbor
vector< vector< OCP_USI > > neighbor
Neighboring information of each bulk: activeGridNum.
Definition: BulkConn.hpp:113
BulkConn::lupblock_Trans
vector< OCP_DBL > lupblock_Trans
last upblock_Trans
Definition: BulkConn.hpp:145
BulkConn::iteratorConn
vector< BulkPair > iteratorConn
All connections (pair of indices) between bulks: numConn.
Definition: BulkConn.hpp:124
BulkConn::lupblock
vector< OCP_USI > lupblock
last upblock
Definition: BulkConn.hpp:143
Bulk
Physical information of each active reservoir bulk.
Definition: Bulk.hpp:104
Bulk::mu
vector< OCP_DBL > mu
Viscosity of phase: numPhase*numBulk.
Definition: Bulk.hpp:320
Bulk::pVnumCom
vector< USI > pVnumCom
num of variable components in the phase
Definition: Bulk.hpp:488
Bulk::dNNR
vector< OCP_DBL > dNNR
Ni change between NR steps.
Definition: Bulk.hpp:411
Bulk::lP
vector< OCP_DBL > lP
last P
Definition: Bulk.hpp:332
Bulk::lS
vector< OCP_DBL > lS
last S
Definition: Bulk.hpp:336
Bulk::NRstep
vector< OCP_DBL > NRstep
NRstep for FIM.
Definition: Bulk.hpp:422
Bulk::rock
vector< Rock * > rock
rock model
Definition: Bulk.hpp:203
Bulk::Pc
vector< OCP_DBL > Pc
Capillary pressure of phase: numPhase*numBulk.
Definition: Bulk.hpp:312
Bulk::flow
vector< FlowUnit * > flow
Vector for capillary pressure, relative perm.
Definition: Bulk.hpp:197
Bulk::ldPcj_dS
vector< OCP_DBL > ldPcj_dS
last Pcj_dS
Definition: Bulk.hpp:385
Bulk::lNt
vector< OCP_DBL > lNt
last Nt
Definition: Bulk.hpp:328
Bulk::numPhase
USI numPhase
Number of phase.
Definition: Bulk.hpp:154
Bulk::lvf
vector< OCP_DBL > lvf
last vf
Definition: Bulk.hpp:330
Bulk::T
vector< OCP_DBL > T
Temperature: numBulk.
Definition: Bulk.hpp:308
Bulk::S
vector< OCP_DBL > S
Saturation of phase: numPhase*numBulk.
Definition: Bulk.hpp:314
Bulk::rho
vector< OCP_DBL > rho
Mass density of phase: numPhase*numBulk.
Definition: Bulk.hpp:318
Bulk::ifComps
OCP_BOOL ifComps
If OCP_TRUE, compositional model will be used.
Definition: Bulk.hpp:219
Bulk::res
OCPRes res
Residual for all equations.
Definition: Bulk.hpp:477
Bulk::dSec_dPri
vector< OCP_DBL > dSec_dPri
d Secondary variable / d Primary variable.
Definition: Bulk.hpp:486
Bulk::lxij
vector< OCP_DBL > lxij
last xij
Definition: Bulk.hpp:339
Bulk::rhoP
vector< OCP_DBL > rhoP
d Rho / d P: numPhase*numBulk.
Definition: Bulk.hpp:352
Bulk::ntg
vector< OCP_DBL > ntg
net to gross of bulk.
Definition: Bulk.hpp:243
Bulk::lrho
vector< OCP_DBL > lrho
last rho
Definition: Bulk.hpp:340
Bulk::PVTNUM
vector< USI > PVTNUM
Identify PVT region in black-oil model: numBulk.
Definition: Bulk.hpp:191
Bulk::lphaseNum
vector< USI > lphaseNum
last phaseNum
Definition: Bulk.hpp:327
Bulk::vfi
vector< OCP_DBL > vfi
d vf / d Ni: numCom*numBulk.
Definition: Bulk.hpp:351
Bulk::lpVnumCom
vector< USI > lpVnumCom
last pVnumCom
Definition: Bulk.hpp:496
Bulk::poro
vector< OCP_DBL > poro
rock porosity * ntg.
Definition: Bulk.hpp:245
Bulk::xij
vector< OCP_DBL > xij
Nij / Nj: numPhase*numCom*numBulk.
Definition: Bulk.hpp:317
Bulk::lrockVp
vector< OCP_DBL > lrockVp
Pore volume: numBulk.
Definition: Bulk.hpp:256
Bulk::kr
vector< OCP_DBL > kr
Relative permeability of phase: numPhase*numBulk.
Definition: Bulk.hpp:321
Bulk::SATNUM
vector< USI > SATNUM
Identify SAT region: numBulk.
Definition: Bulk.hpp:196
Bulk::depth
vector< OCP_DBL > depth
Depth of center of grid cells: activeGridNum.
Definition: Bulk.hpp:241
Bulk::rockVp
vector< OCP_DBL > rockVp
pore volume = Vgrid * ntg * poro.
Definition: Bulk.hpp:246
Bulk::Ni
vector< OCP_DBL > Ni
Moles of component: numCom*numBulk.
Definition: Bulk.hpp:306
Bulk::NRdPmax
OCP_DBL NRdPmax
Max pressure difference in an NR step.
Definition: Bulk.hpp:417
Bulk::Pj
vector< OCP_DBL > Pj
Pressure of phase: numPhase*numBulk.
Definition: Bulk.hpp:311
Bulk::mux
vector< OCP_DBL > mux
d Muj / d xij: numPhase*numCom*numBulk.
Definition: Bulk.hpp:360
Bulk::rhox
vector< OCP_DBL > rhox
d Rhoj / d xij: numPhase*numCom*numBulk.
Definition: Bulk.hpp:354
Bulk::lporoP
vector< OCP_DBL > lporoP
last poroP.
Definition: Bulk.hpp:268
Bulk::index_maxNRdSSP
OCP_USI index_maxNRdSSP
index of grid which has maxNRdSSP
Definition: Bulk.hpp:416
Bulk::lkr
vector< OCP_DBL > lkr
last kr
Definition: Bulk.hpp:343
Bulk::dSNR
vector< OCP_DBL > dSNR
saturation change between NR steps
Definition: Bulk.hpp:406
Bulk::lvfi
vector< OCP_DBL > lvfi
last vfi
Definition: Bulk.hpp:375
Bulk::lporo
vector< OCP_DBL > lporo
last poro.
Definition: Bulk.hpp:255
Bulk::maxNRdSSP
OCP_DBL maxNRdSSP
max difference between dSNR and dSNRP
Definition: Bulk.hpp:415
Bulk::lvj
vector< OCP_DBL > lvj
last vj
Definition: Bulk.hpp:337
Bulk::Nt
vector< OCP_DBL > Nt
Total moles of components in bulks: numBulk.
Definition: Bulk.hpp:305
Bulk::muP
vector< OCP_DBL > muP
d Mu / d P: numPhase*numBulk.
Definition: Bulk.hpp:358
Bulk::dKr_dS
vector< OCP_DBL > dKr_dS
d Krj / d Sk: numPhase * numPhase * bulk.
Definition: Bulk.hpp:362
Bulk::lT
vector< OCP_DBL > lT
last T
Definition: Bulk.hpp:331
Bulk::flashCal
vector< Mixture * > flashCal
Flash calculation class.
Definition: Bulk.hpp:192
Bulk::phaseExist
vector< OCP_BOOL > phaseExist
Existence of phase: numPhase*numBulk.
Definition: Bulk.hpp:313
Bulk::lPj
vector< OCP_DBL > lPj
last Pj
Definition: Bulk.hpp:333
Bulk::xiP
vector< OCP_DBL > xiP
d xi / d P: numPhase*numBulk.
Definition: Bulk.hpp:355
Bulk::poroP
vector< OCP_DBL > poroP
d poro / d P.
Definition: Bulk.hpp:261
Bulk::satcm
vector< vector< OCP_DBL > > satcm
critical saturation when phase becomes mobile / immobile.
Definition: Bulk.hpp:199
Bulk::ldKr_dS
vector< OCP_DBL > ldKr_dS
last dKr_dS
Definition: Bulk.hpp:386
Bulk::lrhoP
vector< OCP_DBL > lrhoP
last rhoP
Definition: Bulk.hpp:376
Bulk::xi
vector< OCP_DBL > xi
Moles density of phase: numPhase*numBulk.
Definition: Bulk.hpp:319
Bulk::lxi
vector< OCP_DBL > lxi
last xi
Definition: Bulk.hpp:341
Bulk::lmux
vector< OCP_DBL > lmux
last mux
Definition: Bulk.hpp:384
Bulk::vf
vector< OCP_DBL > vf
Total fluid volume: numBulk.
Definition: Bulk.hpp:307
Bulk::phaseNum
vector< USI > phaseNum
Num of hydrocarbon phase in each bulk.
Definition: Bulk.hpp:304
Bulk::bRowSizedSdP
vector< USI > bRowSizedSdP
length of dSec_dPri in each bulk
Definition: Bulk.hpp:485
Bulk::vfP
vector< OCP_DBL > vfP
d vf / d P: numBulk.
Definition: Bulk.hpp:349
Bulk::NRdNmax
OCP_DBL NRdNmax
Max Ni difference in an NR step.
Definition: Bulk.hpp:419
Bulk::lres_n
vector< OCP_DBL > lres_n
last res_n
Definition: Bulk.hpp:497
Bulk::lxix
vector< OCP_DBL > lxix
last xix
Definition: Bulk.hpp:381
Bulk::lresPc
vector< OCP_DBL > lresPc
last lresPc;
Definition: Bulk.hpp:498
Bulk::IfUseEoS
OCP_BOOL IfUseEoS() const
Return ifUseEoS.
Definition: Bulk.hpp:215
Bulk::dPNR
vector< OCP_DBL > dPNR
P change between NR steps.
Definition: Bulk.hpp:412
Bulk::wellBulkId
vector< OCP_USI > wellBulkId
Definition: Bulk.hpp:507
Bulk::vj
vector< OCP_DBL > vj
Volume of phase: numPhase*numBulk.
Definition: Bulk.hpp:315
Bulk::lmuP
vector< OCP_DBL > lmuP
last muP
Definition: Bulk.hpp:382
Bulk::lvfP
vector< OCP_DBL > lvfP
last vfP
Definition: Bulk.hpp:373
Bulk::Pb
vector< OCP_DBL > Pb
Bubble point pressure: numBulk.
Definition: Bulk.hpp:310
Bulk::lrhox
vector< OCP_DBL > lrhox
last rhox
Definition: Bulk.hpp:378
Bulk::lphaseExist
vector< OCP_BOOL > lphaseExist
last phaseExist
Definition: Bulk.hpp:335
Bulk::v
vector< OCP_DBL > v
Volume of grids: activeGridNum.
Definition: Bulk.hpp:240
Bulk::numComH
USI numComH
Number of HydroCarbon.
Definition: Bulk.hpp:156
Bulk::dSNRP
vector< OCP_DBL > dSNRP
predicted saturation change between NR steps
Definition: Bulk.hpp:410
Bulk::InitPTSw
void InitPTSw(const USI &tabrow)
Calculate initial equilibrium – hydrostatic equilibration.
Definition: Bulk.cpp:443
Bulk::lxiP
vector< OCP_DBL > lxiP
last xiP
Definition: Bulk.hpp:379
Bulk::lpSderExist
vector< OCP_BOOL > lpSderExist
last pSderExist
Definition: Bulk.hpp:495
Bulk::resPc
vector< OCP_DBL > resPc
a precalculated value
Definition: Bulk.hpp:490
Bulk::lmu
vector< OCP_DBL > lmu
last mu
Definition: Bulk.hpp:342
Bulk::pSderExist
vector< OCP_BOOL > pSderExist
Existence of derivative of phase saturation.
Definition: Bulk.hpp:487
Bulk::numBulk
OCP_USI numBulk
Number of bulks (active grids).
Definition: Bulk.hpp:153
Bulk::NRphaseNum
vector< USI > NRphaseNum
phaseNum in NR step
Definition: Bulk.hpp:423
Bulk::poroInit
vector< OCP_DBL > poroInit
initial rock porosity * ntg.
Definition: Bulk.hpp:244
Bulk::ldSec_dPri
vector< OCP_DBL > ldSec_dPri
last dSec_dPri
Definition: Bulk.hpp:494
Bulk::xix
vector< OCP_DBL > xix
d Xi_j / d xij: numPhase*numCom*numBulk.
Definition: Bulk.hpp:357
Bulk::ROCKNUM
vector< USI > ROCKNUM
index of Rock table for each bulk
Definition: Bulk.hpp:202
Bulk::nj
vector< OCP_DBL > nj
moles number of phase: numPhase*numBulk.
Definition: Bulk.hpp:316
Bulk::dPcj_dS
vector< OCP_DBL > dPcj_dS
d Pcj / d Sk: numPhase * numPhase * bulk.
Definition: Bulk.hpp:361
Bulk::res_n
vector< OCP_DBL > res_n
residual for FIM_n
Definition: Bulk.hpp:489
Bulk::lbRowSizedSdP
vector< USI > lbRowSizedSdP
last bRowSizedSdP
Definition: Bulk.hpp:493
Bulk::xijNR
vector< OCP_DBL > xijNR
store the current NR step's xij in AIM
Definition: Bulk.hpp:522
Bulk::P
vector< OCP_DBL > P
Pressure: numBulk.
Definition: Bulk.hpp:309
Bulk::lNi
vector< OCP_DBL > lNi
last Ni
Definition: Bulk.hpp:329
Bulk::lnj
vector< OCP_DBL > lnj
last nj
Definition: Bulk.hpp:338
Bulk::cfl
vector< OCP_DBL > cfl
CFL number for each bulk.
Definition: Bulk.hpp:465
Bulk::lPc
vector< OCP_DBL > lPc
last Pc
Definition: Bulk.hpp:334
Bulk::maxLendSdP
USI maxLendSdP
length of dSec_dPri.
Definition: Bulk.hpp:484
Bulk::numCom
USI numCom
Number of component.
Definition: Bulk.hpp:155
ControlNR::NRtol
OCP_DBL NRtol
Maximum non-linear convergence error.
Definition: OCPControl.hpp:68
ControlNR::NRdSmin
OCP_DBL NRdSmin
Minimum Saturation change in a Newton iteration.
Definition: OCPControl.hpp:72
ControlNR::maxNRiter
USI maxNRiter
Maximum number of Newton iterations in a time step.
Definition: OCPControl.hpp:67
ControlNR::NRdSmax
OCP_DBL NRdSmax
Maximum Saturation change in a Newton iteration.
Definition: OCPControl.hpp:70
ControlNR::NRdPmin
OCP_DBL NRdPmin
Minimum Pressure change in a Newton iteration.
Definition: OCPControl.hpp:71
ControlTime::cutFacNR
OCP_DBL cutFacNR
Factor by which time step is cut after convergence failure.
Definition: OCPControl.hpp:39
GetWallTime
Get elapsed wall-time in millisecond.
Definition: UtilTiming.hpp:32
GetWallTime::Stop
__inline__ double Stop() const
Stop the timer and return duration from start() in ms.
Definition: UtilTiming.hpp:54
GetWallTime::Start
__inline__ void Start()
Start the timer.
Definition: UtilTiming.hpp:51
IsoT_AIMc::GetSolution
void GetSolution(Reservoir &rs, const vector< OCP_DBL > &u, const OCPControl &ctrl) const
Update P, Ni, BHP after linear system is solved.
Definition: OCPFluidMethod.cpp:4843
IsoT_AIMc::AssembleMatBulks
void AssembleMatBulks(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble linear system for bulks.
Definition: OCPFluidMethod.cpp:4380
IsoT_AIMc::CalFlashEa
void CalFlashEa(Bulk &bk)
Perform flash calculation with Ni for Explicit bulk – Update all properties.
Definition: OCPFluidMethod.cpp:4276
IsoT_AIMc::SolveLinearSystem
void SolveLinearSystem(LinearSystem &ls, Reservoir &rs, OCPControl &ctrl)
Solve the linear system.
Definition: OCPFluidMethod.cpp:4069
IsoT_AIMc::CalFlashEp
void CalFlashEp(Bulk &bk)
Perform flash calculation with Ni for Explicit bulk – Update partial properties.
Definition: OCPFluidMethod.cpp:4257
IsoT_AIMc::CalKrPcI
void CalKrPcI(Bulk &bk)
Calculate relative permeability and capillary pressure for Implicit bulk.
Definition: OCPFluidMethod.cpp:4363
IsoT_AIMc::Setup
void Setup(Reservoir &rs, LinearSystem &ls, const OCPControl &ctrl)
Setup AIMc.
Definition: OCPFluidMethod.cpp:4020
IsoT_AIMc::AssembleMat
void AssembleMat(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble Matrix.
Definition: OCPFluidMethod.cpp:4060
IsoT_AIMc::FinishNR
OCP_BOOL FinishNR(Reservoir &rs, OCPControl &ctrl)
Finish a Newton-Raphson iteration.
Definition: OCPFluidMethod.cpp:4123
IsoT_AIMc::CalFlashI
void CalFlashI(Bulk &bk)
Perform flash calculation with Ni for Implicit bulk.
Definition: OCPFluidMethod.cpp:4296
IsoT_AIMc::SetFIMBulk
void SetFIMBulk(Reservoir &rs)
Determine which bulk are treated Implicit.
Definition: OCPFluidMethod.cpp:4191
IsoT_AIMc::UpdateProperty
OCP_BOOL UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:4098
IsoT_AIMc::UpdateLastTimeStep
void UpdateLastTimeStep(Reservoir &rs) const
Update values of last step for AIMc.
Definition: OCPFluidMethod.cpp:4983
IsoT_AIMc::PassFlashValueEp
void PassFlashValueEp(Bulk &bk, const OCP_USI &n)
Pass flash value needed for Explicit bulk – Update partial properties.
Definition: OCPFluidMethod.cpp:4317
IsoT_AIMc::CalKrPcE
void CalKrPcE(Bulk &bk)
Calculate relative permeability and capillary pressure for Explicit bulk.
Definition: OCPFluidMethod.cpp:4348
IsoT_AIMc::ResetToLastTimeStep
void ResetToLastTimeStep(Reservoir &rs, OCPControl &ctrl)
Reset variables to last time step.
Definition: OCPFluidMethod.cpp:4976
IsoT_AIMc::Prepare
void Prepare(Reservoir &rs, const OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:4043
IsoT_AIMc::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:4028
IsoT_AIMc::FinishStep
void FinishStep(Reservoir &rs, OCPControl &ctrl) const
Finish a time step.
Definition: OCPFluidMethod.cpp:4165
IsoT_AIMc::AllocateReservoir
void AllocateReservoir(Reservoir &rs)
Allocate memory for reservoir.
Definition: OCPFluidMethod.cpp:4174
IsoT_FIM::Setup
void Setup(Reservoir &rs, LinearSystem &ls, const OCPControl &ctrl)
Setup FIM.
Definition: OCPFluidMethod.cpp:828
IsoT_FIM::AssembleMat
void AssembleMat(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble Matrix.
Definition: OCPFluidMethod.cpp:860
IsoT_FIM::UpdateLastTimeStep
void UpdateLastTimeStep(Reservoir &rs) const
Update values of last step for FIM.
Definition: OCPFluidMethod.cpp:2915
IsoT_FIM::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:836
IsoT_FIM::FinishNR
OCP_BOOL FinishNR(Reservoir &rs, OCPControl &ctrl)
Finish a Newton-Raphson iteration.
Definition: OCPFluidMethod.cpp:939
IsoT_FIM::AllocateLinearSystem
void AllocateLinearSystem(LinearSystem &ls, const Reservoir &rs, const OCPControl &ctrl)
Allocate memory for linear system.
Definition: OCPFluidMethod.cpp:1086
IsoT_FIM::Prepare
void Prepare(Reservoir &rs, const OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:852
IsoT_FIM::FinishStep
void FinishStep(Reservoir &rs, OCPControl &ctrl)
Finish a time step.
Definition: OCPFluidMethod.cpp:972
IsoT_FIM::PassFlashValue
void PassFlashValue(Bulk &bk, const OCP_USI &n) const
Pass value needed for FIM from flash to bulk.
Definition: OCPFluidMethod.cpp:1122
IsoT_FIM::SolveLinearSystem
void SolveLinearSystem(LinearSystem &ls, Reservoir &rs, OCPControl &ctrl) const
Solve the linear system.
Definition: OCPFluidMethod.cpp:877
IsoT_FIM::UpdateProperty
OCP_BOOL UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:914
IsoT_FIM::CalRes
void CalRes(Reservoir &rs, const OCP_DBL &dt, const OCP_BOOL &resetRes0) const
Calculate residual.
Definition: OCPFluidMethod.cpp:1201
IsoT_FIM::CalKrPc
void CalKrPc(Bulk &bk) const
Calculate relative permeability and capillary pressure needed for FIM.
Definition: OCPFluidMethod.cpp:1189
IsoT_FIM::ResetToLastTimeStep
void ResetToLastTimeStep(Reservoir &rs, OCPControl &ctrl)
Reset variables to last time step.
Definition: OCPFluidMethod.cpp:2859
IsoT_FIM::AllocateReservoir
void AllocateReservoir(Reservoir &rs)
Allocate memory for reservoir.
Definition: OCPFluidMethod.cpp:980
IsoT_FIM::AssembleMatWells
void AssembleMatWells(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble linear system for wells.
Definition: OCPFluidMethod.cpp:2135
IsoT_FIMn::AssembleMatBulksNew
void AssembleMatBulksNew(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble linear system for bulks.
Definition: OCPFluidMethod.cpp:3198
IsoT_FIMn::FinishStep
void FinishStep(Reservoir &rs, OCPControl &ctrl) const
Finish a time step.
Definition: OCPFluidMethod.cpp:3103
IsoT_FIMn::Prepare
void Prepare(Reservoir &rs, const OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:2991
IsoT_FIMn::AssembleMatProdWellsNew
void AssembleMatProdWellsNew(LinearSystem &ls, const Bulk &bk, const Well &wl, const OCP_DBL &dt) const
Assemble linear system for production wells.
Definition: OCPFluidMethod.cpp:3679
IsoT_FIMn::GetSolution
void GetSolution(Reservoir &rs, const vector< OCP_DBL > &u, const OCPControl &ctrl) const
Update P, Ni, BHP after linear system is solved.
Definition: OCPFluidMethod.cpp:3847
IsoT_FIMn::AllocateReservoir
void AllocateReservoir(Reservoir &rs)
Allocate memory for reservoir.
Definition: OCPFluidMethod.cpp:3111
IsoT_FIMn::ResetToLastTimeStep
void ResetToLastTimeStep(Reservoir &rs, OCPControl &ctrl)
Reset variables to last time step.
Definition: OCPFluidMethod.cpp:3999
IsoT_FIMn::UpdateProperty
OCP_BOOL UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:3042
IsoT_FIMn::AssembleMatInjWellsNew
void AssembleMatInjWellsNew(LinearSystem &ls, const Bulk &bk, const Well &wl, const OCP_DBL &dt) const
Assemble linear system for injection wells.
Definition: OCPFluidMethod.cpp:3524
IsoT_FIMn::AssembleMatWellsNew
void AssembleMatWellsNew(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble linear system for wells.
Definition: OCPFluidMethod.cpp:3503
IsoT_FIMn::InitFlash
void InitFlash(Bulk &bk) const
Perform Flash with Sj and calculate values needed for FIMn.
Definition: OCPFluidMethod.cpp:3129
IsoT_FIMn::PassFlashValue
void PassFlashValue(Bulk &bk, const OCP_USI &n) const
Pass value needed for FIMn from flash to bulk.
Definition: OCPFluidMethod.cpp:3179
IsoT_FIMn::AssembleMat
void AssembleMat(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble Matrix.
Definition: OCPFluidMethod.cpp:2998
IsoT_FIMn::Setup
void Setup(Reservoir &rs, LinearSystem &ls, const OCPControl &ctrl)
Setup FIM.
Definition: OCPFluidMethod.cpp:2967
IsoT_FIMn::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:2975
IsoT_FIMn::CalFlash
void CalFlash(Bulk &bk)
Perform Flash with Ni and calculate values needed for FIMn.
Definition: OCPFluidMethod.cpp:3152
IsoT_FIMn::SolveLinearSystem
void SolveLinearSystem(LinearSystem &ls, Reservoir &rs, OCPControl &ctrl) const
Solve the linear system.
Definition: OCPFluidMethod.cpp:3008
IsoT_FIMn::FinishNR
OCP_BOOL FinishNR(Reservoir &rs, OCPControl &ctrl)
Finish a Newton-Raphson iteration.
Definition: OCPFluidMethod.cpp:3068
IsoT_FIMn::UpdateLastTimeStep
void UpdateLastTimeStep(Reservoir &rs) const
Update values of last step for FIMn.
Definition: OCPFluidMethod.cpp:4008
IsoT_IMPEC::Prepare
void Prepare(Reservoir &rs, OCPControl &ctrl)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:64
IsoT_IMPEC::SolveLinearSystem
void SolveLinearSystem(LinearSystem &ls, Reservoir &rs, OCPControl &ctrl)
Solve the linear system.
Definition: OCPFluidMethod.cpp:79
IsoT_IMPEC::AssembleMat
void AssembleMat(LinearSystem &ls, const Reservoir &rs, const OCP_DBL &dt) const
Assemble Matrix.
Definition: OCPFluidMethod.cpp:71
IsoT_IMPEC::InitFlash
void InitFlash(Bulk &bk) const
Perform Flash with Sj and calculate values needed for FIM.
Definition: OCPFluidMethod.cpp:251
IsoT_IMPEC::Setup
void Setup(Reservoir &rs, LinearSystem &ls, const OCPControl &ctrl)
Setup IMPEC.
Definition: OCPFluidMethod.cpp:38
IsoT_IMPEC::UpdateProperty
OCP_BOOL UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:111
IsoT_IMPEC::PassFlashValue
void PassFlashValue(Bulk &bk, const OCP_USI &n) const
Pass value needed for FIM from flash to bulk.
Definition: OCPFluidMethod.cpp:275
IsoT_IMPEC::CalKrPc
void CalKrPc(Bulk &bk) const
Calculate relative permeability and capillary pressure needed for FIM.
Definition: OCPFluidMethod.cpp:311
IsoT_IMPEC::FinishNR
OCP_BOOL FinishNR(const Reservoir &rs)
Determine if NR iteration finishes.
Definition: OCPFluidMethod.cpp:154
IsoT_IMPEC::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:47
LinearSystem
Linear solvers for discrete systems.
Definition: LinearSystem.hpp:31
LinearSystem::GetNumIters
USI GetNumIters()
Return the number of iterations.
Definition: LinearSystem.hpp:127
LinearSystem::AssembleMatLinearSolver
void AssembleMatLinearSolver()
Assemble Mat for Linear Solver.
Definition: LinearSystem.hpp:57
LinearSystem::AllocateColMem
void AllocateColMem(const vector< USI > &bulk2bulk, const vector< vector< OCP_USI >> well2bulk)
Allocate memory for linear system with max possible number of columns.
Definition: LinearSystem.cpp:26
LinearSystem::AddDiag
void AddDiag(const OCP_USI &n, const OCP_DBL &v)
Add a value at diagonal value.
Definition: LinearSystem.hpp:77
LinearSystem::AssembleRhsCopy
void AssembleRhsCopy(const vector< OCP_DBL > &rhs)
Assign Rhs by Copying.
Definition: LinearSystem.cpp:78
LinearSystem::AssembleRhsAccumulate
void AssembleRhsAccumulate(const vector< OCP_DBL > &rhs)
Assign Rhs by Accumulating.
Definition: LinearSystem.cpp:70
LinearSystem::AllocateRowMem
void AllocateRowMem(const OCP_USI &dimMax, const USI &nb)
Allocate memory for linear system with max possible number of rows.
Definition: LinearSystem.cpp:13
LinearSystem::NewDiag
void NewDiag(const OCP_USI &n, const OCP_DBL &v)
Push back a diagonal val, which is always at the first location.
Definition: LinearSystem.hpp:70
LinearSystem::AssignGuess
void AssignGuess(const OCP_USI &n, const OCP_DBL &v)
Assign an initial value at u[n].
Definition: LinearSystem.hpp:92
LinearSystem::AddRhs
void AddRhs(const OCP_USI &n, const OCP_DBL &v)
Add a value at b[n].
Definition: LinearSystem.hpp:90
LinearSystem::ClearData
void ClearData()
Clear the internal matrix data for scalar-value problems.
Definition: LinearSystem.cpp:57
LinearSystem::NewOffDiag
void NewOffDiag(const OCP_USI &bId, const OCP_USI &eId, const OCP_DBL &v)
Push back a off-diagonal value.
Definition: LinearSystem.hpp:83
LinearSystem::Solve
OCP_INT Solve()
Solve the Linear System.
Definition: LinearSystem.hpp:59
LinearSystem::AddDim
OCP_USI AddDim(const OCP_USI &n)
Setup dimensions.
Definition: LinearSystem.hpp:62
LinearSystem::CheckEquation
void CheckEquation() const
Check whether NAN or INF occurs in equations, used in debug mode.
Definition: LinearSystem.cpp:141
LinearSystem::GetSolution
vector< OCP_DBL > & GetSolution()
Return the solution.
Definition: LinearSystem.hpp:43
LinearSystem::SetupLinearSolver
void SetupLinearSolver(const USI &i, const string &dir, const string &file)
Setup LinearSolver.
Definition: LinearSystem.cpp:171
LinearSystem::CheckSolution
void CheckSolution() const
Check whether NAN or INF occurs in solutions, used in debug mode.
Definition: LinearSystem.cpp:160
OCPControl
All control parameters except for well controllers.
Definition: OCPControl.hpp:94
OCPControl::GetWorkDir
string GetWorkDir() const
Return work dir name.
Definition: OCPControl.hpp:181
OCPControl::GetCurDt
OCP_DBL GetCurDt() const
Return current time step size.
Definition: OCPControl.hpp:133
OCPControl::UpdateIterNR
void UpdateIterNR()
Update the number of Newton iterations.
Definition: OCPControl.hpp:157
OCPControl::ResetIterNRLS
void ResetIterNRLS()
Reset the number of iterations.
Definition: OCPControl.cpp:234
OCPControl::GetLsFile
string GetLsFile() const
Return linear solver file name.
Definition: OCPControl.hpp:184
OCPControl::UpdateIterLS
void UpdateIterLS(const USI &num)
Update the number of linear iterations.
Definition: OCPControl.hpp:154
OCPControl::RecordTimeLS
void RecordTimeLS(const OCP_DBL &t)
Record time used for linear solver.
Definition: OCPControl.hpp:163
OCPRes::maxWellRelRes_mol
OCP_DBL maxWellRelRes_mol
maximum relative residual wrt. total moles for each well
Definition: OCPStructure.hpp:71
OCPRes::maxId_N
OCP_INT maxId_N
index of bulk who has maxRelRes_N
Definition: OCPStructure.hpp:75
OCPRes::maxRelRes_V
OCP_DBL maxRelRes_V
Definition: OCPStructure.hpp:66
OCPRes::maxId_V
OCP_INT maxId_V
index of bulk who has maxRelRes_V
Definition: OCPStructure.hpp:74
OCPRes::resRelN
vector< OCP_DBL > resRelN
Definition: OCPStructure.hpp:60
OCPRes::maxRelRes0_V
OCP_DBL maxRelRes0_V
Definition: OCPStructure.hpp:64
OCPRes::resRelV
vector< OCP_DBL > resRelV
Definition: OCPStructure.hpp:58
OCPRes::resAbs
vector< OCP_DBL > resAbs
residual for all equations for each bulk
Definition: OCPStructure.hpp:57
OCPRes::maxRelRes_N
OCP_DBL maxRelRes_N
maximum relative residual wrt. total moles for each bulk
Definition: OCPStructure.hpp:68
Reservoir::conn
BulkConn conn
Bulk's connection info.
Definition: Reservoir.hpp:76
Reservoir::GetNRdSmax
OCP_DBL GetNRdSmax(OCP_USI &index)
Return NRdSmax.
Definition: Reservoir.hpp:85
Reservoir::optFeatures
OptionalFeatures optFeatures
optional features.
Definition: Reservoir.hpp:77
Reservoir::CalCFL
OCP_DBL CalCFL(const OCP_DBL &dt) const
Calculate the CFL number, including bulks and wells for IMPEC.
Definition: Reservoir.cpp:73
Reservoir::GetWellNum
USI GetWellNum() const
Return the num of Well.
Definition: Reservoir.hpp:68
Reservoir::GetNRdPmax
OCP_DBL GetNRdPmax()
Return NRdPmax.
Definition: Reservoir.hpp:83
Reservoir::CalMaxChange
void CalMaxChange()
Calculate Maximum Change of some reference variables for IMPEC.
Definition: Reservoir.cpp:56
Reservoir::GetComNum
USI GetComNum() const
Return the num of Components.
Definition: Reservoir.hpp:70
Reservoir::bulk
Bulk bulk
Active grid info.
Definition: Reservoir.hpp:74
Reservoir::allWells
AllWells allWells
Wells class info.
Definition: Reservoir.hpp:75
Reservoir::CalIPRT
void CalIPRT(const OCP_DBL &dt)
Calculate num of Injection, Production.
Definition: Reservoir.cpp:64
Reservoir::GetBulkNum
OCP_USI GetBulkNum() const
Return the num of Bulk.
Definition: Reservoir.hpp:66
Well
Definition: Well.hpp:45
Well::CalPerfP
void CalPerfP()
Update pressure in Perforation after well pressure updates.
Definition: Well.hpp:134
Well::IsOpen
OCP_BOOL IsOpen() const
Return the state of the well, Open or Close.
Definition: Well.hpp:149
Well::WellType
USI WellType() const
Return the type of well, Inj or Prod.
Definition: Well.hpp:151
Well::CalProdWeight
void CalProdWeight(const Bulk &myBulk) const
Calculate the production weight.
Definition: Well.cpp:829
Generated on Tue Feb 7 2023 01:49:51 for OpenCAEPoro by doxygen 1.9.1