OpenCAEPoro  v0.5.0
A simulator for multicomponent porous media flow
ThermalMethod.cpp
Go to the documentation of this file.
1 
12 #include "ThermalMethod.hpp"
13 
14 void T_FIM::Setup(Reservoir& rs, LinearSystem& ls, const OCPControl& ctrl)
15 {
16  AllocateReservoir(rs);
17  AllocateLinearSystem(ls, rs, ctrl);
18 }
19 
20 void T_FIM::InitReservoir(Reservoir& rs) const
21 {
22  rs.bulk.InitPTSw(50);
23 
24  InitRock(rs.bulk);
25  CalRock(rs.bulk);
26 
27  InitFlash(rs.bulk);
28  CalKrPc(rs.bulk);
29 
30  CalThermalConduct(rs.conn, rs.bulk);
31 
32  rs.allWells.InitBHP(rs.bulk);
33  UpdateLastTimeStep(rs);
34 }
35 
36 void T_FIM::Prepare(Reservoir& rs, const OCPControl& ctrl)
37 {
38  rs.allWells.PrepareWell(rs.bulk);
39  CalRes(rs, ctrl.GetCurTime() + ctrl.GetCurDt(), ctrl.GetCurDt(), OCP_TRUE);
40 }
41 
42 void T_FIM::AssembleMat(LinearSystem& ls,
43  const Reservoir& rs,
44  const OCP_DBL& t,
45  const OCP_DBL& dt) const
46 {
47  AssembleMatBulks(ls, rs, t, dt);
48  AssembleMatWells(ls, rs, dt);
49  ls.AssembleRhsCopy(rs.bulk.res.resAbs);
50 }
51 
52 void T_FIM::SolveLinearSystem(LinearSystem& ls, Reservoir& rs, OCPControl& ctrl) const
53 {
54 #ifdef DEBUG
55  ls.CheckEquation();
56 #endif // DEBUG
57 
58  // Assemble external linear solver with internal A and b
59  ls.AssembleMatLinearSolver();
60  // Solve linear system
61  GetWallTime Timer;
62  Timer.Start();
63  int status = ls.Solve();
64  if (status < 0) {
65  status = ls.GetNumIters();
66  }
67  // Record time, iterations
68  ctrl.RecordTimeLS(Timer.Stop() / 1000);
69  ctrl.UpdateIterLS(status);
70  ctrl.UpdateIterNR();
71 
72 #ifdef DEBUG
73  // Output A, b, x
74  // ls.OutputLinearSystem("testA_FIMT.out", "testb_FIMT.out");
75  // ls.OutputSolution("testx_FIMT.out");
76  // Check if inf or nan occurs in solution
77  ls.CheckSolution();
78 #endif // DEBUG
79 
80  GetSolution(rs, ls.GetSolution(), ctrl);
81  // rs.GetSolution01FIM(ls.GetSolution());
82  // rs.PrintSolFIM(ctrl.workDir + "testPNi.out");
83  ls.ClearData();
84 }
85 
86 OCP_BOOL T_FIM::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
87 {
88  OCP_DBL& dt = ctrl.current_dt;
89 
90  if (!ctrl.Check(rs, {"BulkNi", "BulkP", "BulkT"})) {
91  ResetToLastTimeStep(rs, ctrl);
92  cout << "Cut time step size and repeat! current dt = " << fixed
93  << setprecision(3) << dt << " days\n";
94  return OCP_FALSE;
95  }
96 
97  // Update reservoir properties
98  CalRock(rs.bulk);
99 
100  CalFlash(rs.bulk);
101  CalKrPc(rs.bulk);
102 
103  CalThermalConduct(rs.conn, rs.bulk);
104  CalHeatLoss(rs.bulk, ctrl.GetCurTime() + dt, dt);
105 
106  rs.allWells.CalTrans(rs.bulk);
107  rs.allWells.CalFlux(rs.bulk);
108 
109  CalRes(rs, ctrl.GetCurTime() + dt, dt, OCP_FALSE);
110 
111  return OCP_TRUE;
112 }
113 
114 OCP_BOOL T_FIM::FinishNR(Reservoir& rs, OCPControl& ctrl)
115 {
116  OCP_USI dSn;
117 
118  const OCP_DBL NRdSmax = rs.GetNRdSmax(dSn);
119  const OCP_DBL NRdPmax = rs.GetNRdPmax();
120  // const OCP_DBL NRdNmax = rs.GetNRdNmax();
121 
122  cout << "### NR : " + to_string(ctrl.iterNR) + " Res: " << setprecision(2)
123  << scientific << rs.bulk.res.maxRelRes0_V << setw(12)
124  << rs.bulk.res.maxRelRes_V << setw(12) << rs.bulk.res.maxRelRes_N << setw(12)
125  << rs.bulk.res.maxWellRelRes_mol << setw(12) << rs.bulk.res.maxRelRes_E
126  << setw(12) << fabs(NRdPmax) << setw(12) << fabs(NRdSmax) << endl;
127 
128  if (((rs.bulk.res.maxRelRes_V <= rs.bulk.res.maxRelRes0_V * ctrl.ctrlNR.NRtol ||
129  rs.bulk.res.maxRelRes_V <= ctrl.ctrlNR.NRtol ||
130  rs.bulk.res.maxRelRes_N <= ctrl.ctrlNR.NRtol) &&
131  rs.bulk.res.maxWellRelRes_mol <= ctrl.ctrlNR.NRtol &&
132  rs.bulk.res.maxRelRes_E <= ctrl.ctrlNR.NRtol) ||
133  (fabs(NRdPmax) <= ctrl.ctrlNR.NRdPmin &&
134  fabs(NRdSmax) <= ctrl.ctrlNR.NRdSmin)) {
135 
136  if (!ctrl.Check(rs, {"WellP"})) {
137  ResetToLastTimeStep(rs, ctrl);
138  return OCP_FALSE;
139  } else {
140  return OCP_TRUE;
141  }
142 
143  } else if (ctrl.iterNR >= ctrl.ctrlNR.maxNRiter) {
144  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
145  ResetToLastTimeStep(rs, ctrl);
146  cout << "### WARNING: NR not fully converged! Cut time step size and repeat! "
147  "current dt = "
148  << fixed << setprecision(3) << ctrl.current_dt << " days\n";
149  return OCP_FALSE;
150  } else {
151  return OCP_FALSE;
152  }
153 }
154 
155 void T_FIM::FinishStep(Reservoir& rs, OCPControl& ctrl)
156 {
157  rs.CalIPRT(ctrl.GetCurDt());
158  rs.CalMaxChange();
159  UpdateLastTimeStep(rs);
160  ctrl.CalNextTimeStep(rs, {"dP", "dS", "iter"});
161  ctrl.UpdateIters();
162 }
163 
164 void T_FIM::AllocateReservoir(Reservoir& rs)
165 {
166  Bulk& bk = rs.bulk;
167  const OCP_USI nb = bk.numBulk;
168  const USI np = bk.numPhase;
169  const USI nc = bk.numCom;
170 
171  // Rock
172  bk.poro.resize(nb);
173  bk.rockVp.resize(nb);
174  bk.vr.resize(nb);
175  bk.Hr.resize(nb);
176 
177  bk.lporo.resize(nb);
178  bk.lrockVp.resize(nb);
179  bk.lvr.resize(nb);
180  bk.lHr.resize(nb);
181 
182  // derivatives
183  bk.poroP.resize(nb);
184  bk.poroT.resize(nb);
185  bk.vrP.resize(nb);
186  bk.vrT.resize(nb);
187  bk.HrT.resize(nb);
188 
189  bk.lporoP.resize(nb);
190  bk.lporoT.resize(nb);
191  bk.lvrP.resize(nb);
192  bk.lvrT.resize(nb);
193  bk.lHrT.resize(nb);
194 
195  // Fluid
196  bk.phaseNum.resize(nb);
197  bk.Nt.resize(nb);
198  bk.Ni.resize(nb * nc);
199  bk.vf.resize(nb);
200  bk.T.resize(nb);
201  bk.P.resize(nb);
202  bk.Pb.resize(nb);
203  bk.Pj.resize(nb * np);
204  bk.Pc.resize(nb * np);
205  bk.phaseExist.resize(nb * np);
206  bk.S.resize(nb * np);
207  bk.xij.resize(nb * np * nc);
208  bk.rho.resize(nb * np);
209  bk.xi.resize(nb * np);
210  bk.mu.resize(nb * np);
211  bk.kr.resize(nb * np);
212  bk.Uf.resize(nb);
213  bk.H.resize(nb * np);
214  bk.kt.resize(nb);
215 
216  bk.lphaseNum.resize(nb);
217  bk.lNt.resize(nb);
218  bk.lNi.resize(nb * nc);
219  bk.lvf.resize(nb);
220  bk.lT.resize(nb);
221  bk.lP.resize(nb);
222  bk.lPj.resize(nb * np);
223  bk.lPc.resize(nb * np);
224  bk.lphaseExist.resize(nb * np);
225  bk.lS.resize(nb * np);
226  bk.lxij.resize(nb * np * nc);
227  bk.lrho.resize(nb * np);
228  bk.lxi.resize(nb * np);
229  bk.lmu.resize(nb * np);
230  bk.lkr.resize(nb * np);
231  bk.lUf.resize(nb);
232  bk.lH.resize(nb * np);
233  bk.lkt.resize(nb);
234 
235  // derivatives
236  bk.vfP.resize(nb);
237  bk.vfT.resize(nb);
238  bk.vfi.resize(nb * nc);
239  bk.rhoP.resize(nb * np);
240  bk.rhoT.resize(nb * np);
241  bk.rhox.resize(nb * nc * np);
242  bk.xiP.resize(nb * np);
243  bk.xiT.resize(nb * np);
244  bk.xix.resize(nb * nc * np);
245  bk.muP.resize(nb * np);
246  bk.muT.resize(nb * np);
247  bk.mux.resize(nb * nc * np);
248  bk.dPcj_dS.resize(nb * np * np);
249  bk.dKr_dS.resize(nb * np * np);
250  bk.UfP.resize(nb);
251  bk.UfT.resize(nb);
252  bk.Ufi.resize(nb * nc);
253  bk.HT.resize(nb * np);
254  bk.Hx.resize(nb * np * nc);
255  bk.ktP.resize(nb);
256  bk.ktT.resize(nb);
257  bk.ktS.resize(nb * np);
258 
259  bk.lvfP.resize(nb);
260  bk.lvfT.resize(nb);
261  bk.lvfi.resize(nb * nc);
262  bk.lrhoP.resize(nb * np);
263  bk.lrhoT.resize(nb * np);
264  bk.lrhox.resize(nb * nc * np);
265  bk.lxiP.resize(nb * np);
266  bk.lxiT.resize(nb * np);
267  bk.lxix.resize(nb * nc * np);
268  bk.lmuP.resize(nb * np);
269  bk.lmuT.resize(nb * np);
270  bk.lmux.resize(nb * nc * np);
271  bk.ldPcj_dS.resize(nb * np * np);
272  bk.ldKr_dS.resize(nb * np * np);
273  bk.UfP.resize(nb);
274  bk.UfT.resize(nb);
275  bk.Ufi.resize(nb * nc);
276  bk.HT.resize(nb * np);
277  bk.Hx.resize(nb * np * nc);
278  bk.ktP.resize(nb);
279  bk.ktT.resize(nb);
280  bk.ktS.resize(nb * np);
281 
282  // FIM-Specified
283  bk.maxLendSdP = (nc + 2) * (nc + 1) * np;
284  bk.dSec_dPri.resize(nb * bk.maxLendSdP);
285  bk.bRowSizedSdP.resize(nb);
286  bk.pSderExist.resize(nb * np);
287  bk.pVnumCom.resize(nb * np);
288 
289  bk.ldSec_dPri.resize(nb * bk.maxLendSdP);
290  bk.lbRowSizedSdP.resize(nb);
291  bk.lpSderExist.resize(nb * np);
292  bk.lpVnumCom.resize(nb * np);
293 
294  // Allocate Residual
295  bk.res.SetupT(nb, rs.allWells.numWell, nc);
296 
297  // NR
298  bk.NRstep.resize(nb);
299  bk.NRphaseNum.resize(nb);
300  bk.dSNR.resize(nb * np);
301  bk.dSNRP.resize(nb * np);
302  bk.dNNR.resize(nb * nc);
303  bk.dPNR.resize(nb);
304  bk.dTNR.resize(nb);
305 
306  // BulkConn
307  BulkConn& conn = rs.conn;
308  const OCP_USI& numConn = conn.numConn;
309 
310  conn.upblock.resize(numConn * np);
311  conn.upblock_Rho.resize(numConn * np);
312  conn.upblock_Velocity.resize(numConn * np);
313  conn.Adkt.resize(numConn);
314  conn.AdktP.resize(numConn * 2);
315  conn.AdktT.resize(numConn * 2);
316  conn.AdktS.resize(numConn * 2 * np);
317 }
318 
319 void T_FIM::AllocateLinearSystem(LinearSystem& ls,
320  const Reservoir& rs,
321  const OCPControl& ctrl)
322 {
323  ls.AllocateRowMem(rs.GetBulkNum() + rs.GetWellNum(), rs.GetComNum() + 2);
324  ls.AllocateColMem(rs.conn.GetNeighborNum(), rs.allWells.GetWell2Bulk());
325  ls.SetupLinearSolver(VECTORFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
326 }
327 
328 void T_FIM::InitRock(Bulk& bk) const
329 {
330 
331  for (OCP_USI n = 0; n < bk.numBulk; n++) {
332  if (bk.bType[n] == 0) {
333  // non fluid bulk
334  bk.poroInit[n] = 0;
335  bk.poro[n] = 0;
336  bk.rockVp[n] = 0;
337  bk.vr[n] = bk.v[n];
338  } else {
339  // fluid bulk
340  bk.poroInit[n] *= bk.ntg[n];
341  }
342  }
343 }
344 
345 void T_FIM::CalRock(Bulk& bk) const
346 {
347  const OCP_USI nb = bk.numBulk;
348  OCP_DBL vr, vrP, vrT;
349 
350  for (OCP_USI n = 0; n < nb; n++) {
351  if (bk.bType[n] > 0) {
352  // fluid exists
353  bk.rock[bk.ROCKNUM[n]]->CalPoroT(bk.P[n], bk.T[n], bk.poroInit[n],
354  bk.poro[n], bk.poroP[n], bk.poroT[n], vr,
355  vrP, vrT);
356 
357  bk.rockVp[n] = bk.v[n] * bk.poro[n];
358  bk.vr[n] = bk.v[n] * vr;
359  bk.vrP[n] = bk.v[n] * vrP;
360  bk.vrT[n] = bk.v[n] * vrT;
361  }
362  bk.rock[bk.ROCKNUM[n]]->CalRockHT(bk.T[n], bk.Hr[n], bk.HrT[n]);
363  }
364 }
365 
366 void T_FIM::InitFlash(Bulk& bk) const
367 {
368  const OCP_USI nb = bk.numBulk;
369  const OCP_USI np = bk.numPhase;
370  const OCP_USI nc = bk.numCom;
371 
372  for (OCP_USI n = 0; n < nb; n++) {
373  if (bk.bType[n] > 0) {
374  bk.flashCal[bk.PVTNUM[n]]->InitFlashFIM(bk.P[n], bk.Pb[n], bk.T[n],
375  &bk.S[n * np], bk.rockVp[n],
376  &bk.Ni[n * nc], n);
377  for (USI i = 0; i < nc; i++) {
378  bk.Ni[n * nc + i] = bk.flashCal[bk.PVTNUM[n]]->GetNi(i);
379  }
380  PassFlashValue(bk, n);
381  }
382  }
383 }
384 
385 void T_FIM::CalFlash(Bulk& bk) const
386 {
387  const OCP_USI nb = bk.numBulk;
388  const OCP_USI np = bk.numPhase;
389  const OCP_USI nc = bk.numCom;
390 
391  for (OCP_USI n = 0; n < nb; n++) {
392  if (bk.bType[n] > 0) {
393  bk.flashCal[bk.PVTNUM[n]]->FlashFIM(bk.P[n], bk.T[n], &bk.Ni[n * nc],
394  &bk.S[n * np], bk.phaseNum[n],
395  &bk.xij[n * np * nc], n);
396  PassFlashValue(bk, n);
397  }
398  }
399 }
400 
401 void T_FIM::PassFlashValue(Bulk& bk, const OCP_USI& n) const
402 {
403  const USI np = bk.numPhase;
404  const USI nc = bk.numCom;
405  const OCP_USI bIdp = n * np;
406  const USI pvtnum = bk.PVTNUM[n];
407 
408  bk.phaseNum[n] = 0;
409  bk.Nt[n] = bk.flashCal[pvtnum]->GetNt();
410  bk.vf[n] = bk.flashCal[pvtnum]->GetVf();
411  bk.Uf[n] = bk.flashCal[pvtnum]->GetUf();
412 
413  for (USI j = 0; j < np; j++) {
414  // Important! Saturation must be passed no matter if the phase exists. This is
415  // because it will be used to calculate relative permeability and capillary
416  // pressure at each time step. Make sure that all saturations are updated at
417  // each step!
418  bk.S[bIdp + j] = bk.flashCal[pvtnum]->GetS(j);
419  bk.dSNR[bIdp + j] = bk.S[bIdp + j] - bk.dSNR[bIdp + j];
420  if (bk.phaseExist[bIdp + j]) {
421  if (fabs(bk.maxNRdSSP) < fabs(bk.dSNR[bIdp + j] - bk.dSNRP[bIdp + j])) {
422  bk.maxNRdSSP = bk.dSNR[bIdp + j] - bk.dSNRP[bIdp + j];
423  bk.index_maxNRdSSP = n;
424  }
425  }
426  bk.phaseExist[bIdp + j] = bk.flashCal[pvtnum]->GetPhaseExist(j);
427  if (bk.phaseExist[bIdp + j]) {
428  bk.phaseNum[n]++;
429  bk.rho[bIdp + j] = bk.flashCal[pvtnum]->GetRho(j);
430  bk.xi[bIdp + j] = bk.flashCal[pvtnum]->GetXi(j);
431  bk.mu[bIdp + j] = bk.flashCal[pvtnum]->GetMu(j);
432  bk.H[bIdp + j] = bk.flashCal[pvtnum]->GetH(j);
433 
434  // Derivatives
435  bk.rhoP[bIdp + j] = bk.flashCal[pvtnum]->GetRhoP(j);
436  bk.rhoT[bIdp + j] = bk.flashCal[pvtnum]->GetRhoT(j);
437  bk.xiP[bIdp + j] = bk.flashCal[pvtnum]->GetXiP(j);
438  bk.xiT[bIdp + j] = bk.flashCal[pvtnum]->GetXiT(j);
439  bk.muP[bIdp + j] = bk.flashCal[pvtnum]->GetMuP(j);
440  bk.muT[bIdp + j] = bk.flashCal[pvtnum]->GetMuT(j);
441  bk.HT[bIdp + j] = bk.flashCal[pvtnum]->GetHT(j);
442 
443  for (USI i = 0; i < nc; i++) {
444  bk.xij[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetXij(j, i);
445  bk.rhox[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetRhoX(j, i);
446  bk.xix[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetXiX(j, i);
447  bk.mux[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetMuX(j, i);
448  bk.Hx[bIdp * nc + j * nc + i] = bk.flashCal[pvtnum]->GetHx(j, i);
449  }
450  }
451  }
452  bk.vfP[n] = bk.flashCal[pvtnum]->GetVfP();
453  bk.vfT[n] = bk.flashCal[pvtnum]->GetVfT();
454  bk.UfP[n] = bk.flashCal[pvtnum]->GetUfP();
455  bk.UfT[n] = bk.flashCal[pvtnum]->GetUfT();
456 
457  for (USI i = 0; i < nc; i++) {
458  bk.vfi[n * nc + i] = bk.flashCal[pvtnum]->GetVfi(i);
459  bk.Ufi[n * nc + i] = bk.flashCal[pvtnum]->GetUfi(i);
460  }
461 
462  Dcopy(bk.maxLendSdP, &bk.dSec_dPri[n * bk.maxLendSdP],
463  &bk.flashCal[pvtnum]->GetDXsDXp()[0]);
464 }
465 
466 void T_FIM::CalKrPc(Bulk& bk) const
467 {
468  const USI& np = bk.numPhase;
469 
470  for (OCP_USI n = 0; n < bk.numBulk; n++) {
471  if (bk.bType[n] > 0) {
472  OCP_USI bId = n * np;
473  bk.flow[bk.SATNUM[n]]->CalKrPcDeriv(&bk.S[bId], &bk.kr[bId], &bk.Pc[bId],
474  &bk.dKr_dS[bId * np],
475  &bk.dPcj_dS[bId * np], n);
476  for (USI j = 0; j < np; j++)
477  bk.Pj[n * np + j] = bk.P[n] + bk.Pc[n * np + j];
478  }
479  }
480 }
481 
482 void T_FIM::CalThermalConduct(BulkConn& conn, Bulk& bk) const
483 {
484  const OCP_USI nb = bk.numBulk;
485  const OCP_USI np = bk.numPhase;
486 
487  for (OCP_USI n = 0; n < nb; n++) {
488  if (bk.bType[n] > 0) {
489  // fluid bulk
490  OCP_DBL tmp = 0;
491  for (USI j = 0; j < np; j++) {
492  tmp += bk.S[n * np + j] * bk.thconp[j];
493  bk.ktS[n * np + j] = bk.poro[n] * bk.thconp[j];
494  }
495  bk.kt[n] = bk.poro[n] * tmp + (1 - bk.poro[n]) * bk.thconr[n];
496  bk.ktP[n] = bk.poroP[n] * (tmp - bk.thconr[n]);
497  bk.ktT[n] = bk.poroT[n] * (tmp - bk.thconr[n]);
498  } else {
499  // non fluid bulk
500  bk.kt[n] = bk.thconr[n];
501  }
502  }
503 
504  OCP_USI bId, eId;
505  OCP_DBL areaB, areaE, T1, T2;
506  OCP_DBL tmpB, tmpE;
507 
508  for (OCP_USI c = 0; c < conn.numConn; c++) {
509  bId = conn.iteratorConn[c].BId();
510  eId = conn.iteratorConn[c].EId();
511  if (bk.bType[bId] > 0 && bk.bType[eId] > 0) {
512  // fluid bulk connections
513  areaB = conn.iteratorConn[c].AreaB();
514  areaE = conn.iteratorConn[c].AreaE();
515  T1 = bk.kt[bId] * areaB;
516  T2 = bk.kt[eId] * areaE;
517  conn.Adkt[c] = 1 / (1 / T1 + 1 / T2);
518 
519  tmpB = pow(conn.Adkt[c], 2) / pow(T1, 2) * areaB;
520  tmpE = pow(conn.Adkt[c], 2) / pow(T2, 2) * areaE;
521  conn.AdktP[c * 2 + 0] = tmpB * bk.ktP[bId];
522  conn.AdktP[c * 2 + 1] = tmpE * bk.ktP[eId];
523  conn.AdktT[c * 2 + 0] = tmpB * bk.ktT[bId];
524  conn.AdktT[c * 2 + 1] = tmpE * bk.ktT[eId];
525  for (USI j = 0; j < np; j++) {
526  conn.AdktS[c * np * 2 + j] = tmpB * bk.ktS[bId * np + j];
527  conn.AdktS[c * np * 2 + np + j] = tmpE * bk.ktS[eId * np + j];
528  }
529  }
530  }
531 }
532 
533 void T_FIM::CalHeatLoss(Bulk& bk, const OCP_DBL& t, const OCP_DBL& dt) const
534 {
535  bk.hLoss.CalHeatLoss(bk.bLocation, bk.T, bk.lT, bk.initT, t, dt);
536 }
537 
538 void T_FIM::ResetToLastTimeStep(Reservoir& rs, OCPControl& ctrl)
539 {
540  // Bulk
541  Bulk& bk = rs.bulk;
542 
543  // Rock
544  bk.poro = bk.lporo;
545  bk.rockVp = bk.lrockVp;
546  bk.vr = bk.lvr;
547  bk.Hr = bk.lHr;
548  // derivatives
549  bk.poroP = bk.lporoP;
550  bk.poroT = bk.lporoT;
551  bk.vrP = bk.lvrP;
552  bk.vrT = bk.lvrT;
553  bk.HrT = bk.lHrT;
554 
555  // Fluid
556  bk.phaseNum = bk.lphaseNum;
557  bk.Nt = bk.lNt;
558  bk.Ni = bk.lNi;
559  bk.vf = bk.lvf;
560  bk.T = bk.lT;
561  bk.P = bk.lP;
562  bk.Pj = bk.lPj;
563  bk.Pc = bk.lPc;
564  bk.phaseExist = bk.lphaseExist;
565  bk.S = bk.lS;
566  bk.xij = bk.lxij;
567  bk.rho = bk.lrho;
568  bk.xi = bk.lxi;
569  bk.mu = bk.lmu;
570  bk.kr = bk.lkr;
571  bk.Uf = bk.lUf;
572  bk.H = bk.lH;
573  bk.kt = bk.lkt;
574  // derivatives
575  bk.vfP = bk.lvfP;
576  bk.vfT = bk.lvfT;
577  bk.vfi = bk.lvfi;
578  bk.rhoP = bk.lrhoP;
579  bk.rhoT = bk.lrhoT;
580  bk.rhox = bk.lrhox;
581  bk.xiP = bk.lxiP;
582  bk.xiT = bk.lxiT;
583  bk.xix = bk.lxix;
584  bk.muP = bk.lmuP;
585  bk.muT = bk.lmuT;
586  bk.mux = bk.lmux;
587  bk.dPcj_dS = bk.ldPcj_dS;
588  bk.dKr_dS = bk.ldKr_dS;
589  bk.UfP = bk.lUfP;
590  bk.UfT = bk.lUfT;
591  bk.Ufi = bk.lUfi;
592  bk.HT = bk.lHT;
593  bk.Hx = bk.lHx;
594  bk.ktP = bk.lktP;
595  bk.ktT = bk.lktT;
596  bk.ktS = bk.lktS;
597  bk.dSec_dPri = bk.ldSec_dPri;
598 
599  bk.hLoss.ResetToLastTimeStep();
600 
601  BulkConn& conn = rs.conn;
602 
603  conn.Adkt = conn.lAdkt;
604  conn.AdktP = conn.lAdktP;
605  conn.AdktT = conn.lAdktT;
606  conn.AdktS = conn.lAdktS;
607 
608  // Wells
609  rs.allWells.ResetBHP();
610  rs.allWells.CalTrans(bk);
611  rs.allWells.CaldG(bk);
612  rs.allWells.CalFlux(bk);
613 
614  rs.optFeatures.ResetToLastTimeStep();
615 
616  // Iters
617  ctrl.ResetIterNRLS();
618 
619  CalRes(rs, ctrl.GetCurTime() + ctrl.GetCurDt(), ctrl.GetCurDt(), OCP_TRUE);
620 }
621 
622 void T_FIM::UpdateLastTimeStep(Reservoir& rs) const
623 {
624  // Bulk
625  Bulk& bk = rs.bulk;
626 
627  // Rock
628  bk.lporo = bk.poro;
629  bk.lrockVp = bk.rockVp;
630  bk.lvr = bk.vr;
631  bk.lHr = bk.Hr;
632  // derivatives
633  bk.lporoP = bk.poroP;
634  bk.lporoT = bk.poroT;
635  bk.lvrP = bk.vrP;
636  bk.lvrT = bk.vrT;
637  bk.lHrT = bk.HrT;
638 
639  // Fluid
640  bk.lphaseNum = bk.phaseNum;
641  bk.lNt = bk.Nt;
642  bk.lNi = bk.Ni;
643  bk.lvf = bk.vf;
644  bk.lT = bk.T;
645  bk.lP = bk.P;
646  bk.lPj = bk.Pj;
647  bk.lPc = bk.Pc;
648  bk.lphaseExist = bk.phaseExist;
649  bk.lS = bk.S;
650  bk.lxij = bk.xij;
651  bk.lrho = bk.rho;
652  bk.lxi = bk.xi;
653  bk.lmu = bk.mu;
654  bk.lkr = bk.kr;
655  bk.lUf = bk.Uf;
656  bk.lH = bk.H;
657  bk.lkt = bk.kt;
658  // derivatives
659  bk.lvfP = bk.vfP;
660  bk.lvfT = bk.vfT;
661  bk.lvfi = bk.vfi;
662  bk.lrhoP = bk.rhoP;
663  bk.lrhoT = bk.rhoT;
664  bk.lrhox = bk.rhox;
665  bk.lxiP = bk.xiP;
666  bk.lxiT = bk.xiT;
667  bk.lxix = bk.xix;
668  bk.lmuP = bk.muP;
669  bk.lmuT = bk.muT;
670  bk.lmux = bk.mux;
671  bk.ldPcj_dS = bk.dPcj_dS;
672  bk.ldKr_dS = bk.dKr_dS;
673  bk.lUfP = bk.UfP;
674  bk.lUfT = bk.UfT;
675  bk.lUfi = bk.Ufi;
676  bk.lHT = bk.HT;
677  bk.lHx = bk.Hx;
678  bk.lktP = bk.ktP;
679  bk.lktT = bk.ktT;
680  bk.lktS = bk.ktS;
681  bk.ldSec_dPri = bk.dSec_dPri;
682 
683  bk.hLoss.UpdateLastTimeStep();
684 
685  BulkConn& conn = rs.conn;
686 
687  conn.lAdkt = conn.Adkt;
688  conn.lAdktP = conn.AdktP;
689  conn.lAdktT = conn.AdktT;
690  conn.lAdktS = conn.AdktS;
691 
692  rs.allWells.UpdateLastTimeStepBHP();
693  rs.optFeatures.UpdateLastTimeStep();
694 }
695 
696 void T_FIM::CalRes(Reservoir& rs,
697  const OCP_DBL& t,
698  const OCP_DBL& dt,
699  const OCP_BOOL& resetRes0)
700 {
701  const Bulk& bk = rs.bulk;
702  const USI nb = bk.numBulk;
703  const USI np = bk.numPhase;
704  const USI nc = bk.numCom;
705  const USI len = nc + 2;
706  OCPRes& Res = bk.res;
707  BulkConn& conn = rs.conn;
708 
709  Res.SetZero();
710 
711  // Bulk to Bulk
712 
713  OCP_USI bId, eId, uId, bIdb;
714  // Accumalation Term
715  for (OCP_USI n = 0; n < nb; n++) {
716  if (bk.bType[n] > 0) {
717  // Fluid bulk
718  bId = n * len;
719  bIdb = n * nc;
720  // Volume Conservation
721  Res.resAbs[bId] = bk.rockVp[n] - bk.vf[n];
722  // Mass Conservation
723  for (USI i = 0; i < nc; i++) {
724  Res.resAbs[n * len + 1 + i] = bk.Ni[bIdb + i] - bk.lNi[bIdb + i];
725  }
726  // Energy Conservation
727  Res.resAbs[n * len + nc + 1] =
728  (bk.vf[n] * bk.Uf[n] + bk.vr[n] * bk.Hr[n]) -
729  (bk.lvf[n] * bk.lUf[n] + bk.lvr[n] * bk.lHr[n]);
730  } else {
731  // Non fluid bulk
732  Res.resAbs[n * len + nc + 1] = bk.vr[n] * bk.Hr[n] - bk.lvr[n] * bk.lHr[n];
733  }
734 
735  // Heat Loss
736  if (bk.hLoss.IfHeatLoss() && bk.bLocation[n] > 0) {
737  const OCP_DBL lambda = bk.bLocation[n] == 1 ? bk.hLoss.obD : bk.hLoss.ubD;
738  const OCP_DBL kappa = bk.bLocation[n] == 1 ? bk.hLoss.obK : bk.hLoss.ubK;
739  // dT
740  Res.resAbs[n * len + nc + 1] +=
741  dt * kappa *
742  (2 * (bk.T[n] - bk.initT[n]) / sqrt(lambda * t) - bk.hLoss.p[n]) *
743  bk.dx[n] * bk.dy[n];
744  }
745  }
746 
747  OCP_USI bId_np_j, eId_np_j, uId_np_j;
748  OCP_DBL rho, dP, dT, dNi, Akd;
749 
750  for (OCP_USI c = 0; c < conn.numConn; c++) {
751  bId = conn.iteratorConn[c].BId();
752  eId = conn.iteratorConn[c].EId();
753 
754  // Thermal conductive term all exist
755  dT = bk.T[bId] - bk.T[eId];
756  Res.resAbs[bId * len + 1 + nc] += conn.Adkt[c] * dT * dt;
757  Res.resAbs[eId * len + 1 + nc] -= conn.Adkt[c] * dT * dt;
758 
759  if (bk.bType[bId] > 0 && bk.bType[eId] > 0) {
760  // Fluid Bulk Connection
761  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
762 
763  for (USI j = 0; j < np; j++) {
764  bId_np_j = bId * np + j;
765  eId_np_j = eId * np + j;
766 
767  const OCP_BOOL exbegin = bk.phaseExist[bId_np_j];
768  const OCP_BOOL exend = bk.phaseExist[eId_np_j];
769 
770  if ((exbegin) && (exend)) {
771  rho = (bk.rho[bId_np_j] + bk.rho[eId_np_j]) / 2;
772  } else if (exbegin && (!exend)) {
773  rho = bk.rho[bId_np_j];
774  } else if ((!exbegin) && (exend)) {
775  rho = bk.rho[eId_np_j];
776  } else {
777  conn.upblock[c * np + j] = bId;
778  conn.upblock_Rho[c * np + j] = 0;
779  continue;
780  }
781 
782  uId = bId;
783  dP = (bk.Pj[bId_np_j] - GRAVITY_FACTOR * rho * bk.depth[bId]) -
784  (bk.Pj[eId_np_j] - GRAVITY_FACTOR * rho * bk.depth[eId]);
785  if (dP < 0) {
786  uId = eId;
787  }
788  conn.upblock_Rho[c * np + j] = rho;
789  conn.upblock[c * np + j] = uId;
790  uId_np_j = uId * np + j;
791 
792  if (bk.phaseExist[uId_np_j]) {
793  conn.upblock_Velocity[c * np + j] =
794  Akd * bk.kr[uId_np_j] / bk.mu[uId_np_j] * dP;
795  } else {
796  conn.upblock_Velocity[c * np + j] = 0;
797  continue;
798  }
799 
800  const OCP_DBL tmpV =
801  dt * conn.upblock_Velocity[c * np + j] * bk.xi[uId_np_j];
802 
803  for (USI i = 0; i < nc; i++) {
804  dNi = tmpV * bk.xij[uId_np_j * nc + i];
805  Res.resAbs[bId * len + 1 + i] += dNi;
806  Res.resAbs[eId * len + 1 + i] -= dNi;
807  }
808  Res.resAbs[bId * len + 1 + nc] += tmpV * bk.H[uId_np_j];
809  Res.resAbs[eId * len + 1 + nc] -= tmpV * bk.H[uId_np_j];
810  }
811  }
812  }
813 
814  // Well to Bulk
815  USI wId = nb * len;
816  for (const auto& wl : rs.allWells.wells) {
817  if (wl.IsOpen()) {
818  // Well to Bulk
819 
820  for (USI p = 0; p < wl.PerfNum(); p++) {
821  const OCP_USI k = wl.PerfLocation(p);
822  // Mass Conservation
823  for (USI i = 0; i < nc; i++) {
824  Res.resAbs[k * len + 1 + i] += wl.PerfQi_lbmol(p, i) * dt;
825  }
826  }
827 
828  if (wl.WellType() == INJ) {
829  // Energy Conservation -- Well to Bulk
830  const OCP_DBL Hw =
831  bk.flashCal[0]->CalInjWellEnthalpy(wl.InjTemp(), &wl.InjZi()[0]);
832  for (USI p = 0; p < wl.PerfNum(); p++) {
833  const OCP_USI k = wl.PerfLocation(p);
834  Res.resAbs[k * len + 1 + nc] +=
835  wl.PerfInjQt_ft3(p) * wl.PerfXi(p) * Hw * dt;
836  }
837 
838  // Well Self -- Injection
839  switch (wl.OptMode()) {
840  case BHP_MODE:
841  Res.resAbs[wId] = wl.BHP() - wl.MaxBHP();
842  break;
843  case RATE_MODE:
844  case ORATE_MODE:
845  case GRATE_MODE:
846  case WRATE_MODE:
847  case LRATE_MODE:
848  Res.resAbs[wId] = wl.MaxRate();
849  for (USI i = 0; i < nc; i++) {
850  Res.resAbs[wId] += wl.Qi_lbmol(i);
851  }
852  // if (opt.reInj) {
853  // for (auto& w : opt.connWell) {
854  // OCP_DBL tmp = 0;
855  // for (USI i = 0; i < numCom; i++) {
856  // tmp += allWell[w].qi_lbmol[i];
857  // }
858  // tmp *= opt.reInjFactor;
859  // Res.resAbs[bId] += tmp;
860  // }
861  // }
862  Res.maxWellRelRes_mol =
863  max(Res.maxWellRelRes_mol,
864  fabs(Res.resAbs[wId] / wl.MaxRate()));
865  break;
866  default:
867  OCP_ABORT("Wrong well opt mode!");
868  break;
869  }
870  } else {
871  // Energy Conservation -- Bulk to Well
872  for (USI p = 0; p < wl.PerfNum(); p++) {
873  const OCP_USI k = wl.PerfLocation(p);
874  // Mass Conservation
875  for (USI j = 0; j < np; j++) {
876  Res.resAbs[k * len + 1 + nc] += wl.PerfProdQj_ft3(p, j) *
877  bk.xi[k * np + j] *
878  bk.H[k * np + j] * dt;
879  }
880  }
881 
882  // Well Self -- Production
883  switch (wl.OptMode()) {
884  case BHP_MODE:
885  Res.resAbs[wId] = wl.BHP() - wl.MinBHP();
886  break;
887  case RATE_MODE:
888  case ORATE_MODE:
889  case GRATE_MODE:
890  case WRATE_MODE:
891  case LRATE_MODE:
892  wl.CalProdWeight(bk);
893  Res.resAbs[wId] = -wl.MaxRate();
894  for (USI i = 0; i < nc; i++) {
895  Res.resAbs[wId] += wl.Qi_lbmol(i) * wl.ProdWeight(i);
896  }
897  Res.maxWellRelRes_mol =
898  max(Res.maxWellRelRes_mol,
899  fabs(Res.resAbs[wId] / wl.MaxRate()));
900  break;
901  default:
902  OCP_ABORT("Wrong well opt mode!");
903  break;
904  }
905  }
906  wId += len;
907  }
908  }
909 
910  // Calculate RelRes
911  OCP_DBL tmp;
912  for (OCP_USI n = 0; n < nb; n++) {
913 
914  // Energy equations always exist
915  OCP_DBL eT = bk.vf[n] * bk.Uf[n] + bk.vr[n] * bk.Hr[n];
916  tmp = fabs(Res.resAbs[n * len + nc + 1] / eT);
917  if (Res.maxRelRes_E < tmp) {
918  Res.maxRelRes_E = tmp;
919  Res.maxId_E = n;
920  }
921 
922  if (bk.bType[n] > 0) {
923  // Fluid Bulk
924  for (USI i = 0; i < len - 1; i++) {
925  tmp = fabs(Res.resAbs[n * len + i] / bk.rockVp[n]);
926  if (Res.maxRelRes_V < tmp) {
927  Res.maxRelRes_V = tmp;
928  Res.maxId_V = n;
929  }
930  }
931 
932  for (USI i = 1; i < len - 1; i++) {
933  tmp = fabs(Res.resAbs[n * len + i] / bk.Nt[n]);
934  if (Res.maxRelRes_N < tmp) {
935  Res.maxRelRes_N = tmp;
936  Res.maxId_N = n;
937  }
938  }
939  }
940  }
941 
942  Dscalar(Res.resAbs.size(), -1.0, Res.resAbs.data());
943  if (resetRes0) Res.SetInitRes();
944 }
945 
946 void T_FIM::AssembleMatBulks(LinearSystem& ls,
947  const Reservoir& rs,
948  const OCP_DBL& t,
949  const OCP_DBL& dt) const
950 {
951 
952  const Bulk& bk = rs.bulk;
953  const BulkConn& conn = rs.conn;
954  const OCP_USI nb = bk.numBulk;
955  const USI np = bk.numPhase;
956  const USI nc = bk.numCom;
957  const USI ncol = nc + 2;
958  const USI ncol2 = np * nc + np;
959  const USI bsize = ncol * ncol;
960  const USI bsize2 = ncol * ncol2;
961 
962  ls.AddDim(nb);
963 
964  vector<OCP_DBL> bmat(bsize, 0);
965  // Accumulation term
966  for (USI i = 1; i < nc + 1; i++) {
967  // Mass consevation
968  bmat[i * ncol + i] = 1;
969  }
970  vector<OCP_DBL> bmatR(bmat);
971  bmatR[0] = 1;
972  for (OCP_USI n = 0; n < nb; n++) {
973  if (bk.bType[n] > 0) {
974  // Fluid Bulk
975  // Volume consevation
976  // dP
977  bmat[0] = bk.v[n] * bk.poroP[n] - bk.vfP[n];
978  // dNi
979  for (USI i = 0; i < nc; i++) {
980  bmat[i + 1] = -bk.vfi[n * nc + i];
981  }
982  // dT
983  bmat[nc + 1] = bk.v[n] * bk.poroT[n] - bk.vfT[n];
984  // Energy consevation
985  // dP
986  bmat[ncol * (ncol - 1)] =
987  bk.vfP[n] * bk.Uf[n] + bk.vf[n] * bk.UfP[n] + bk.vrP[n] * bk.Hr[n];
988  // dNi
989  for (USI i = 0; i < nc; i++) {
990  bmat[ncol * (ncol - 1) + i + 1] =
991  bk.vfi[n * nc + i] * bk.Uf[n] + bk.vf[n] * bk.Ufi[n * nc + i];
992  }
993  // dT
994  bmat[ncol * ncol - 1] = bk.vfT[n] * bk.Uf[n] + bk.vf[n] * bk.UfT[n] +
995  bk.vrT[n] * bk.Hr[n] + bk.vr[n] * bk.HrT[n];
996 
997  // Heat Loss iterm
998  if (bk.hLoss.IfHeatLoss() && bk.bLocation[n] > 0) {
999  const OCP_DBL lambda =
1000  bk.bLocation[n] == 1 ? bk.hLoss.obD : bk.hLoss.ubD;
1001  const OCP_DBL kappa =
1002  bk.bLocation[n] == 1 ? bk.hLoss.obK : bk.hLoss.ubK;
1003  // dT
1004  bmat[ncol * ncol - 1] += dt * kappa *
1005  (2 / sqrt(lambda * t) - bk.hLoss.pT[n]) *
1006  bk.dx[n] * bk.dy[n];
1007  }
1008 
1009  ls.NewDiag(n, bmat);
1010  } else {
1011  // Non Fluid Bulk
1012  // Energy consevation
1013  // dT
1014  bmatR[ncol * ncol - 1] = bk.vrT[n] * bk.Hr[n] + bk.vr[n] * bk.HrT[n];
1015 
1016  // Heat Loss iterm
1017  if (bk.hLoss.IfHeatLoss() && bk.bLocation[n] > 0) {
1018  const OCP_DBL lambda =
1019  bk.bLocation[n] == 1 ? bk.hLoss.obD : bk.hLoss.ubD;
1020  const OCP_DBL kappa =
1021  bk.bLocation[n] == 1 ? bk.hLoss.obK : bk.hLoss.ubK;
1022  // dT
1023  bmatR[ncol * ncol - 1] += dt * kappa *
1024  (2 / sqrt(lambda * t) - bk.hLoss.pT[n]) *
1025  bk.dx[n] * bk.dy[n];
1026  }
1027 
1028  ls.NewDiag(n, bmatR);
1029  }
1030  }
1031 
1032  // flux term
1033  OCP_DBL Akd;
1034  OCP_DBL transJ, transIJ, transH;
1035  vector<OCP_DBL> dFdXpB(bsize, 0); // begin bulk: dF / dXp
1036  vector<OCP_DBL> dFdXpE(bsize, 0); // end bulk: dF / dXp
1037  vector<OCP_DBL> dFdXsB(bsize2, 0); // begin bulk: dF / dXs
1038  vector<OCP_DBL> dFdXsE(bsize2, 0); // end bulk: dF / dXs
1039  OCP_DBL* dFdXpU; // up bulk: dF / dXp
1040  OCP_DBL* dFdXpD; // down bulk: dF / dXp
1041  OCP_DBL* dFdXsU; // up bulk: dF / dXs
1042  OCP_DBL* dFdXsD; // down bulk: dF / dXs
1043 
1044  OCP_USI bId, eId, uId;
1045  OCP_USI bId_np_j, eId_np_j, uId_np_j, dId_np_j;
1046  OCP_BOOL phaseExistBj, phaseExistEj, phaseExistDj;
1047  OCP_DBL xi, xij, kr, mu, rhox, xiP, xiT, xix, muP, muT, mux, H, HT, Hx;
1048  OCP_DBL dP, dT, dGamma;
1049  OCP_DBL rhoWghtU, rhoWghtD;
1050  OCP_DBL tmp;
1051 
1052  for (OCP_USI c = 0; c < conn.numConn; c++) {
1053 
1054  fill(dFdXpB.begin(), dFdXpB.end(), 0.0);
1055  fill(dFdXpE.begin(), dFdXpE.end(), 0.0);
1056  fill(dFdXsB.begin(), dFdXsB.end(), 0.0);
1057  fill(dFdXsE.begin(), dFdXsE.end(), 0.0);
1058 
1059  bId = conn.iteratorConn[c].BId();
1060  eId = conn.iteratorConn[c].EId();
1061 
1062  if (bk.bType[bId] > 0 && bk.bType[eId] > 0) {
1063  // Fluid Bulk Connection
1064  Akd = CONV1 * CONV2 * conn.iteratorConn[c].Area();
1065  dGamma = GRAVITY_FACTOR * (bk.depth[bId] - bk.depth[eId]);
1066 
1067  for (USI j = 0; j < np; j++) {
1068  uId = conn.upblock[c * np + j];
1069  uId_np_j = uId * np + j;
1070  if (!bk.phaseExist[uId_np_j]) continue;
1071  bId_np_j = bId * np + j;
1072  eId_np_j = eId * np + j;
1073  phaseExistBj = bk.phaseExist[bId_np_j];
1074  phaseExistEj = bk.phaseExist[eId_np_j];
1075 
1076  if (bId == uId) {
1077  dFdXpU = &dFdXpB[0];
1078  dFdXpD = &dFdXpE[0];
1079  dFdXsU = &dFdXsB[0];
1080  dFdXsD = &dFdXsE[0];
1081  phaseExistDj = phaseExistEj;
1082  dId_np_j = eId_np_j;
1083  } else {
1084  dFdXpU = &dFdXpE[0];
1085  dFdXpD = &dFdXpB[0];
1086  dFdXsU = &dFdXsE[0];
1087  dFdXsD = &dFdXsB[0];
1088  phaseExistDj = phaseExistBj;
1089  dId_np_j = bId_np_j;
1090  }
1091  if (phaseExistDj) {
1092  rhoWghtU = 0.5;
1093  rhoWghtD = 0.5;
1094  } else {
1095  rhoWghtU = 1;
1096  rhoWghtD = 0;
1097  }
1098 
1099  dP = bk.Pj[bId_np_j] - bk.Pj[eId_np_j] -
1100  conn.upblock_Rho[c * np + j] * dGamma;
1101  dT = bk.T[bId] - bk.T[eId];
1102  xi = bk.xi[uId_np_j];
1103  kr = bk.kr[uId_np_j];
1104  mu = bk.mu[uId_np_j];
1105  xiP = bk.xiP[uId_np_j];
1106  xiT = bk.xiT[uId_np_j];
1107  muP = bk.muP[uId_np_j];
1108  muT = bk.muT[uId_np_j];
1109  H = bk.H[uId_np_j];
1110  HT = bk.HT[uId_np_j];
1111  transJ = Akd * kr / mu;
1112 
1113  // Mass Conservation
1114  for (USI i = 0; i < nc; i++) {
1115  xij = bk.xij[uId_np_j * nc + i];
1116  transIJ = transJ * xi * xij;
1117 
1118  // dP
1119  dFdXpB[(i + 1) * ncol] += transIJ;
1120  dFdXpE[(i + 1) * ncol] -= transIJ;
1121 
1122  tmp = transJ * xiP * xij * dP;
1123  tmp += -transIJ * muP / mu * dP;
1124  dFdXpU[(i + 1) * ncol] +=
1125  (tmp - transIJ * rhoWghtU * bk.rhoP[uId_np_j] * dGamma);
1126  dFdXpD[(i + 1) * ncol] +=
1127  -transIJ * rhoWghtD * bk.rhoP[dId_np_j] * dGamma;
1128 
1129  // dT
1130  tmp = transJ * xiT * xij * dP;
1131  tmp += -transIJ * muT / mu * dP;
1132  dFdXpU[(i + 2) * ncol - 1] +=
1133  (tmp - transIJ * rhoWghtU * bk.rhoT[uId_np_j] * dGamma);
1134  dFdXpD[(i + 2) * ncol - 1] +=
1135  -transIJ * rhoWghtD * bk.rhoT[dId_np_j] * dGamma;
1136 
1137  // dS
1138  for (USI k = 0; k < np; k++) {
1139  dFdXsB[(i + 1) * ncol2 + k] +=
1140  transIJ * bk.dPcj_dS[bId_np_j * np + k];
1141  dFdXsE[(i + 1) * ncol2 + k] -=
1142  transIJ * bk.dPcj_dS[eId_np_j * np + k];
1143  dFdXsU[(i + 1) * ncol2 + k] +=
1144  Akd * bk.dKr_dS[uId_np_j * np + k] / mu * xi * xij * dP;
1145  }
1146  // dxij
1147  for (USI k = 0; k < nc; k++) {
1148  rhox = bk.rhox[uId_np_j * nc + k];
1149  xix = bk.xix[uId_np_j * nc + k];
1150  mux = bk.mux[uId_np_j * nc + k];
1151  tmp = -transIJ * rhoWghtU * rhox * dGamma;
1152  tmp += transJ * xix * xij * dP;
1153  tmp += -transIJ * mux / mu * dP;
1154  dFdXsU[(i + 1) * ncol2 + np + j * nc + k] += tmp;
1155  dFdXsD[(i + 1) * ncol2 + np + j * nc + k] +=
1156  -transIJ * rhoWghtD * bk.rhox[dId_np_j * nc + k] * dGamma;
1157  }
1158  dFdXsU[(i + 1) * ncol2 + np + j * nc + i] += transJ * xi * dP;
1159  }
1160 
1161  // Energy Conservation
1162  transH = transJ * xi * H;
1163  // dP
1164  dFdXpB[(ncol - 1) * ncol] += transH;
1165  dFdXpE[(ncol - 1) * ncol] -= transH;
1166 
1167  tmp = transJ * xiP * H * dP;
1168  tmp += -transJ * xi * muP / mu * dP * H;
1169  dFdXpU[(ncol - 1) * ncol] +=
1170  (tmp - transH * rhoWghtU * bk.rhoP[uId_np_j] * dGamma);
1171  dFdXpD[(ncol - 1) * ncol] +=
1172  -transH * rhoWghtD * bk.rhoP[dId_np_j] * dGamma;
1173 
1174  // dT
1175  tmp = transJ * xiT * H * dP;
1176  tmp += transJ * xi * HT * dP;
1177  tmp += -transH * muT / mu * dP;
1178  dFdXpU[ncol * ncol - 1] +=
1179  (tmp - transH * rhoWghtU * bk.rhoT[uId_np_j] * dGamma);
1180  dFdXpD[ncol * ncol - 1] +=
1181  -transH * rhoWghtD * bk.rhoT[dId_np_j] * dGamma;
1182 
1183  // dS
1184  for (USI k = 0; k < np; k++) {
1185  dFdXsB[(nc + 1) * ncol2 + k] +=
1186  transH * bk.dPcj_dS[bId_np_j * np + k];
1187  dFdXsE[(nc + 1) * ncol2 + k] -=
1188  transH * bk.dPcj_dS[eId_np_j * np + k];
1189  dFdXsU[(nc + 1) * ncol2 + k] +=
1190  Akd * bk.dKr_dS[uId_np_j * np + k] / mu * xi * H * dP;
1191  }
1192  // dxij
1193  for (USI k = 0; k < nc; k++) {
1194  rhox = bk.rhox[uId_np_j * nc + k];
1195  xix = bk.xix[uId_np_j * nc + k];
1196  mux = bk.mux[uId_np_j * nc + k];
1197  Hx = bk.Hx[uId_np_j * nc + k];
1198  tmp = -transH * rhoWghtU * rhox * dGamma;
1199  tmp += transJ * xix * H * dP;
1200  tmp += transJ * xi * Hx * dP;
1201  tmp += -transH * mux / mu * dP;
1202  dFdXsU[(nc + 1) * ncol2 + np + j * nc + k] += tmp;
1203  dFdXsD[(nc + 1) * ncol2 + np + j * nc + k] +=
1204  -transH * rhoWghtD * bk.rhox[dId_np_j * nc + k] * dGamma;
1205  }
1206  }
1207  }
1208 
1209  // Thermal Conduction always exist
1210  // dP
1211  dFdXpB[(ncol - 1) * ncol] += conn.AdktP[c * 2 + 0] * dT;
1212  dFdXpE[(ncol - 1) * ncol] += conn.AdktP[c * 2 + 1] * dT;
1213  // dT
1214  dFdXpB[ncol * ncol - 1] += conn.Adkt[c] + conn.AdktT[c * 2 + 0] * dT;
1215  dFdXpE[ncol * ncol - 1] += -conn.Adkt[c] + conn.AdktT[c * 2 + 1] * dT;
1216  // dS
1217  for (OCP_USI j = 0; j < np; j++) {
1218  dFdXsB[(nc + 1) * ncol2 + j] += conn.AdktS[c * np + j] * dT;
1219  dFdXsE[(nc + 1) * ncol2 + j] += conn.AdktS[c * np + np + j] * dT;
1220  }
1221 
1222  // Assemble
1223  bmat = dFdXpB;
1224  DaABpbC(ncol, ncol, ncol2, 1, dFdXsB.data(), &bk.dSec_dPri[bId * bsize2], 1,
1225  bmat.data());
1226  Dscalar(bsize, dt, bmat.data());
1227  // Begin - Begin -- add
1228  ls.AddDiag(bId, bmat);
1229  // End - Begin -- insert
1230  Dscalar(bsize, -1, bmat.data());
1231  ls.NewOffDiag(eId, bId, bmat);
1232 
1233 #ifdef OCP_NANCHECK
1234  if (!CheckNan(bmat.size(), &bmat[0])) {
1235  OCP_ABORT("INF or INF in bmat !");
1236  }
1237 #endif
1238  // End
1239  bmat = dFdXpE;
1240  DaABpbC(ncol, ncol, ncol2, 1, dFdXsE.data(), &bk.dSec_dPri[eId * bsize2], 1,
1241  bmat.data());
1242  Dscalar(bsize, dt, bmat.data());
1243  // Begin - End -- insert
1244  ls.NewOffDiag(bId, eId, bmat);
1245  // End - End -- add
1246  Dscalar(bsize, -1, bmat.data());
1247  ls.AddDiag(eId, bmat);
1248 
1249 #ifdef OCP_NANCHECK
1250  if (!CheckNan(bmat.size(), &bmat[0])) {
1251  OCP_ABORT("INF or INF in bmat !");
1252  }
1253 #endif
1254  }
1255 }
1256 
1257 void T_FIM::AssembleMatWells(LinearSystem& ls,
1258  const Reservoir& rs,
1259  const OCP_DBL& dt) const
1260 {
1261  for (auto& wl : rs.allWells.wells) {
1262  if (wl.IsOpen()) {
1263 
1264  switch (wl.WellType()) {
1265  case INJ:
1266  AssembleMatInjWells(ls, rs.bulk, wl, dt);
1267  break;
1268  case PROD:
1269  AssembleMatProdWells(ls, rs.bulk, wl, dt);
1270  break;
1271  default:
1272  OCP_ABORT("Wrong well type");
1273  }
1274  }
1275  }
1276 }
1277 
1278 void T_FIM::AssembleMatInjWells(LinearSystem& ls,
1279  const Bulk& bk,
1280  const Well& wl,
1281  const OCP_DBL& dt) const
1282 {
1283 
1284  const USI nc = bk.numCom;
1285  const USI np = bk.numPhase;
1286  const USI ncol = nc + 2;
1287  const USI ncol2 = np * nc + np;
1288  const USI bsize = ncol * ncol;
1289  const USI bsize2 = ncol * ncol2;
1290  OCP_USI n_np_j;
1291 
1292  vector<OCP_DBL> bmat(bsize, 0);
1293  vector<OCP_DBL> bmat2(bsize, 0);
1294  vector<OCP_DBL> dQdXpB(bsize, 0);
1295  vector<OCP_DBL> dQdXpW(bsize, 0);
1296  vector<OCP_DBL> dQdXsB(bsize2, 0);
1297 
1298  OCP_DBL mu, muP, muT, dP, Hw;
1299  OCP_DBL transJ, transIJ;
1300 
1301  const OCP_USI wId = ls.AddDim(1) - 1;
1302  ls.NewDiag(wId, bmat);
1303 
1304  Hw = bk.flashCal[0]->CalInjWellEnthalpy(wl.InjTemp(), &wl.InjZi()[0]);
1305 
1306  for (USI p = 0; p < wl.PerfNum(); p++) {
1307  const OCP_USI n = wl.PerfLocation(p);
1308  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
1309  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
1310  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
1311 
1312  dP = bk.P[n] - wl.BHP() - wl.DG(p);
1313 
1314  for (USI j = 0; j < np; j++) {
1315  n_np_j = n * np + j;
1316  if (!bk.phaseExist[n_np_j]) continue;
1317 
1318  mu = bk.mu[n_np_j];
1319  muP = bk.muP[n_np_j];
1320  muT = bk.muT[n_np_j];
1321 
1322  for (USI i = 0; i < nc; i++) {
1323 
1324  // Mass Conservation
1325  if (!wl.IfUseUnweight()) {
1326  transIJ = wl.PerfTransj(p, j) * wl.PerfXi(p) * wl.InjZi(i);
1327  // dQ / dP
1328  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu);
1329  dQdXpW[(i + 1) * ncol] += -transIJ;
1330 
1331  // dQ / dT
1332  dQdXpB[(i + 2) * ncol - 1] += transIJ * (-dP * muT / mu);
1333  dQdXpW[(i + 2) * ncol - 1] += 0;
1334  } else {
1335  // dQ / dP
1336  transIJ = wl.PerfTransInj(p) * wl.PerfXi(p) * wl.InjZi(i);
1337  dQdXpB[(i + 1) * ncol] += transIJ;
1338  dQdXpW[(i + 1) * ncol] += -transIJ;
1339  }
1340 
1341  if (!wl.IfUseUnweight()) {
1342  // dQ / dS
1343  for (USI k = 0; k < np; k++) {
1344  dQdXsB[(i + 1) * ncol2 + k] +=
1345  CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * wl.PerfXi(p) *
1346  wl.InjZi(i) * bk.dKr_dS[n_np_j * np + k] * dP / mu;
1347  }
1348  // dQ / dxij
1349  for (USI k = 0; k < nc; k++) {
1350  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] +=
1351  -transIJ * dP / mu * bk.mux[n_np_j * nc + k];
1352  }
1353  }
1354  }
1355 
1356  // Energy Conservation
1357  if (!wl.IfUseUnweight()) {
1358  transJ = wl.PerfTransj(p, j) * wl.PerfXi(p);
1359  // dQ / dP
1360  dQdXpB[(nc + 1) * ncol] += transJ * Hw * (1 - dP * muP / mu);
1361  dQdXpW[(nc + 1) * ncol] += -transJ * Hw;
1362 
1363  // dQ / dT
1364  dQdXpB[(nc + 2) * ncol - 1] += transJ * Hw * (-dP * muT / mu);
1365  dQdXpW[(nc + 2) * ncol - 1] += 0;
1366 
1367  // dQ / dS
1368  for (USI k = 0; k < np; k++) {
1369  dQdXsB[(nc + 1) * ncol2 + k] +=
1370  CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * wl.PerfXi(p) *
1371  bk.dKr_dS[n_np_j * np + k] * dP / mu * Hw;
1372  }
1373  // dQ / dxij
1374  for (USI k = 0; k < nc; k++) {
1375  dQdXsB[(nc + 1) * ncol2 + np + j * nc + k] +=
1376  -transJ * dP / mu * bk.mux[n_np_j * nc + k] * Hw;
1377  }
1378  } else {
1379  transJ = wl.PerfTransInj(p) * wl.PerfXi(p);
1380  // dQ / dP
1381  dQdXpB[(nc + 1) * ncol] += transJ * Hw;
1382  dQdXpW[(nc + 1) * ncol] += -transJ * Hw;
1383  }
1384 
1385  if (wl.IfUseUnweight()) break;
1386  }
1387 
1388  // Bulk to Well
1389  bmat = dQdXpB;
1390  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2], 1,
1391  bmat.data());
1392  Dscalar(bsize, dt, bmat.data());
1393  // Add
1394  ls.AddDiag(n, bmat);
1395 
1396  // Insert
1397  bmat = dQdXpW;
1398  Dscalar(bsize, dt, bmat.data());
1399  ls.NewOffDiag(n, wId, bmat);
1400 
1401  // Well
1402  switch (wl.OptMode()) {
1403  case RATE_MODE:
1404  case ORATE_MODE:
1405  case GRATE_MODE:
1406  case WRATE_MODE:
1407  case LRATE_MODE:
1408  // Diag
1409  fill(bmat.begin(), bmat.end(), 0.0);
1410  for (USI i = 0; i < nc; i++) {
1411  bmat[0] += dQdXpW[(i + 1) * ncol];
1412  bmat[(i + 1) * ncol + i + 1] = 1;
1413  }
1414  bmat[ncol * ncol - 1] = 1;
1415  ls.AddDiag(wId, bmat);
1416 
1417  // OffDiag
1418  bmat = dQdXpB;
1419  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2],
1420  1, bmat.data());
1421  fill(bmat2.begin(), bmat2.end(), 0.0);
1422  for (USI i = 0; i < nc; i++) {
1423  Daxpy(ncol, 1.0, bmat.data() + (i + 1) * ncol, bmat2.data());
1424  }
1425  ls.NewOffDiag(wId, n, bmat2);
1426  break;
1427 
1428  case BHP_MODE:
1429  // Diag
1430  fill(bmat.begin(), bmat.end(), 0.0);
1431  for (USI i = 0; i < ncol; i++) {
1432  bmat[i * ncol + i] = 1;
1433  }
1434  // Add
1435  ls.AddDiag(wId, bmat);
1436  // OffDiag
1437  fill(bmat.begin(), bmat.end(), 0.0);
1438  // Insert
1439  ls.NewOffDiag(wId, n, bmat);
1440  break;
1441 
1442  default:
1443  OCP_ABORT("Wrong Well Opt mode!");
1444  break;
1445  }
1446  }
1447 }
1448 
1449 void T_FIM::AssembleMatProdWells(LinearSystem& ls,
1450  const Bulk& bk,
1451  const Well& wl,
1452  const OCP_DBL& dt) const
1453 {
1454  const USI np = bk.numPhase;
1455  const USI nc = bk.numCom;
1456  const USI ncol = nc + 2;
1457  const USI ncol2 = np * nc + np;
1458  const USI bsize = ncol * ncol;
1459  const USI bsize2 = ncol * ncol2;
1460  OCP_USI n_np_j;
1461 
1462  vector<OCP_DBL> bmat(bsize, 0);
1463  vector<OCP_DBL> bmat2(bsize, 0);
1464  vector<OCP_DBL> dQdXpB(bsize, 0);
1465  vector<OCP_DBL> dQdXpW(bsize, 0);
1466  vector<OCP_DBL> dQdXsB(bsize2, 0);
1467 
1468  OCP_DBL xij, xi, xiP, xiT, mu, muP, muT, dP, transIJ, transJ, H, HT, Hx, tmp;
1469 
1470  const OCP_USI wId = ls.AddDim(1) - 1;
1471  ls.NewDiag(wId, bmat);
1472 
1473  // Set Prod Weight
1474  if (wl.OptMode() != BHP_MODE) wl.CalProdWeight(bk);
1475 
1476  for (USI p = 0; p < wl.PerfNum(); p++) {
1477  const OCP_USI n = wl.PerfLocation(p);
1478  fill(dQdXpB.begin(), dQdXpB.end(), 0.0);
1479  fill(dQdXpW.begin(), dQdXpW.end(), 0.0);
1480  fill(dQdXsB.begin(), dQdXsB.end(), 0.0);
1481 
1482  for (USI j = 0; j < np; j++) {
1483  n_np_j = n * np + j;
1484  if (!bk.phaseExist[n_np_j]) continue;
1485 
1486  dP = bk.Pj[n_np_j] - wl.BHP() - wl.DG(p);
1487  xi = bk.xi[n_np_j];
1488  mu = bk.mu[n_np_j];
1489  xiP = bk.xiP[n_np_j];
1490  xiT = bk.xiT[n_np_j];
1491  muP = bk.muP[n_np_j];
1492  muT = bk.muT[n_np_j];
1493  H = bk.H[n_np_j];
1494  HT = bk.HT[n_np_j];
1495 
1496  // Mass Conservation
1497  for (USI i = 0; i < nc; i++) {
1498  xij = bk.xij[n_np_j * nc + i];
1499  Hx = bk.Hx[n_np_j * nc + i];
1500  // dQ / dP
1501  transIJ = wl.PerfTransj(p, j) * xi * xij;
1502  dQdXpB[(i + 1) * ncol] += transIJ * (1 - dP * muP / mu) +
1503  dP * wl.PerfTransj(p, j) * xij * xiP;
1504  dQdXpW[(i + 1) * ncol] += -transIJ;
1505 
1506  // dQ / dT
1507  dQdXpB[(i + 2) * ncol - 1] +=
1508  transIJ * (-dP * muT / mu) + dP * wl.PerfTransj(p, j) * xij * xiT;
1509  dQdXpW[(i + 2) * ncol - 1] += 0;
1510 
1511  // dQ / dS
1512  for (USI k = 0; k < np; k++) {
1513  tmp = CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * dP / mu * xi *
1514  xij * bk.dKr_dS[n_np_j * np + k];
1515  // capillary pressure
1516  tmp += transIJ * bk.dPcj_dS[n_np_j * np + k];
1517  dQdXsB[(i + 1) * ncol2 + k] += tmp;
1518  }
1519  // dQ / dxij
1520  for (USI k = 0; k < nc; k++) {
1521  tmp = dP * wl.PerfTransj(p, j) * xij *
1522  (bk.xix[n_np_j * nc + k] - xi / mu * bk.mux[n_np_j * nc + k]);
1523  dQdXsB[(i + 1) * ncol2 + np + j * nc + k] += tmp;
1524  }
1525  dQdXsB[(i + 1) * ncol2 + np + j * nc + i] +=
1526  wl.PerfTransj(p, j) * xi * dP;
1527  }
1528 
1529  // Energy Conservation
1530  transJ = wl.PerfTransj(p, j) * xi;
1531  // dQ / dP
1532  dQdXpB[(nc + 1) * ncol] +=
1533  transJ * (1 - dP * muP / mu) * H + dP * wl.PerfTransj(p, j) * xiP * H;
1534  dQdXpW[(nc + 1) * ncol] += -transJ * H;
1535 
1536  // dQ / dT
1537  dQdXpB[(nc + 2) * ncol - 1] += transJ * (-dP * muT / mu) * H +
1538  dP * wl.PerfTransj(p, j) * xiT * H +
1539  transJ * dP * HT;
1540  dQdXpW[(nc + 2) * ncol - 1] += 0;
1541 
1542  // dQ / dS
1543  for (USI k = 0; k < np; k++) {
1544  tmp = CONV1 * wl.PerfWI(p) * wl.PerfMultiplier(p) * dP / mu * xi *
1545  bk.dKr_dS[n_np_j * np + k] * H;
1546  // capillary pressure
1547  tmp += transJ * bk.dPcj_dS[n_np_j * np + k] * H;
1548  dQdXsB[(nc + 1) * ncol2 + k] += tmp;
1549  }
1550 
1551  // dQ / dxij
1552  for (USI k = 0; k < nc; k++) {
1553  tmp =
1554  dP * wl.PerfTransj(p, j) *
1555  (bk.xix[n_np_j * nc + k] - xi / mu * bk.mux[n_np_j * nc + k]) *
1556  H +
1557  transJ * dP * Hx;
1558  dQdXsB[(nc + 1) * ncol2 + np + j * nc + k] += tmp;
1559  }
1560  }
1561 
1562  // Bulk to Well
1563  bmat = dQdXpB;
1564  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2], 1,
1565  bmat.data());
1566  Dscalar(bsize, dt, bmat.data());
1567  // Add
1568  ls.AddDiag(n, bmat);
1569  // Insert
1570  bmat = dQdXpW;
1571  Dscalar(bsize, dt, bmat.data());
1572  ls.NewOffDiag(n, wId, bmat);
1573 
1574  // Well
1575  switch (wl.OptMode()) {
1576  case RATE_MODE:
1577  case ORATE_MODE:
1578  case GRATE_MODE:
1579  case WRATE_MODE:
1580  case LRATE_MODE:
1581  // Diag
1582  fill(bmat.begin(), bmat.end(), 0.0);
1583  for (USI i = 0; i < nc; i++) {
1584  bmat[0] += dQdXpW[(i + 1) * ncol] * wl.ProdWeight(i);
1585  bmat[(i + 1) * ncol + i + 1] = 1;
1586  }
1587  bmat[ncol * ncol - 1] = 1;
1588  ls.AddDiag(wId, bmat);
1589 
1590  // OffDiag
1591  bmat = dQdXpB;
1592  DaABpbC(ncol, ncol, ncol2, 1, dQdXsB.data(), &bk.dSec_dPri[n * bsize2],
1593  1, bmat.data());
1594  fill(bmat2.begin(), bmat2.end(), 0.0);
1595  for (USI i = 0; i < nc; i++) {
1596  Daxpy(ncol, wl.ProdWeight(i), bmat.data() + (i + 1) * ncol,
1597  bmat2.data());
1598  }
1599  ls.NewOffDiag(wId, n, bmat2);
1600  break;
1601 
1602  case BHP_MODE:
1603  // Diag
1604  fill(bmat.begin(), bmat.end(), 0.0);
1605  for (USI i = 0; i < ncol; i++) {
1606  bmat[i * ncol + i] = 1;
1607  }
1608  // Add
1609  ls.AddDiag(wId, bmat);
1610  // OffDiag
1611  fill(bmat.begin(), bmat.end(), 0.0);
1612  // Insert
1613  ls.NewOffDiag(wId, n, bmat);
1614  break;
1615 
1616  default:
1617  OCP_ABORT("Wrong Well Opt mode!");
1618  break;
1619  }
1620  }
1621 }
1622 
1623 void T_FIM::GetSolution(Reservoir& rs,
1624  const vector<OCP_DBL>& u,
1625  const OCPControl& ctrl) const
1626 {
1627 
1628  // Bulk
1629  const OCP_DBL dSmaxlim = ctrl.ctrlNR.NRdSmax;
1630  // const OCP_DBL dPmaxlim = ctrl.ctrlNR.NRdPmax;
1631 
1632  Bulk& bk = rs.bulk;
1633  const OCP_USI nb = bk.numBulk;
1634  const USI np = bk.numPhase;
1635  const USI nc = bk.numCom;
1636  const USI row = np * (nc + 1);
1637  const USI col = nc + 2;
1638  vector<OCP_DBL> dtmp(row, 0);
1639  OCP_DBL chopmin = 1;
1640  OCP_DBL choptmp = 0;
1641 
1642  bk.dSNR = bk.S;
1643  bk.NRphaseNum = bk.phaseNum;
1644  bk.NRdPmax = 0;
1645  bk.NRdNmax = 0;
1646  bk.NRdTmax = 0;
1647 
1648  for (OCP_USI n = 0; n < nb; n++) {
1649  // const vector<OCP_DBL>& scm = satcm[SATNUM[n]];
1650 
1651  if (bk.bType[n] > 0) {
1652  // Fluid Bulk
1653 
1654  chopmin = 1;
1655  // compute the chop
1656  fill(dtmp.begin(), dtmp.end(), 0.0);
1657  DaAxpby(np, col, 1, &bk.dSec_dPri[n * bk.maxLendSdP], u.data() + n * col, 1,
1658  dtmp.data());
1659 
1660  for (USI j = 0; j < np; j++) {
1661  choptmp = 1;
1662  if (fabs(dtmp[j]) > dSmaxlim) {
1663  choptmp = dSmaxlim / fabs(dtmp[j]);
1664  } else if (bk.S[n * np + j] + dtmp[j] < 0.0) {
1665  choptmp = 0.9 * bk.S[n * np + j] / fabs(dtmp[j]);
1666  }
1667  chopmin = min(chopmin, choptmp);
1668  }
1669 
1670  // dS
1671  for (USI j = 0; j < np; j++) {
1672  bk.dSNRP[n * np + j] = chopmin * dtmp[j];
1673  bk.S[n * np + j] += bk.dSNRP[n * np + j];
1674  }
1675 
1676  // dP
1677  OCP_DBL dP = u[n * col];
1678  if (fabs(bk.NRdPmax) < fabs(dP)) bk.NRdPmax = dP;
1679  bk.P[n] += dP;
1680  bk.dPNR[n] = dP;
1681 
1682  // dNi
1683  bk.NRstep[n] = chopmin;
1684  for (USI i = 0; i < nc; i++) {
1685  bk.dNNR[n * nc + i] = u[n * col + 1 + i] * chopmin;
1686  if (fabs(bk.NRdNmax) < fabs(bk.dNNR[n * nc + i]) / bk.Nt[n])
1687  bk.NRdNmax = bk.dNNR[n * nc + i] / bk.Nt[n];
1688 
1689  bk.Ni[n * nc + i] += bk.dNNR[n * nc + i];
1690  }
1691  }
1692 
1693  // dT
1694  OCP_DBL dT = u[n * col + col - 1];
1695  if (fabs(bk.NRdTmax) < fabs(dT)) bk.NRdTmax = dT;
1696  bk.T[n] += dT;
1697  bk.dTNR[n] = dT;
1698  }
1699 
1700  // Well
1701  OCP_USI wId = nb * col;
1702  for (auto& wl : rs.allWells.wells) {
1703  if (wl.IsOpen()) {
1704  wl.SetBHP(wl.BHP() + u[wId]);
1705  wId += col;
1706  }
1707  }
1708 }
1709 
1710 /*----------------------------------------------------------------------------*/
1711 /* Brief Change History of This File */
1712 /*----------------------------------------------------------------------------*/
1713 /* Author Date Actions */
1714 /*----------------------------------------------------------------------------*/
1715 /* Shizhe Li Nov/10/2022 Create file */
1716 /*----------------------------------------------------------------------------*/
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
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
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
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
ThermalMethod.hpp
ThermalMethod class declaration.
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::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::AdktP
vector< OCP_DBL > AdktP
d Adkt / d P, order: connections -> bId.P -> eId.P
Definition: BulkConn.hpp:150
BulkConn::lAdkt
vector< OCP_DBL > lAdkt
last Adkt
Definition: BulkConn.hpp:147
BulkConn::AdktS
vector< OCP_DBL > AdktS
d Adkt / d S, order: connections -> bId.phase -> eId.phase
Definition: BulkConn.hpp:153
BulkConn::upblock_Velocity
vector< OCP_DBL > upblock_Velocity
Definition: BulkConn.hpp:138
BulkConn::AdktT
vector< OCP_DBL > AdktT
d Adkt / d T, order: connections -> bId.T -> eId.T
Definition: BulkConn.hpp:151
BulkConn::lAdktS
vector< OCP_DBL > lAdktS
last AdktS
Definition: BulkConn.hpp:158
BulkConn::Adkt
vector< OCP_DBL > Adkt
Thermal conductivity between neighbors.
Definition: BulkConn.hpp:140
BulkConn::iteratorConn
vector< BulkPair > iteratorConn
All connections (pair of indices) between bulks: numConn.
Definition: BulkConn.hpp:124
BulkConn::lAdktP
vector< OCP_DBL > lAdktP
last AdktP
Definition: BulkConn.hpp:156
BulkConn::lAdktT
vector< OCP_DBL > lAdktT
last AdktT
Definition: BulkConn.hpp:157
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::lkt
vector< OCP_DBL > lkt
last kt
Definition: Bulk.hpp:346
Bulk::UfP
vector< OCP_DBL > UfP
d Uf / d P: numbulk
Definition: Bulk.hpp:363
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::vrP
vector< OCP_DBL > vrP
d vr / d p, numbulk
Definition: Bulk.hpp:263
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::bLocation
vector< USI > bLocation
Location of bulk: top, bottom, side.
Definition: Bulk.hpp:206
Bulk::lporoT
vector< OCP_DBL > lporoT
last poroT.
Definition: Bulk.hpp:269
Bulk::rho
vector< OCP_DBL > rho
Mass density of phase: numPhase*numBulk.
Definition: Bulk.hpp:318
Bulk::kt
vector< OCP_DBL > kt
Coefficient of thermal diffusivity: activeGridNum.
Definition: Bulk.hpp:324
Bulk::ktP
vector< OCP_DBL > ktP
d kt / d P: numbulk
Definition: Bulk.hpp:368
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::dy
vector< OCP_DBL > dy
Size of cell in y-direction: activeGridNum.
Definition: Bulk.hpp:238
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::lHrT
vector< OCP_DBL > lHrT
Last HrT.
Definition: Bulk.hpp:272
Bulk::lHx
vector< OCP_DBL > lHx
last Hx
Definition: Bulk.hpp:391
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::Ufi
vector< OCP_DBL > Ufi
d Uf / d Ni: numCom * numBulk
Definition: Bulk.hpp:365
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::lHr
vector< OCP_DBL > lHr
Last Hr: activeGridNum.
Definition: Bulk.hpp:258
Bulk::lporoP
vector< OCP_DBL > lporoP
last poroP.
Definition: Bulk.hpp:268
Bulk::lUfP
vector< OCP_DBL > lUfP
last UfP
Definition: Bulk.hpp:387
Bulk::index_maxNRdSSP
OCP_USI index_maxNRdSSP
index of grid which has maxNRdSSP
Definition: Bulk.hpp:416
Bulk::lrhoT
vector< OCP_DBL > lrhoT
last rhoT
Definition: Bulk.hpp:377
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::vfT
vector< OCP_DBL > vfT
d vf / d T, numBulk
Definition: Bulk.hpp:350
Bulk::lHT
vector< OCP_DBL > lHT
last HT
Definition: Bulk.hpp:390
Bulk::hLoss
HeatLoss hLoss
Heat loss iterm.
Definition: Bulk.hpp:207
Bulk::lvfT
vector< OCP_DBL > lvfT
last vfT
Definition: Bulk.hpp:374
Bulk::lmuT
vector< OCP_DBL > lmuT
last muT
Definition: Bulk.hpp:383
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::Hx
vector< OCP_DBL > Hx
d Hj / d xij: numPhase * numCom * numbulk
Definition: Bulk.hpp:367
Bulk::xiT
vector< OCP_DBL > xiT
d xij / d T, numPhase * numbulk
Definition: Bulk.hpp:356
Bulk::poroT
vector< OCP_DBL > poroT
d poro / d T.
Definition: Bulk.hpp:262
Bulk::thconr
vector< OCP_DBL > thconr
Rock ifThermal conductivity: activeGridNum.
Definition: Bulk.hpp:250
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::lUfT
vector< OCP_DBL > lUfT
last UfT
Definition: Bulk.hpp:388
Bulk::UfT
vector< OCP_DBL > UfT
d Uf / d T: numbulk
Definition: Bulk.hpp:364
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::lktS
vector< OCP_DBL > lktS
last ktS
Definition: Bulk.hpp:394
Bulk::lPj
vector< OCP_DBL > lPj
last Pj
Definition: Bulk.hpp:333
Bulk::lUf
vector< OCP_DBL > lUf
last Uf
Definition: Bulk.hpp:344
Bulk::lktT
vector< OCP_DBL > lktT
last ktT
Definition: Bulk.hpp:393
Bulk::xiP
vector< OCP_DBL > xiP
d xi / d P: numPhase*numBulk.
Definition: Bulk.hpp:355
Bulk::rhoT
vector< OCP_DBL > rhoT
d rhoj / d T: numPhase * numbulk
Definition: Bulk.hpp:353
Bulk::poroP
vector< OCP_DBL > poroP
d poro / d P.
Definition: Bulk.hpp:261
Bulk::lvrP
vector< OCP_DBL > lvrP
last vrp.
Definition: Bulk.hpp:270
Bulk::lktP
vector< OCP_DBL > lktP
last ktP
Definition: Bulk.hpp:392
Bulk::vrT
vector< OCP_DBL > vrT
dvr / dT: activeGridNum.
Definition: Bulk.hpp:264
Bulk::ldKr_dS
vector< OCP_DBL > ldKr_dS
last dKr_dS
Definition: Bulk.hpp:386
Bulk::Uf
vector< OCP_DBL > Uf
Internal energy of fluid: numBulk.
Definition: Bulk.hpp:322
Bulk::lUfi
vector< OCP_DBL > lUfi
last Ufi
Definition: Bulk.hpp:389
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::NRdTmax
OCP_DBL NRdTmax
Max temperature difference in an NR step.
Definition: Bulk.hpp:418
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::lxiT
vector< OCP_DBL > lxiT
last xiT
Definition: Bulk.hpp:380
Bulk::lxix
vector< OCP_DBL > lxix
last xix
Definition: Bulk.hpp:381
Bulk::bType
vector< USI > bType
Indicate bulk type, 0: rock, 1: rock and fluid.
Definition: Bulk.hpp:205
Bulk::dPNR
vector< OCP_DBL > dPNR
P change between NR steps.
Definition: Bulk.hpp:412
Bulk::lmuP
vector< OCP_DBL > lmuP
last muP
Definition: Bulk.hpp:382
Bulk::ktT
vector< OCP_DBL > ktT
d kt / d T: activeGridNum.
Definition: Bulk.hpp:369
Bulk::lvrT
vector< OCP_DBL > lvrT
Last vrT.
Definition: Bulk.hpp:271
Bulk::ktS
vector< OCP_DBL > ktS
d kt / d S: numPhase * numbulk
Definition: Bulk.hpp:370
Bulk::muT
vector< OCP_DBL > muT
d muj / d T: numPhase * numbulk
Definition: Bulk.hpp:359
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::Hr
vector< OCP_DBL > Hr
Enthalpy of rock: activeGridNum.
Definition: Bulk.hpp:252
Bulk::v
vector< OCP_DBL > v
Volume of grids: activeGridNum.
Definition: Bulk.hpp:240
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::HT
vector< OCP_DBL > HT
d Hj / d T: numPhase * numbulk
Definition: Bulk.hpp:366
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::thconp
vector< OCP_DBL > thconp
Phase thermal conductivity: numPhase.
Definition: Bulk.hpp:172
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::dx
vector< OCP_DBL > dx
Size of cell in x-direction: activeGridNum.
Definition: Bulk.hpp:237
Bulk::H
vector< OCP_DBL > H
Enthalpy of phase: numPhase*numBulk.
Definition: Bulk.hpp:323
Bulk::dTNR
vector< OCP_DBL > dTNR
T change between NR steps.
Definition: Bulk.hpp:413
Bulk::dPcj_dS
vector< OCP_DBL > dPcj_dS
d Pcj / d Sk: numPhase * numPhase * bulk.
Definition: Bulk.hpp:361
Bulk::HrT
vector< OCP_DBL > HrT
dHr / dT: activeGridNum.
Definition: Bulk.hpp:265
Bulk::initT
vector< OCP_DBL > initT
Initial temperature of each bulk: numBulk.
Definition: Bulk.hpp:169
Bulk::lbRowSizedSdP
vector< USI > lbRowSizedSdP
last bRowSizedSdP
Definition: Bulk.hpp:493
Bulk::vr
vector< OCP_DBL > vr
Volume of rock: activeGridNum.
Definition: Bulk.hpp:251
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::lPc
vector< OCP_DBL > lPc
last Pc
Definition: Bulk.hpp:334
Bulk::lvr
vector< OCP_DBL > lvr
Last vr: activeGridNum.
Definition: Bulk.hpp:257
Bulk::maxLendSdP
USI maxLendSdP
length of dSec_dPri.
Definition: Bulk.hpp:484
Bulk::numCom
USI numCom
Number of component.
Definition: Bulk.hpp:155
Bulk::lH
vector< OCP_DBL > lH
last H
Definition: Bulk.hpp:345
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
HeatLoss::pT
vector< OCP_DBL > pT
dP / dT
Definition: Bulk.hpp:90
HeatLoss::ubD
OCP_DBL ubD
Thermal diffusivity of underburden rock.
Definition: Bulk.hpp:87
HeatLoss::ubK
OCP_DBL ubK
Thermal conductivity of underburden rock.
Definition: Bulk.hpp:83
HeatLoss::p
vector< OCP_DBL > p
Auxiliary variable.
Definition: Bulk.hpp:89
HeatLoss::obD
OCP_DBL obD
Thermal diffusivity of overburden rock.
Definition: Bulk.hpp:86
HeatLoss::obK
OCP_DBL obK
Thermal conductivity of overburden rock.
Definition: Bulk.hpp:81
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::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::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::UpdateIters
void UpdateIters()
Update the number of iterations.
Definition: OCPControl.cpp:225
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::GetCurTime
OCP_DBL GetCurTime() const
Return the current time.
Definition: OCPControl.hpp:130
OCPControl::RecordTimeLS
void RecordTimeLS(const OCP_DBL &t)
Record time used for linear solver.
Definition: OCPControl.hpp:163
OCPRes::maxId_E
OCP_INT maxId_E
index of bulk who has maxRelRes_E
Definition: OCPStructure.hpp:76
OCPRes::maxRelRes_E
OCP_DBL maxRelRes_E
maximum relative residual wrt. total energy for each bulk
Definition: OCPStructure.hpp:69
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::maxRelRes0_V
OCP_DBL maxRelRes0_V
Definition: OCPStructure.hpp:64
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::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::IsOpen
OCP_BOOL IsOpen() const
Return the state of the well, Open or Close.
Definition: Well.hpp:149
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