OpenCAEPoro  0.2.0 Sep/22/2022
A simulator for multicomponent porous media flow
OCPFluidMethod.cpp
Go to the documentation of this file.
1 
12 #include "OCPFluidMethod.hpp"
13 
15 // OCP_IMPEC
17 
18 void OCP_IMPEC::Setup(Reservoir& rs, LinearSystem& myLS, const OCPControl& ctrl)
19 {
20  // Allocate Memory of auxiliary variables for IMPEC
21  rs.AllocateAuxIMPEC();
22  // Allocate Memory of Matrix for IMPEC
23  rs.AllocateMatIMPEC(myLS);
24 
25  myLS.SetupLinearSolver(SCALARFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
26 }
27 
29 void OCP_IMPEC::InitReservoir(Reservoir& rs) const
30 {
31  rs.InitIMPEC();
32 }
33 
34 void OCP_IMPEC::Prepare(Reservoir& rs, OCP_DBL& dt)
35 {
36  rs.PrepareWell();
37  OCP_DBL cfl = rs.CalCFL(dt);
38  if (cfl > 1) dt /= (cfl + 1);
39 }
40 
41 void OCP_IMPEC::SolveLinearSystem(LinearSystem& myLS, Reservoir& rs, OCPControl& ctrl)
42 {
43 #ifdef _DEBUG
44  myLS.CheckEquation();
45 #endif // DEBUG
46 
47  myLS.AssembleMatLinearSolver();
48 
49 #ifdef _DEBUG
50  myLS.OutputLinearSystem("testA.out", "testb.out");
51 #endif // _DEBUG
52 
53  GetWallTime Timer;
54  Timer.Start();
55  int status = myLS.Solve();
56  if (status < 0) {
57  status = myLS.GetNumIters();
58  }
59 
60  ctrl.UpdateTimeLS(Timer.Stop() / 1000);
61  ctrl.UpdateIterLS(status);
62  ctrl.UpdateIterNR();
63 
64 #ifdef DEBUG
65  myLS.OutputSolution("testx.out");
66 #endif // DEBUG
67 
68  rs.GetSolutionIMPEC(myLS.GetSolution());
69  myLS.ClearData();
70 }
71 
72 bool OCP_IMPEC::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
73 {
74  OCP_DBL& dt = ctrl.current_dt;
75 
76  // first check : Pressure check
77  OCP_INT flagCheck = rs.CheckP(true, true);
78  switch (flagCheck) {
79  case 1:
80  // Negative Bulk P, well P, or Perforation P
81  dt /= 2;
82  return false;
83  case 2:
84  // Switch Well opt Mode, or close the crossflow perforation
85  dt /= 1;
86  return false;
87  default:
88  // All right
89  break;
90  }
91 
92  // Calculate Flux between bulks and between bulks and wells
93  rs.CalFLuxIMPEC();
94 
95  // second check : CFL check
96  OCP_DBL cfl = rs.CalCFL(dt);
97  if (cfl > 1) {
98  dt /= 2;
99  rs.ResetVal01IMPEC();
100  cout << "CFL is too big" << endl;
101  return false;
102  }
103 
104  rs.MassConseveIMPEC(dt);
105 
106  // third check: Ni check
107  if (!rs.CheckNi()) {
108  dt /= 2;
109  rs.ResetVal02IMPEC();
110  cout << "Negative Ni occurs\n";
111  return false;
112  }
113 
114  rs.CalVpore();
115  rs.CalFlashIMPEC();
116 
117  // fouth check: Volume error check
118  if (!rs.CheckVe(0.01)) {
119  // cout << ctrl.GetCurTime() << "Days" << "=======" << endl;
120  dt /= 2;
121  rs.ResetVal03IMPEC();
122  return false;
123  }
124 
125  rs.CalKrPc();
126  rs.CalConnFluxIMPEC();
127  // rs.allWells.ShowWellStatus(rs.bulk);
128 
129  return true;
130 }
131 
132 
133 bool OCP_IMPEC::FinishNR(const Reservoir& rs)
134 {
135  //for (USI j = 0; j < rs.bulk.numPhase; j++)
136  // cout << rs.bulk.totalPhaseNum[j] << " ";
137  //cout << endl;
138  return true;
139 }
140 
141 
142 bool OCP_IMPEC::UpdateProperty01(Reservoir& rs, OCPControl& ctrl)
143 {
144  OCP_DBL& dt = ctrl.current_dt;
145 
146  // first check : Pressure check
147  OCP_INT flagCheck = rs.CheckP();
148  switch (flagCheck) {
149  case 1:
150  // Negative Bulk P, well P, or Perforation P
151  dt /= 2;
152  return false;
153  case 2:
154  // Switch Well opt Mode, or close the crossflow perforation
155  dt /= 1;
156  return false;
157  default:
158  // All right
159  break;
160  }
161 
162  rs.CalFLuxIMPEC();
163 
164  // second check : CFL check
165  OCP_DBL cfl = rs.CalCFL(dt);
166  if (cfl > 1) {
167  dt /= 2;
168  rs.ResetVal01IMPEC();
169  cout << "CFL is too big" << endl;
170  return false;
171  }
172 
173  rs.MassConseveIMPEC(dt);
174 
175  // third check: Ni check
176  if (!rs.CheckNi()) {
177  dt /= 2;
178  rs.ResetVal02IMPEC();
179  cout << "Negative Ni occurs\n";
180  return false;
181  }
182 
183  rs.CalVpore();
184  rs.CalFlashIMPEC();
185  rs.CalKrPc();
186  rs.CalConnFluxIMPEC();
187 
188  return true;
189 }
190 
191 bool OCP_IMPEC::FinishNR01(Reservoir& rs, OCPControl& ctrl)
192 {
193  OCP_DBL& dt = ctrl.current_dt;
194  // fouth check: Volume error check
195  if (!rs.CheckVe(0.01)) {
196  // continue NR
197  rs.bulk.ResetNi();
198  rs.allWells.CalTrans(rs.bulk);
199  rs.allWells.CalFlux(rs.bulk);
200  rs.allWells.CalProdWeight(rs.bulk);
201  rs.allWells.CaldG(rs.bulk);
202  return false;
203  }
204  return true;
205 }
206 
207 void OCP_IMPEC::FinishStep(Reservoir& rs, OCPControl& ctrl)
208 {
209  rs.CalIPRT(ctrl.GetCurDt());
210  rs.CalMaxChange();
211  rs.UpdateLastStepIMPEC();
212  ctrl.CalNextTstepIMPEC(rs);
213  ctrl.UpdateIters();
214 }
215 
217 // OCP_FIM
219 
220 void OCP_FIM::Setup(Reservoir& rs, LinearSystem& myLS, const OCPControl& ctrl)
221 {
222  // Allocate Bulk and BulkConn Memory
223  rs.AllocateAuxFIM();
224  // Allocate memory for internal matrix structure
225  rs.AllocateMatFIM(myLS);
226  // Allocate memory for resiual of FIM
227  OCP_USI num = (rs.GetBulkNum() + rs.GetWellNum()) * (rs.GetComNum() + 1);
228  resFIM.res.resize(num);
229 
230  myLS.SetupLinearSolver(VECTORFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
231 }
232 
233 void OCP_FIM::InitReservoir(Reservoir& rs) const
234 {
235  rs.InitFIM();
236 }
237 
238 void OCP_FIM::Prepare(Reservoir& rs, OCP_DBL& dt)
239 {
240  rs.PrepareWell();
241  rs.CalWellFlux();
242  rs.CalResFIM(resFIM, dt);
243  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
244 }
245 
246 void OCP_FIM::AssembleMat(LinearSystem& myLS, const Reservoir& rs,
247  const OCP_DBL& dt) const
248 {
249  rs.AssembleMatFIM(myLS, dt);
250  myLS.AssembleRhs(resFIM.res);
251 }
252 
253 void OCP_FIM::SolveLinearSystem(LinearSystem& myLS, Reservoir& rs, OCPControl& ctrl) const
254 {
255 #ifdef _DEBUG
256  myLS.CheckEquation();
257 #endif // DEBUG
258 
259  myLS.AssembleMatLinearSolver();
260 
261  GetWallTime Timer;
262  Timer.Start();
263  int status = myLS.Solve();
264  if (status < 0) {
265  status = myLS.GetNumIters();
266  }
267  // cout << "LS step = " << status << endl;
268 
269 #ifdef DEBUG
270  myLS.OutputLinearSystem("testA_FIM.out", "testb_FIM.out");
271  myLS.OutputSolution("testx_FIM.out");
272  myLS.CheckSolution();
273 #endif // DEBUG
274 
275  ctrl.UpdateTimeLS(Timer.Stop() / 1000);
276  ctrl.UpdateIterLS(status);
277  ctrl.UpdateIterNR();
278 
279  rs.GetSolutionFIM(myLS.GetSolution(), ctrl.ctrlNR.NRdPmax, ctrl.ctrlNR.NRdSmax);
280  // rs.GetSolution01FIM(myLS.GetSolution());
281  // rs.PrintSolFIM(ctrl.workDir + "testPNi.out");
282  myLS.ClearData();
283 }
284 
285 bool OCP_FIM::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
286 {
287  OCP_DBL& dt = ctrl.current_dt;
288 
289 
290  // Second check: Ni check and bulk Pressure check
291  if (!rs.CheckNi() || rs.CheckP(true, false) != 0) {
292  dt *= ctrl.ctrlTime.cutFacNR;
293  rs.ResetFIM(false);
294  rs.CalResFIM(resFIM, dt);
295  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
296  cout << "Cut time step size and repeat! current dt = " << fixed << setprecision(3) << dt << " days\n";
297  return false;
298  }
299 
300  // Update reservoir properties
301  rs.CalFlashDerivFIM();
302  rs.CalKrPcDerivFIM();
303  rs.CalVpore();
304  rs.CalWellTrans();
305  rs.CalWellFlux();
306  rs.CalResFIM(resFIM, dt);
307 
308  //if (rs.bulk.NRdPmax < 1E-4) {
309  // // correct
310  // cout << "correct" << endl;
311  // rs.bulk.CorrectNi(resFIM.res);
312  // // Update reservoir properties
313  // rs.CalFlashDerivFIM();
314  // rs.CalKrPcDerivFIM();
315  // rs.CalVpore();
316  // rs.CalWellTrans();
317  // rs.CalWellFlux();
318  // rs.CalResFIM(resFIM, dt);
319  //}
320 
321  return true;
322 }
323 
324 bool OCP_FIM::FinishNR(Reservoir& rs, OCPControl& ctrl)
325 {
326  OCP_USI dSn;
327 
328  const OCP_DBL NRdSmax = rs.GetNRdSmax(dSn);
329  OCP_DBL NRdPmax = rs.GetNRdPmax();
330  const OCP_DBL NRdNmax = rs.GetNRdNmax();
331  OCP_DBL NRdSmaxP = rs.GetNRdSmaxP();
332 
333 
334 
335  if (ctrl.printLevel > 1) {
336 
337  if (true) {
338  vector<OCP_INT> totalPhaseNum(3, 0);
339  vector<OCP_INT> ltotalPhaseNum(3, 0);
340  for (OCP_USI n = 0; n < rs.GetBulkNum(); n++) {
341  if (rs.bulk.NRphaseNum[n] == 0) ltotalPhaseNum[0]++;
342  if (rs.bulk.NRphaseNum[n] == 1) ltotalPhaseNum[1]++;
343  if (rs.bulk.NRphaseNum[n] == 2) ltotalPhaseNum[2]++;
344  if (rs.bulk.phaseNum[n] == 0) totalPhaseNum[0]++;
345  if (rs.bulk.phaseNum[n] == 1) totalPhaseNum[1]++;
346  if (rs.bulk.phaseNum[n] == 2) totalPhaseNum[2]++;
347  }
348  cout << to_string(totalPhaseNum[0]) + " (" + to_string((totalPhaseNum[0] - ltotalPhaseNum[0]) * 100.0 / rs.GetBulkNum())
349  + "%) " << to_string(totalPhaseNum[1]) + " (" + to_string((totalPhaseNum[1] - ltotalPhaseNum[1]) * 100.0 / rs.GetBulkNum())
350  + "%) " << to_string(totalPhaseNum[2]) + " (" + to_string((totalPhaseNum[2] - ltotalPhaseNum[2]) * 100.0 / rs.GetBulkNum())
351  + "%) " ;
352  }
353 
354  if (true) {
355  OCP_USI eVnum = 0;
356  OCP_USI eNnum = 0;
357  for (OCP_USI n = 0; n < rs.GetBulkNum(); n++) {
358  if (rs.bulk.eV[n] < ctrl.ctrlNR.NRtol)
359  eVnum++;
360  if (rs.bulk.eN[n] < ctrl.ctrlNR.NRtol)
361  eNnum++;
362  }
363  cout << "eVnum : " << setprecision(ceil(log(rs.GetBulkNum()) / log(10))) << fixed << setw(10) << (1 - eVnum * 1.0 / rs.GetBulkNum()) * 100.0 << "% eNnum : "
364  << setw(10) << (1 - eNnum * 1.0 / rs.GetBulkNum()) * 100.0 << "%" << endl;
365  }
366 
367  const USI sp = rs.grid.GetNumDigitIJK();
368  USI tmpI, tmpJ, tmpK;
369  string tmps;
370 
371  rs.grid.GetIJKBulk(tmpI, tmpJ, tmpK, rs.bulk.index_maxNRdSSP);
372  tmps = GetIJKformat(to_string(tmpI), to_string(tmpJ), to_string(tmpK), sp);
373  cout << "### NR : " + to_string(ctrl.iterNR) + " Res: " << setprecision(2) << scientific << resFIM.maxRelRes0_v << setw(12)
374  << resFIM.maxRelRes_v << setw(12) << resFIM.maxRelRes_mol << setw(11) << resFIM.maxWellRelRes_mol
375  << setw(20) << NRdPmax << setw(11) << NRdNmax << setw(11) << NRdSmax
376  << setw(40) << sqrt(rs.bulk.NRdSSP) / rs.bulk.numBulk << setw(12) << rs.bulk.maxNRdSSP << " "
377  << rs.bulk.phaseNum[rs.bulk.index_maxNRdSSP] << " " << rs.bulk.NRphaseNum[rs.bulk.index_maxNRdSSP] << " "
378  << tmps << " " << rs.bulk.ePEC[rs.bulk.index_maxNRdSSP] << endl << endl;
379 
380 
381  rs.grid.GetIJKBulk(tmpI, tmpJ, tmpK, dSn);
382  tmps = GetIJKformat(to_string(tmpI), to_string(tmpJ), to_string(tmpK), sp);
383  cout << "S: " << tmps << setw(12) << rs.bulk.eV[dSn] << setw(12) << rs.bulk.eN[dSn]
384  << setw(12) << rs.bulk.ePEC[dSn] << setw(5) << rs.bulk.phaseNum[dSn] << setw(5) << rs.bulk.NRphaseNum[dSn] << setw(12) << rs.bulk.dPNR[dSn]
385  << setw(8) << " |" << setw(12) << rs.bulk.S[dSn * 3] << setw(12) << rs.bulk.dSNR[dSn * 3] << setw(12) << rs.bulk.dSNRP[dSn * 3]
386  << setw(12) << rs.bulk.S[dSn * 3 + 1] << setw(12) << rs.bulk.dSNR[dSn * 3 + 1] << setw(12) << rs.bulk.dSNRP[dSn * 3 + 1] << endl;
387 
388  rs.grid.GetIJKBulk(tmpI, tmpJ, tmpK, resFIM.maxId_v);
389  tmps = GetIJKformat(to_string(tmpI), to_string(tmpJ), to_string(tmpK), sp);
390  cout << "V: " << tmps << setw(12) << rs.bulk.eV[resFIM.maxId_v] << setw(12) << rs.bulk.eN[resFIM.maxId_v]
391  << setw(12) << rs.bulk.ePEC[resFIM.maxId_v] << setw(5) << rs.bulk.phaseNum[resFIM.maxId_v] << setw(5) << rs.bulk.NRphaseNum[resFIM.maxId_v] << setw(12) << rs.bulk.dPNR[resFIM.maxId_v]
392  << setw(8) << " |" << setw(12) << rs.bulk.S[resFIM.maxId_v * 3] << setw(12) << rs.bulk.dSNR[resFIM.maxId_v * 3] << setw(12) << rs.bulk.dSNRP[resFIM.maxId_v * 3]
393  << setw(12) << rs.bulk.S[resFIM.maxId_v * 3 + 1] << setw(12) << rs.bulk.dSNR[resFIM.maxId_v * 3 + 1] << setw(12) << rs.bulk.dSNRP[resFIM.maxId_v * 3 + 1] << endl;
394 
395  rs.grid.GetIJKBulk(tmpI, tmpJ, tmpK, resFIM.maxId_mol);
396  tmps = GetIJKformat(to_string(tmpI), to_string(tmpJ), to_string(tmpK), sp);
397  cout << "N: " << tmps << setw(12) << rs.bulk.eV[resFIM.maxId_mol] << setw(12) << rs.bulk.eN[resFIM.maxId_mol]
398  << setw(12) << rs.bulk.ePEC[resFIM.maxId_mol] << setw(5) << rs.bulk.phaseNum[resFIM.maxId_mol] << setw(5) << rs.bulk.NRphaseNum[resFIM.maxId_mol] << setw(12) << rs.bulk.dPNR[resFIM.maxId_mol]
399  << setw(8) << " |" << setw(12) << rs.bulk.S[resFIM.maxId_mol * 3] << setw(12) << rs.bulk.dSNR[resFIM.maxId_mol * 3] << setw(12) << rs.bulk.dSNRP[resFIM.maxId_mol * 3]
400  << setw(12) << rs.bulk.S[resFIM.maxId_mol * 3 + 1] << setw(12) << rs.bulk.dSNR[resFIM.maxId_mol * 3 + 1] << setw(12) << rs.bulk.dSNRP[resFIM.maxId_mol * 3 + 1] << endl;
401  // cout << "Res2 " << Dnorm2(resFIM.res.size(), &resFIM.res[0]) / resFIM.res.size() << " ";
402  // cout << "SdP " << Dnorm1(rs.bulk.dPNR.size(), &rs.bulk.dPNR[0]) / rs.bulk.dPNR.size() << " ";
403  // cout << "SdNi " << Dnorm1(rs.bulk.dNNR.size(), &rs.bulk.dNNR[0]) / rs.bulk.dNNR.size() << " ";
404  cout << endl;
405  // rs.ShowRes(resFIM.res);
406  }
407 
408 
409  if (ctrl.iterNR > ctrl.ctrlNR.maxNRiter) {
410  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
411  rs.ResetFIM(false);
412  rs.CalResFIM(resFIM, ctrl.current_dt);
413  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
414  ctrl.ResetIterNRLS();
415  cout << "### WARNING: NR not fully converged! Cut time step size and repeat! current dt = "
416  << fixed << setprecision(3) << ctrl.current_dt << " days\n";
417  return false;
418  }
419 
420  if (((resFIM.maxRelRes_v <= resFIM.maxRelRes0_v * ctrl.ctrlNR.NRtol ||
421  resFIM.maxRelRes_v <= ctrl.ctrlNR.NRtol ||
422  resFIM.maxRelRes_mol <= ctrl.ctrlNR.NRtol)
423  && resFIM.maxWellRelRes_mol <= ctrl.ctrlNR.NRtol) ||
424  (fabs(NRdPmax) <= ctrl.ctrlNR.NRdPmin && fabs(NRdSmax) <= ctrl.ctrlNR.NRdSmin)) {
425 
426  OCP_INT flagCheck = rs.CheckP(false, true);
427 #if DEBUG
428  if (flagCheck > 0) {
429  cout << ">> Switch well constraint: Case " << flagCheck << endl;
430  }
431 #endif
432 
433  switch (flagCheck) {
434  case 1:
435  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
436  rs.ResetFIM(true);
437  rs.CalResFIM(resFIM, ctrl.current_dt);
438  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
439  ctrl.ResetIterNRLS();
440  cout << "-----" << endl;
441  return false;
442  case 2:
443  ctrl.current_dt /= 1;
444  rs.ResetFIM(true);
445  rs.CalResFIM(resFIM, ctrl.current_dt);
446  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
447  ctrl.ResetIterNRLS();
448  cout << "-----" << endl;
449  return false;
450  default:
451  return true;
452  break;
453  }
454  } else {
455  return false;
456  }
457 }
458 
459 void OCP_FIM::FinishStep(Reservoir& rs, OCPControl& ctrl) const
460 {
461  rs.CalIPRT(ctrl.GetCurDt());
462  rs.CalMaxChange();
463  rs.UpdateLastStepFIM();
464  ctrl.CalNextTstepFIM(rs);
465  ctrl.UpdateIters();
466 }
467 
468 
470 // OCP_FIMn
472 
473 
474 void OCP_FIMn::InitReservoir(Reservoir& rs)const
475 {
476  rs.InitFIM_n();
477 }
478 
479 
481 void OCP_FIMn::AssembleMat(LinearSystem& myLS, const Reservoir& rs, const OCP_DBL& dt) const
482 {
483  rs.AssembleMatFIM_n(myLS, dt);
484  myLS.AssembleRhs(resFIM.res);
485 }
486 
488 void OCP_FIMn::SolveLinearSystem(LinearSystem& myLS, Reservoir& rs, OCPControl& ctrl) const
489 {
490 #ifdef _DEBUG
491  myLS.CheckEquation();
492 #endif // DEBUG
493 
494  myLS.AssembleMatLinearSolver();
495 
496  GetWallTime Timer;
497  Timer.Start();
498  int status = myLS.Solve();
499  if (status < 0) {
500  status = myLS.GetNumIters();
501  }
502  // cout << "LS step = " << status << endl;
503 
504 #ifdef DEBUG
505  myLS.OutputLinearSystem("testA.out", "testb.out");
506  myLS.OutputSolution("testx.out");
507  myLS.CheckSolution();
508 #endif // DEBUG
509 
510  ctrl.UpdateTimeLS(Timer.Stop() / 1000);
511  ctrl.UpdateIterLS(status);
512  ctrl.UpdateIterNR();
513 
514  rs.GetSolutionFIM_n(myLS.GetSolution(), ctrl.ctrlNR.NRdPmax, ctrl.ctrlNR.NRdSmax);
515  // rs.GetSolution01FIM(myLS.GetSolution());
516  // rs.PrintSolFIM(ctrl.workDir + "testPNi.out");
517  myLS.ClearData();
518 }
519 
521 bool OCP_FIMn::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
522 {
523  OCP_DBL& dt = ctrl.current_dt;
524 
525 
526  // Second check: Ni check and bulk Pressure check
527  if (!rs.CheckNi() || rs.CheckP(true, false) != 0) {
528  dt *= ctrl.ctrlTime.cutFacNR;
529  rs.ResetFIM(false);
530  rs.CalResFIM(resFIM, dt);
531  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
532  cout << "Cut time stepsize and repeat!\n";
533  return false;
534  }
535 
536  // Update reservoir properties
537  rs.CalFlashDerivFIM_n();
538  rs.CalKrPcDerivFIM();
539  rs.CalVpore();
540  rs.CalWellTrans();
541  rs.CalWellFlux();
542  rs.CalResFIM(resFIM, dt);
543 
544  return true;
545 }
546 
547 
549 // OCP_AIMc
551 
552 void OCP_AIMc::Setup(Reservoir& rs, LinearSystem& myLS, const OCPControl& ctrl)
553 {
554  // Allocate Bulk and BulkConn Memory
555  rs.AllocateAuxAIMc();
556  // Allocate memory for internal matrix structure
557  rs.AllocateMatFIM(myLS);
558  // Allocate memory for resiual of FIM
559  OCP_USI num = (rs.GetBulkNum() + rs.GetWellNum()) * (rs.GetComNum() + 1);
560  resFIM.res.resize(num);
561 
562  myLS.SetupLinearSolver(VECTORFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
563 
564 }
565 
566 void OCP_AIMc::InitReservoir(Reservoir& rs) const
567 {
568  rs.InitAIMc();
569 }
570 
571 void OCP_AIMc::Prepare(Reservoir& rs, OCP_DBL& dt)
572 {
573  rs.PrepareWell();
574  rs.CalWellFlux();
575  rs.CalResAIMc(resFIM, dt);
576  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
577 
578  // Set FIM Bulk
579  rs.CalCFL(dt);
580  rs.SetupWellBulk();
581  rs.SetupFIMBulk();
582  //rs.bulk.FIMBulk.resize(rs.bulk.numBulk, 0);
583  //for (OCP_USI n = 0; n < rs.bulk.numBulk; n++) {
584  // rs.bulk.FIMBulk[n] = n;
585  // rs.bulk.map_Bulk2FIM[n] = n;
586  //}
587  //rs.bulk.numFIMBulk = rs.bulk.numBulk;
588  // Calculate FIM Bulk properties
589  rs.CalFlashDerivAIMc();
590  rs.CalKrPcDerivAIMc();
591  // rs.bulk.CheckDiff();
592  rs.UpdateLastStepFIM();
593 
594  rs.bulk.ShowFIMBulk(false);
595 }
596 
597 void OCP_AIMc::AssembleMat(LinearSystem& myLS, const Reservoir& rs, const OCP_DBL& dt) const
598 {
599  rs.AssembleMatAIMc(myLS, dt);
600  myLS.AssembleRhs(resFIM.res);
601 }
602 
603 void OCP_AIMc::SolveLinearSystem(LinearSystem& myLS, Reservoir& rs, OCPControl& ctrl)
604 {
605 #ifdef _DEBUG
606  myLS.CheckEquation();
607 #endif // DEBUG
608 
609  myLS.AssembleMatLinearSolver();
610 
611  GetWallTime Timer;
612  Timer.Start();
613  int status = myLS.Solve();
614  if (status < 0) {
615  status = myLS.GetNumIters();
616  }
617  // cout << "LS step = " << status << endl;
618 
619 #ifdef DEBUG
620  myLS.OutputLinearSystem("testA.out", "testb.out");
621  myLS.OutputSolution("testx.out");
622  myLS.CheckSolution();
623 #endif // DEBUG
624 
625  ctrl.UpdateTimeLS(Timer.Stop() / 1000);
626  ctrl.UpdateIterLS(status);
627  ctrl.UpdateIterNR();
628 
629  rs.GetSolutionAIMc(myLS.GetSolution(), ctrl.ctrlNR.NRdPmax, ctrl.ctrlNR.NRdSmax);
630  // rs.GetSolution01FIM(myLS.GetSolution());
631  // rs.PrintSolFIM(ctrl.workDir + "testPNi.out");
632  myLS.ClearData();
633 }
634 
635 bool OCP_AIMc::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
636 {
637  OCP_DBL& dt = ctrl.current_dt;
638 
639  // Second check: Ni check and bulk Pressure check
640  if (!rs.CheckNi() || rs.CheckP(true, false) != 0) {
641  dt *= ctrl.ctrlTime.cutFacNR;
642  rs.ResetFIM(false);
643  rs.CalResAIMc(resFIM, dt);
644  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
645  cout << "Cut time stepsize and repeat!\n";
646  return false;
647  }
648 
649  rs.CalFlashDerivAIMc();
650  rs.CalKrPcDerivAIMc();
651  rs.CalFlashAIMc();
652  // Important, Pj must be updated with current and last Pc for IMPEC Bulk
653  rs.UpdatePj();
654 
655  rs.CalVpore();
656  rs.CalWellTrans();
657  rs.CalWellFlux();
658  rs.CalResAIMc(resFIM, dt);
659  return true;
660 }
661 
662 bool OCP_AIMc::FinishNR(Reservoir& rs, OCPControl& ctrl)
663 {
664  const OCP_DBL NRdPmax = rs.GetNRdPmax();
665  const OCP_DBL NRdNmax = rs.GetNRdNmax();
666  const OCP_DBL NRdSmax = rs.GetNRdSmaxP();
667 
668 //#ifdef _DEBUG
669  cout << "### DEBUG: Residuals = " << setprecision(3) << scientific << resFIM.maxRelRes0_v << " "
670  << resFIM.maxRelRes_v << " " << resFIM.maxRelRes_mol << " " << NRdPmax
671  << " " << NRdSmax << endl;
672  //for (OCP_USI n = 0; n < resFIM.res.size(); n++) {
673  // cout << resFIM.res[n] << endl;
674  //}
675 //#endif
676 
677  if (ctrl.iterNR > ctrl.ctrlNR.maxNRiter) {
678  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
679  rs.ResetFIM(false);
680  rs.CalResAIMc(resFIM, ctrl.current_dt);
681  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
682  ctrl.ResetIterNRLS();
683  cout << "### WARNING: NR not fully converged! Cut time step size and repeat!\n";
684  return false;
685  }
686 
687  if (((resFIM.maxRelRes_v <= resFIM.maxRelRes0_v * ctrl.ctrlNR.NRtol ||
688  resFIM.maxRelRes_v <= ctrl.ctrlNR.NRtol ||
689  resFIM.maxRelRes_mol <= ctrl.ctrlNR.NRtol)
690  && resFIM.maxWellRelRes_mol <= ctrl.ctrlNR.NRtol) ||
691  (NRdPmax <= ctrl.ctrlNR.NRdPmin && NRdSmax <= ctrl.ctrlNR.NRdSmin)) {
692 
693  OCP_INT flagCheck = rs.CheckP(false, true);
694 #if DEBUG
695  if (flagCheck > 0) {
696  cout << ">> Switch well constraint: Case " << flagCheck << endl;
697  }
698 #endif
699 
700  switch (flagCheck) {
701  case 1:
702  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
703  rs.ResetFIM(true);
704  rs.CalResAIMc(resFIM, ctrl.current_dt);
705  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
706  ctrl.ResetIterNRLS();
707  return false;
708  case 2:
709  // ctrl.current_dt /= 1;
710  rs.ResetFIM(true);
711  rs.CalResAIMc(resFIM, ctrl.current_dt);
712  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
713  ctrl.ResetIterNRLS();
714  return false;
715  default:
716  // Update IMPEC Bulk Properties
717  rs.CalFlashAIMc01();
718  rs.CalKrPcAIMc();
719  return true;
720  break;
721  }
722  }
723  else {
725  //rs.CalCFL(ctrl.current_dt);
726  //rs.SetupWellBulk();
727  //rs.SetupFIMBulk(true);
729  //rs.CalFlashDerivAIMc();
730  //rs.CalKrPcDerivAIMc();
732  //rs.bulk.ShowFIMBulk(false);
733  return false;
734  }
735 }
736 
737 
739 // OCP_AIMs
741 
742 void OCP_AIMs::Setup(Reservoir& rs, LinearSystem& myLS, const OCPControl& ctrl)
743 {
744  // Allocate Memory of auxiliary variables for AIMt
745  rs.AllocateAuxAIMs();
746  // Allocate Memory of Matrix for FIM
747  rs.AllocateMatFIM(myLS);
748  // Allocate memory for resiual of FIM
749  OCP_USI num = (rs.GetBulkNum() + rs.GetWellNum()) * (rs.GetComNum() + 1);
750  resFIM.res.resize(num);
751 
752  myLS.SetupLinearSolver(VECTORFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
753 }
754 
755 void OCP_AIMs::Prepare(Reservoir& rs, OCP_DBL& dt)
756 {
757  // for Well
758  rs.PrepareWell();
759  rs.CalWellFlux();
760 
761  // Set FIM Bulk
762  rs.SetupWellBulk();
763  rs.SetupFIMBulk();
764  rs.SetupFIMBulkBoundAIMs();
765 
766  rs.bulk.ShowFIMBulk();
767 
768  // Calculate Resiual
769  rs.CalResAIMs(resFIM, dt);
770  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
771 
772  // Calculat property of FIM Bulk
773  rs.CalFlashDerivAIM(true);
774  rs.CalKrPcDerivAIM(true);
775 
776  // Store particular property of FIM Bulk
777  rs.bulk.UpdateLastStepAIM();
778 }
779 
780 void OCP_AIMs::AssembleMat(LinearSystem& myLS, const Reservoir& rs, const OCP_DBL& dt)
781 {
782  rs.AssembleMatAIMs(myLS, resFIM.res, dt);
783  myLS.AssembleRhs(resFIM.res);
784 }
785 
786 void OCP_AIMs::SolveLinearSystem(LinearSystem& myLS, Reservoir& rs, OCPControl& ctrl)
787 {
788 #ifdef _DEBUG
789  myLS.CheckEquation();
790 #endif // DEBUG
791 
792  myLS.AssembleMatLinearSolver();
793 
794  GetWallTime Timer;
795  Timer.Start();
796  int status = myLS.Solve();
797  if (status < 0) {
798  status = myLS.GetNumIters();
799  }
800  // cout << "LS step = " << status << endl;
801 
802 #ifdef DEBUG
803  myLS.OutputLinearSystem("testA.out", "testb.out");
804  myLS.OutputSolution("testx.out");
805  myLS.CheckSolution();
806 #endif // DEBUG
807 
808  ctrl.UpdateTimeLS(Timer.Stop() / 1000);
809  ctrl.UpdateIterLS(status);
810  ctrl.UpdateIterNR();
811 
812  rs.GetSolutionAIMs(myLS.GetSolution(), ctrl.ctrlNR.NRdPmax, ctrl.ctrlNR.NRdSmax);
813  // rs.GetSolution01FIM(myLS.GetSolution());
814  // rs.PrintSolFIM(ctrl.workDir + "testPNi.out");
815  myLS.ClearData();
816 }
817 
818 bool OCP_AIMs::UpdateProperty(Reservoir& rs, OCPControl& ctrl)
819 {
820  OCP_DBL& dt = ctrl.current_dt;
821 
822  // Second check: Ni check and bulk Pressure check
823  if (!rs.CheckNi() || rs.CheckP(true, false) != 0) {
824  dt *= ctrl.ctrlTime.cutFacNR;
825  rs.ResetValAIM();
826  rs.CalResAIMs(resFIM, dt);
827  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
828  cout << "Cut time stepsize and repeat! -- 01\n";
829  return false;
830  }
831 
832  // Update reservoir properties
833  rs.CalFlashIMPEC();
834  rs.CalKrPc();
835  rs.CalFlashDerivAIM(true);
836  rs.CalKrPcDerivAIM(true);
837  rs.CalVpore();
838  rs.CalFLuxIMPEC();
839  rs.CalWellTrans();
840  rs.CalWellFlux();
841  rs.CalConnFluxIMPEC();
842  rs.CalResAIMs(resFIM, dt);
843 
844  return true;
845 }
846 
847 bool OCP_AIMs::FinishNR(Reservoir& rs, OCPControl& ctrl)
848 {
849  OCP_DBL NRdPmax = rs.GetNRdPmax();
850  OCP_DBL NRdSmax = rs.GetNRdSmaxP();
851 
852 #ifdef _DEBUG
853  cout << "### DEBUG: Residuals = " << setprecision(3) << scientific << resFIM.maxRelRes0_v << " "
854  << resFIM.maxRelRes_v << " " << resFIM.maxRelRes_mol << " " << NRdSmax
855  << " " << NRdPmax << endl;
856 #endif
857 
858  if (ctrl.iterNR > ctrl.ctrlNR.maxNRiter) {
859  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
860  rs.ResetValAIM();
861  rs.CalResAIMs(resFIM, ctrl.current_dt);
862  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
863  ctrl.ResetIterNRLS();
864  cout << "### WARNING: NR not fully converged! Cut time stepsize and repeat!\n";
865  return false;
866  }
867 
868  if (resFIM.maxRelRes_v <= resFIM.maxRelRes0_v * ctrl.ctrlNR.NRtol ||
869  resFIM.maxRelRes_v <= ctrl.ctrlNR.NRtol ||
870  resFIM.maxRelRes_mol <= ctrl.ctrlNR.NRtol ||
871  (NRdPmax <= ctrl.ctrlNR.NRdPmin && NRdSmax <= ctrl.ctrlNR.NRdSmin)) {
872 
873  OCP_INT flagCheck = rs.CheckP(false, true);
874 #if DEBUG
875  if (flagCheck > 0) {
876  cout << ">> Switch well constraint: Case " << flagCheck << endl;
877  }
878 #endif
879 
880  switch (flagCheck) {
881  case 1:
882  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
883  rs.ResetValAIM();
884  rs.CalResAIMs(resFIM, ctrl.current_dt);
885  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
886  ctrl.ResetIterNRLS();
887  return false;
888  case 2:
889  ctrl.current_dt /= 1;
890  rs.ResetValAIM();
891  rs.CalResAIMs(resFIM, ctrl.current_dt);
892  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
893  ctrl.ResetIterNRLS();
894  return false;
895  default:
896  break;
897  }
898  // Mass Conserve
899  rs.bulk.InFIMNi();
900  rs.conn.MassConserveIMPEC(rs.bulk, ctrl.current_dt);
901  rs.bulk.OutFIMNi();
902  if (!rs.CheckNi()) {
903  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
904  rs.ResetValAIM();
905  rs.CalResAIMs(resFIM, ctrl.current_dt);
906  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
907  cout << "Cut time stepsize and repeat!\n";
908  return false;
909  }
910  if (!rs.CheckVe(0.01)) {
911  // cout << ctrl.GetCurTime() << "Days" << "=======" << endl;
912 
913  rs.AddFIMBulk();
914  rs.SetupFIMBulkBoundAIMs();
915 
916  rs.bulk.ShowFIMBulk();
917 
918  // Calculate Resiual
919  rs.CalResAIMs(resFIM, ctrl.current_dt);
920  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
921 
922  // Calculat property of FIM Bulk
923  rs.CalFlashDerivAIM(true);
924  rs.CalKrPcDerivAIM(true);
925 
926  // Store particular property of FIM Bulk
927  rs.bulk.UpdateLastStepAIM();
928 
929 
930 
931  // ctrl.current_dt /= 2;
932  rs.ResetValAIM();
933  ctrl.ResetIterNRLS();
934 
935  cout << "Cut time stepsize and repeat! -- 02\n";
936  return false;
937  }
938  OCP_DBL cfl = rs.CalCFLAIM(ctrl.current_dt);
939  if (cfl > 1) {
940  cout << "CFL is too big" << endl;
941  ctrl.current_dt /= 2;
942  rs.ResetValAIM();
943  ctrl.ResetIterNRLS();
944  return false;
945  }
946  return true;
947  }
948  else {
949  return false;
950  }
951 }
952 
953 void OCP_AIMs::FinishStep(Reservoir& rs, OCPControl& ctrl)
954 {
955  // rs.GetNTQT(ctrl.current_dt);
956 
957  rs.CalIPRT(ctrl.GetCurDt());
958  rs.CalMaxChange();
959  rs.UpdateLastStepAIM();
960  ctrl.CalNextTstepIMPEC(rs);
961  ctrl.UpdateIters();
962 
963 
964 }
965 
967 // OCP_AIMt
969 
970 void OCP_AIMt::Setup(Reservoir& rs, LinearSystem& myLS, LinearSystem& myAuxLS, const OCPControl& ctrl)
971 {
972  // Allocate Memory of auxiliary variables for AIMt
973  rs.AllocateAuxAIMt();
974  // Allocate Memory of Matrix for IMPEC
975  rs.AllocateMatIMPEC(myLS);
976  // Allocate memory for internal matrix structure for local FIM
977  rs.AllocateMatAIMt(myAuxLS);
978  // Allocate memory for resiual of FIM
979  OCP_USI num = (rs.GetMaxFIMBulk() + rs.GetWellNum()) * (rs.GetComNum() + 1);
980  resFIM.res.resize(num);
981 
982  myLS.SetupLinearSolver(SCALARFASP, ctrl.GetWorkDir(), ctrl.GetLsFile());
983  myAuxLS.SetupLinearSolver(VECTORFASP, ctrl.GetWorkDir(), "./bsr.fasp");
984 }
985 
986 void OCP_AIMt::Prepare(Reservoir& rs, OCP_DBL& dt)
987 {
988  rs.PrepareWell();
989  OCP_DBL cfl = rs.CalCFL(dt);
990  if (cfl > 1) dt /= (cfl + 1);
991 
992  // setup WellbulkId
993  rs.SetupWellBulk();
994 }
995 
996 
997 bool OCP_AIMt::UpdateProperty(Reservoir& rs, OCPControl& ctrl, LinearSystem& myAuxLS)
998 {
999  OCP_DBL& dt = ctrl.current_dt;
1000 
1001  rs.CalFLuxIMPEC();
1002  rs.CalCFL(dt);
1003  rs.MassConseveIMPEC(dt);
1004 
1005  // third check: Ni check
1006  if (!rs.CheckNi()) {
1007  dt /= 2;
1008  rs.ResetVal03IMPEC();
1009  cout << "Negative Ni occurs\n";
1010  return false;
1011  }
1012 
1013  rs.CalVpore();
1014  rs.CalFlashIMPEC();
1015 
1016  // Perform FIM in local grid
1017  // Init
1018  rs.SetupFIMBulk();
1019 
1020  // cout << "FIM Bulk : " << rs.bulk.numFIMBulk << endl;
1021 
1022  //for (USI i = 0; i < rs.bulk.numFIMBulk; i++) {
1023  // cout << rs.bulk.P[rs.bulk.FIMBulk[i]] << " ";
1024  //}
1025  //cout << endl << endl;
1026 
1027  rs.bulk.FlashDerivAIM(false);
1028  rs.CalKrPcDerivAIM(false);
1029 
1030  rs.CalResAIMt(resFIM, dt);
1031  resFIM.maxRelRes0_v = resFIM.maxRelRes_v;
1032 
1033  //for (USI i = 0; i < rs.bulk.numFIMBulk; i++) {
1034  // cout << rs.bulk.P[rs.bulk.FIMBulk[i]] << " ";
1035  //}
1036  //cout << endl << endl;
1037  ctrl.iterNR = 0;
1038  // cout << ctrl.iterNR << " " << resFIM.maxRelRes0_v << endl;
1039  while (true) {
1040  rs.AssembleMatAIMt(myAuxLS, dt);
1041  myAuxLS.AssembleRhs(resFIM.res);
1042  myAuxLS.AssembleMatLinearSolver();
1043  int status = myAuxLS.Solve();
1044 
1045  rs.GetSolutionAIMt(myAuxLS.GetSolution(), ctrl.ctrlNR.NRdPmax, ctrl.ctrlNR.NRdSmax);
1046  myAuxLS.ClearData();
1047 
1048  // third check: Ni check
1049  if (!rs.CheckNi()) {
1050  dt /= 2;
1051  rs.ResetVal03IMPEC();
1052  cout << "Negative Ni occurs\n";
1053  return false;
1054  }
1055 
1056  rs.bulk.FlashDerivAIM(false);
1057  rs.CalKrPcDerivAIM(false);
1058  rs.CalVpore();
1059  rs.CalWellTrans();
1060  rs.CalWellFlux();
1061  rs.CalResAIMt(resFIM, dt);
1062 
1063  //for (USI i = 0; i < rs.bulk.numFIMBulk; i++) {
1064  // cout << fixed << rs.bulk.P[rs.bulk.FIMBulk[i]] << " ";
1065  //}
1066  //cout << scientific << resFIM.maxRelRes_v;
1067  //cout << endl << endl;
1068  OCP_DBL NRdPmax = rs.GetNRdPmax();
1069  OCP_DBL NRdSmax = rs.GetNRdSmaxP();
1070  ctrl.iterNR++;
1071  //cout << ctrl.iterNR << " " << resFIM.maxRelRes_v << " "
1072  // << resFIM.maxRelRes_mol << " " << NRdPmax << " "
1073  // << NRdSmaxP << " " << endl;
1074 
1075  if (resFIM.maxRelRes_v <= resFIM.maxRelRes0_v * ctrl.ctrlNR.NRtol ||
1076  resFIM.maxRelRes_v <= ctrl.ctrlNR.NRtol ||
1077  resFIM.maxRelRes_mol <= ctrl.ctrlNR.NRtol ||
1078  (NRdPmax <= ctrl.ctrlNR.NRdPmin && NRdSmax <= ctrl.ctrlNR.NRdSmin)) {
1079  break;
1080  }
1081  if (ctrl.iterNR > ctrl.ctrlNR.maxNRiter) {
1082  ctrl.current_dt *= ctrl.ctrlTime.cutFacNR;
1083  rs.ResetVal03IMPEC();
1084  cout << "Local FIM Failed!" << endl;
1085  return false;
1086  }
1087  }
1088 
1089  // Pressure check
1090  OCP_INT flagCheck = rs.CheckP(true, true);
1091  switch (flagCheck) {
1092  case 1:
1093  // Negative Bulk P, well P, or Perforation P
1094  dt /= 2;
1095  return false;
1096  case 2:
1097  // Switch Well opt Mode, or close the crossflow perforation
1098  dt /= 1;
1099  return false;
1100  default:
1101  // All right
1102  break;
1103  }
1104 
1105  // fouth check: Volume error check
1106  if (!rs.CheckVe(0.01)) {
1107  // cout << ctrl.GetCurTime() << "Days" << "=======" << endl;
1108  dt /= 2;
1109  rs.ResetVal03IMPEC();
1110  return false;
1111  }
1112 
1113  rs.CalKrPc();
1114  rs.CalConnFluxIMPEC();
1115 
1116  return true;
1117 }
1118 
1119 
1120 /*----------------------------------------------------------------------------*/
1121 /* Brief Change History of This File */
1122 /*----------------------------------------------------------------------------*/
1123 /* Author Date Actions */
1124 /*----------------------------------------------------------------------------*/
1125 /* Shizhe Li Nov/01/2021 Create file */
1126 /* Chensong Zhang Jan/08/2022 Update output */
1127 /*----------------------------------------------------------------------------*/
VECTORFASP
const USI VECTORFASP
Use vector linear solver in Fasp.
Definition: OCPConst.hpp:79
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:22
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:26
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:24
SCALARFASP
const USI SCALARFASP
Use scalar linear solver in Fasp.
Definition: OCPConst.hpp:78
OCP_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:25
OCPFluidMethod.hpp
Declaration of solution methods for fluid part in OpenCAEPoro.
AllWells::CalFlux
void CalFlux(const Bulk &myBulk)
Calculate volume flow rate and moles flow rate of each perforation.
Definition: AllWells.cpp:203
AllWells::CalTrans
void CalTrans(const Bulk &myBulk)
Calculate Transmissibility of Wells.
Definition: AllWells.cpp:192
AllWells::CaldG
void CaldG(const Bulk &myBulk)
Calculate dG.
Definition: AllWells.cpp:225
AllWells::CalProdWeight
void CalProdWeight(const Bulk &myBulk)
Calculate Prodweight.
Definition: AllWells.cpp:214
BulkConn::MassConserveIMPEC
void MassConserveIMPEC(Bulk &myBulk, const OCP_DBL &dt) const
Update mole composition of each bulk according to mass conservation for IMPEC.
Definition: BulkConn.cpp:426
Bulk::OutFIMNi
void OutFIMNi()
FIMNi -> Ni in FIM Bulk.
Definition: Bulk.cpp:3614
Bulk::ResetNi
void ResetNi()
Reset Ni to the ones of the last time step.
Definition: Bulk.hpp:201
Bulk::FlashDerivAIM
void FlashDerivAIM(const bool &IfAIMs)
Perform flash calculation with Ni and calculate derivatives.
Definition: Bulk.cpp:3173
Bulk::InFIMNi
void InFIMNi()
Ni in FIM Bulk -> FIMNi.
Definition: Bulk.cpp:3602
ControlNR::NRtol
OCP_DBL NRtol
Maximum non-linear convergence error.
Definition: OCPControl.hpp:67
ControlNR::NRdSmin
OCP_DBL NRdSmin
Minimum Saturation change in a Newton iteration.
Definition: OCPControl.hpp:71
ControlNR::NRdPmax
OCP_DBL NRdPmax
Maximum Pressure change in a Newton iteration.
Definition: OCPControl.hpp:68
ControlNR::maxNRiter
USI maxNRiter
Maximum number of Newton iterations in a timestep.
Definition: OCPControl.hpp:66
ControlNR::NRdSmax
OCP_DBL NRdSmax
Maximum Saturation change in a Newton iteration.
Definition: OCPControl.hpp:69
ControlNR::NRdPmin
OCP_DBL NRdPmin
Minimum Pressure change in a Newton iteration.
Definition: OCPControl.hpp:70
ControlTime::cutFacNR
OCP_DBL cutFacNR
Factor by which timestep 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
Grid::GetIJKBulk
void GetIJKBulk(USI &i, USI &j, USI &k, const OCP_USI &n) const
Return the 3D coordinate for object grid with bulk(active grids) index.
Definition: Grid.cpp:356
LinearSystem
Linear solvers for discrete systems.
Definition: LinearSystem.hpp:31
LinearSystem::OutputLinearSystem
void OutputLinearSystem(const string &fileA, const string &fileb) const
Output the mat and rhs to fileA and fileb. // TODO: output to some obj?
Definition: LinearSystem.cpp:65
LinearSystem::GetNumIters
USI GetNumIters()
Return the Max Iters.
Definition: LinearSystem.hpp:71
LinearSystem::AssembleMatLinearSolver
void AssembleMatLinearSolver()
Assemble Mat for Linear Solver.
Definition: LinearSystem.hpp:66
LinearSystem::ClearData
void ClearData()
Clear the internal matrix data for scalar-value problems.
Definition: LinearSystem.cpp:44
LinearSystem::Solve
OCP_INT Solve()
Solve the Linear System.
Definition: LinearSystem.hpp:68
LinearSystem::AssembleRhs
void AssembleRhs(const vector< OCP_DBL > &rhs)
Assign Rhs — used for FIM now.
Definition: LinearSystem.cpp:57
LinearSystem::CheckEquation
void CheckEquation() const
Check whether NAN or INF occurs in equations, used in debug mode.
Definition: LinearSystem.cpp:123
LinearSystem::GetSolution
vector< OCP_DBL > & GetSolution()
Return the solution.
Definition: LinearSystem.hpp:52
LinearSystem::SetupLinearSolver
void SetupLinearSolver(const USI &i, const string &dir, const string &file)
Setup LinearSolver.
Definition: LinearSystem.cpp:153
LinearSystem::OutputSolution
void OutputSolution(const string &filename) const
Output the solution to a disk file name.
Definition: LinearSystem.cpp:112
LinearSystem::CheckSolution
void CheckSolution() const
Check whether NAN or INF occurs in solutions, used in debug mode.
Definition: LinearSystem.cpp:142
OCPControl
All control parameters except for well controlers.
Definition: OCPControl.hpp:94
OCPControl::CalNextTstepIMPEC
void CalNextTstepIMPEC(const Reservoir &reservoir)
Calculate the next time step according to max change of some variables.
Definition: OCPControl.cpp:247
OCPControl::GetCurDt
OCP_DBL & GetCurDt()
Return current dt.
Definition: OCPControl.hpp:133
OCPControl::GetWorkDir
string GetWorkDir() const
Return work dir name.
Definition: OCPControl.hpp:179
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:324
OCPControl::ResetIterNRLS
void ResetIterNRLS()
Reset the number of iterations.
Definition: OCPControl.cpp:333
OCPControl::GetLsFile
string GetLsFile() const
Return linear solver file name.
Definition: OCPControl.hpp:182
OCPControl::UpdateIterLS
void UpdateIterLS(const USI &num)
Update the number of linear iterations.
Definition: OCPControl.hpp:154
OCPControl::UpdateTimeLS
void UpdateTimeLS(const OCP_DBL &t)
Update time used for linear solver.
Definition: OCPControl.hpp:163
OCP_AIMc::UpdateProperty
bool UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:635
OCP_AIMc::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:566
OCP_AIMc::SolveLinearSystem
void SolveLinearSystem(LinearSystem &myLS, Reservoir &rs, OCPControl &ctrl)
Solve the linear system.
Definition: OCPFluidMethod.cpp:603
OCP_AIMc::Prepare
void Prepare(Reservoir &rs, OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:571
OCP_AIMc::FinishNR
bool FinishNR(Reservoir &rs, OCPControl &ctrl)
Finish a Newton-Raphson iteration.
Definition: OCPFluidMethod.cpp:662
OCP_AIMc::AssembleMat
void AssembleMat(LinearSystem &myLS, const Reservoir &rs, const OCP_DBL &dt) const
Assemble Matrix.
Definition: OCPFluidMethod.cpp:597
OCP_AIMc::Setup
void Setup(Reservoir &rs, LinearSystem &myLS, const OCPControl &ctrl)
Setup AIMc.
Definition: OCPFluidMethod.cpp:552
OCP_AIMs::AssembleMat
void AssembleMat(LinearSystem &myLS, const Reservoir &rs, const OCP_DBL &dt)
Assemble Matrix.
Definition: OCPFluidMethod.cpp:780
OCP_AIMs::Setup
void Setup(Reservoir &rs, LinearSystem &myLS, const OCPControl &ctrl)
Setup AIMs.
Definition: OCPFluidMethod.cpp:742
OCP_AIMs::SolveLinearSystem
void SolveLinearSystem(LinearSystem &myLS, Reservoir &rs, OCPControl &ctrl)
Solve the linear system.
Definition: OCPFluidMethod.cpp:786
OCP_AIMs::Prepare
void Prepare(Reservoir &rs, OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:755
OCP_AIMs::FinishStep
void FinishStep(Reservoir &rs, OCPControl &ctrl)
Finish a time step.
Definition: OCPFluidMethod.cpp:953
OCP_AIMs::UpdateProperty
bool UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:818
OCP_AIMs::FinishNR
bool FinishNR(Reservoir &rs, OCPControl &ctrl)
Finish a Newton-Raphson iteration.
Definition: OCPFluidMethod.cpp:847
OCP_AIMt::Prepare
void Prepare(Reservoir &rs, OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:986
OCP_AIMt::UpdateProperty
bool UpdateProperty(Reservoir &rs, OCPControl &ctrl, LinearSystem &myAuxLS)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:997
OCP_AIMt::Setup
void Setup(Reservoir &rs, LinearSystem &myLS, LinearSystem &myAuxLS, const OCPControl &ctrl)
Setup AIMt.
Definition: OCPFluidMethod.cpp:970
OCP_FIM::AssembleMat
void AssembleMat(LinearSystem &myLS, const Reservoir &rs, const OCP_DBL &dt) const
Assemble Matrix.
Definition: OCPFluidMethod.cpp:246
OCP_FIM::Prepare
void Prepare(Reservoir &rs, OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:238
OCP_FIM::SolveLinearSystem
void SolveLinearSystem(LinearSystem &myLS, Reservoir &rs, OCPControl &ctrl) const
Solve the linear system.
Definition: OCPFluidMethod.cpp:253
OCP_FIM::Setup
void Setup(Reservoir &rs, LinearSystem &myLS, const OCPControl &ctrl)
Setup FIM.
Definition: OCPFluidMethod.cpp:220
OCP_FIM::FinishStep
void FinishStep(Reservoir &rs, OCPControl &ctrl) const
Finish a time step.
Definition: OCPFluidMethod.cpp:459
OCP_FIM::resFIM
ResFIM resFIM
Resiual for FIM.
Definition: OCPFluidMethod.hpp:81
OCP_FIM::UpdateProperty
bool UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:285
OCP_FIM::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:233
OCP_FIM::FinishNR
bool FinishNR(Reservoir &rs, OCPControl &ctrl)
Finish a Newton-Raphson iteration.
Definition: OCPFluidMethod.cpp:324
OCP_FIMn::SolveLinearSystem
void SolveLinearSystem(LinearSystem &myLS, Reservoir &rs, OCPControl &ctrl) const
Solve the linear system.
Definition: OCPFluidMethod.cpp:488
OCP_FIMn::UpdateProperty
bool UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:521
OCP_FIMn::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:474
OCP_FIMn::AssembleMat
void AssembleMat(LinearSystem &myLS, const Reservoir &rs, const OCP_DBL &dt) const
Assemble Matrix.
Definition: OCPFluidMethod.cpp:481
OCP_IMPEC::InitReservoir
void InitReservoir(Reservoir &rs) const
Init.
Definition: OCPFluidMethod.cpp:29
OCP_IMPEC::Setup
void Setup(Reservoir &rs, LinearSystem &myLS, const OCPControl &ctrl)
Setup IMPEC.
Definition: OCPFluidMethod.cpp:18
OCP_IMPEC::SolveLinearSystem
void SolveLinearSystem(LinearSystem &myLS, Reservoir &rs, OCPControl &ctrl)
Solve the linear system.
Definition: OCPFluidMethod.cpp:41
OCP_IMPEC::Prepare
void Prepare(Reservoir &rs, OCP_DBL &dt)
Prepare for Assembling matrix.
Definition: OCPFluidMethod.cpp:34
OCP_IMPEC::UpdateProperty
bool UpdateProperty(Reservoir &rs, OCPControl &ctrl)
Update properties of fluids.
Definition: OCPFluidMethod.cpp:72
OCP_IMPEC::FinishNR
bool FinishNR(const Reservoir &rs)
Determine if NR iteration finishes.
Definition: OCPFluidMethod.cpp:133
Reservoir::CheckVe
bool CheckVe(const OCP_DBL &Vlim) const
Check error between Fluids and Pores.
Definition: Reservoir.cpp:124
Reservoir::AllocateMatIMPEC
void AllocateMatIMPEC(LinearSystem &myLS) const
Allocate Maxmimum memory for internal Matirx for IMPEC.
Definition: Reservoir.cpp:222
Reservoir::InitFIM
void InitFIM()
Initialize the properties of Reservoir for FIM.
Definition: Reservoir.cpp:308
Reservoir::AssembleMatAIMs
void AssembleMatAIMs(LinearSystem &myLS, vector< OCP_DBL > &res, const OCP_DBL &dt) const
Assemble Matrix for AIMs.
Definition: Reservoir.cpp:601
Reservoir::AssembleMatFIM
void AssembleMatFIM(LinearSystem &myLS, const OCP_DBL &dt) const
Assemble Matrix for FIM.
Definition: Reservoir.cpp:384
Reservoir::AllocateAuxAIMt
void AllocateAuxAIMt()
Allocate memory for auxiliary variables used for AIMt.
Definition: Reservoir.cpp:517
Reservoir::UpdateLastStepFIM
void UpdateLastStepFIM()
Update value of last step for FIM.
Definition: Reservoir.cpp:365
Reservoir::AllocateAuxFIM
void AllocateAuxFIM()
Allocate memory for auxiliary variables used for FIM.
Definition: Reservoir.cpp:300
Reservoir::ResetVal01IMPEC
void ResetVal01IMPEC()
Reset Capillary Pressure, Flux for IMPEC.
Definition: Reservoir.cpp:259
Reservoir::CalKrPc
void CalKrPc()
Calculate Relative Permeability and Capillary for each Bulk.
Definition: Reservoir.cpp:74
Reservoir::UpdateLastStepAIM
void UpdateLastStepAIM()
Update value of last step for IMPEC.
Definition: Reservoir.cpp:640
Reservoir::CalFlashIMPEC
void CalFlashIMPEC()
Calculate Flash For IMPEC.
Definition: Reservoir.cpp:205
Reservoir::MassConseveIMPEC
void MassConseveIMPEC(const OCP_DBL &dt)
Calculate Ni according to Flux.
Definition: Reservoir.cpp:197
Reservoir::GetSolutionFIM
void GetSolutionFIM(const vector< OCP_DBL > &u, const OCP_DBL &dPmax, const OCP_DBL &dSmax)
Definition: Reservoir.cpp:409
Reservoir::CalFlashDerivAIM
void CalFlashDerivAIM(const bool &IfAIMs)
Calculate Flash for local FIM, some derivatives are needed.
Definition: Reservoir.cpp:537
Reservoir::GetNRdSmax
OCP_DBL GetNRdSmax(OCP_USI &index)
Return NRdSmax.
Definition: Reservoir.hpp:171
Reservoir::GetSolutionIMPEC
void GetSolutionIMPEC(const vector< OCP_DBL > &u)
Return the Solution to Reservoir Pressure for IMPEC.
Definition: Reservoir.cpp:241
Reservoir::CheckNi
bool CheckNi()
Check if abnormal Pressure occurs.
Definition: Reservoir.cpp:117
Reservoir::GetNRdSmaxP
OCP_DBL GetNRdSmaxP()
Return NRdSmaxP.
Definition: Reservoir.hpp:175
Reservoir::CalVpore
void CalVpore()
Calculate pore of Bulks.
Definition: Reservoir.cpp:67
Reservoir::PrepareWell
void PrepareWell()
Calculate Well Properties at the beginning of each time step.
Definition: Reservoir.cpp:46
Reservoir::AllocateAuxIMPEC
void AllocateAuxIMPEC()
Allocate memory for auxiliary variables used for IMPEC.
Definition: Reservoir.cpp:142
Reservoir::InitIMPEC
void InitIMPEC()
Initialize the properties of Reservoir for IMPEC.
Definition: Reservoir.cpp:150
Reservoir::GetMaxFIMBulk
OCP_USI GetMaxFIMBulk() const
Return MaxNUMFIMBulk.
Definition: Reservoir.hpp:80
Reservoir::ResetVal03IMPEC
void ResetVal03IMPEC()
Definition: Reservoir.cpp:275
Reservoir::GetWellNum
USI GetWellNum() const
Return the num of Well.
Definition: Reservoir.hpp:82
Reservoir::CalKrPcDerivAIM
void CalKrPcDerivAIM(const bool &IfAIMs)
Calculate Relative Permeability and Capillary and some derivatives for each Bulk.
Definition: Reservoir.cpp:543
Reservoir::GetSolutionAIMt
void GetSolutionAIMt(const vector< OCP_DBL > &u, const OCP_DBL &dPmax, const OCP_DBL &dSmax)
Definition: Reservoir.cpp:570
Reservoir::UpdateLastStepIMPEC
void UpdateLastStepIMPEC()
Update value of last step for IMPEC.
Definition: Reservoir.cpp:212
Reservoir::GetNRdPmax
OCP_DBL GetNRdPmax()
Return NRdPmax.
Definition: Reservoir.hpp:169
Reservoir::CalMaxChange
void CalMaxChange()
Calculate Maximum Change of some reference variables for IMPEC.
Definition: Reservoir.cpp:81
Reservoir::AllocateMatFIM
void AllocateMatFIM(LinearSystem &myLS) const
Allocate Maxmimum memory for internal Matirx for FIM.
Definition: Reservoir.cpp:373
Reservoir::CheckP
OCP_INT CheckP(const bool &bulkCheck=true, const bool &wellCheck=true)
Check if abnormal Pressure occurs.
Definition: Reservoir.cpp:98
Reservoir::AssembleMatAIMc
void AssembleMatAIMc(LinearSystem &myLS, const OCP_DBL &dt) const
Assemble Matrix for AIMc.
Definition: Reservoir.cpp:659
Reservoir::ResetFIM
void ResetFIM(const bool &flag)
Reset FIM.
Definition: Reservoir.cpp:471
Reservoir::GetComNum
USI GetComNum() const
Return the num of Components.
Definition: Reservoir.hpp:84
Reservoir::CalFlashDerivAIMc
void CalFlashDerivAIMc()
Calculate Flash for local FIM, some derivatives are needed.
Definition: Reservoir.cpp:699
Reservoir::CalConnFluxIMPEC
void CalConnFluxIMPEC()
Calculate flux between bulks.
Definition: Reservoir.cpp:190
Reservoir::CalWellTrans
void CalWellTrans()
Calculate Trans of Wells.
Definition: Reservoir.cpp:60
Reservoir::AllocateMatAIMt
void AllocateMatAIMt(LinearSystem &myLS) const
Allocate Maxmimum memory for internal Matirx for local FIM.
Definition: Reservoir.cpp:527
Reservoir::AllocateAuxAIMc
void AllocateAuxAIMc()
Allocate memory for auxiliary variables used for FIM.
Definition: Reservoir.cpp:650
Reservoir::AssembleMatAIMt
void AssembleMatAIMt(LinearSystem &myLS, const OCP_DBL &dt) const
Assemble Matrix for AIMt -— local FIM here.
Definition: Reservoir.cpp:563
Reservoir::SetupFIMBulk
void SetupFIMBulk(const bool &NRflag=false)
Setup FIMBulk.
Definition: Reservoir.hpp:189
Reservoir::AllocateAuxAIMs
void AllocateAuxAIMs()
Allocate memory for auxiliary variables used for AIMs.
Definition: Reservoir.cpp:578
Reservoir::CalIPRT
void CalIPRT(const OCP_DBL &dt)
Calculate num of Injection, Production.
Definition: Reservoir.cpp:89
Reservoir::CalCFL
OCP_DBL CalCFL(const OCP_DBL &dt)
Calcluate the CFL number, including bulks and wells for IMPEC.
Definition: Reservoir.cpp:170
Reservoir::GetNRdNmax
OCP_DBL GetNRdNmax()
Return NRdNmax.
Definition: Reservoir.hpp:173
Reservoir::ResetVal02IMPEC
void ResetVal02IMPEC()
Reset Capillary Pressure, Moles of Componnets, Flux for IMPEC.
Definition: Reservoir.cpp:266
Reservoir::CalKrPcDerivAIMc
void CalKrPcDerivAIMc()
Calculate Relative Permeability and Capillary and some derivatives for each Bulk.
Definition: Reservoir.cpp:706
Reservoir::CalKrPcDerivFIM
void CalKrPcDerivFIM()
Calculate Relative Permeability and Capillary and some derivatives for each Bulk.
Definition: Reservoir.cpp:358
Reservoir::CalFlashDerivFIM
void CalFlashDerivFIM()
Calculate Flash for FIM, some derivatives are needed.
Definition: Reservoir.cpp:344
Reservoir::CalFLuxIMPEC
void CalFLuxIMPEC()
Calculate flux between bulks, bulks and wells.
Definition: Reservoir.cpp:182
Reservoir::GetBulkNum
OCP_USI GetBulkNum() const
Return the num of Bulk.
Definition: Reservoir.hpp:78
Reservoir::CalResAIMs
void CalResAIMs(ResFIM &resFIM, const OCP_DBL &dt)
Calculate the Resiual for AIMs, it's also RHS of Linear System.
Definition: Reservoir.cpp:588
Reservoir::CalResAIMc
void CalResAIMc(ResFIM &resFIM, const OCP_DBL &dt)
Calculate the Resiual for FIM, it's also RHS of Linear System.
Definition: Reservoir.cpp:668
Reservoir::CalWellFlux
void CalWellFlux()
Calculate Flux between Bulk and Wells.
Definition: Reservoir.cpp:53
Reservoir::CalResFIM
void CalResFIM(ResFIM &resFIM, const OCP_DBL &dt)
Calculate the Resiual for FIM, it's also RHS of Linear System.
Definition: Reservoir.cpp:434
Reservoir::CalResAIMt
void CalResAIMt(ResFIM &resFIM, const OCP_DBL &dt)
Calculate the Resiual for local FIM, it's also RHS of Linear System.
Definition: Reservoir.cpp:549
Generated on Sat Oct 15 2022 00:08:52 for OpenCAEPoro by doxygen 1.9.1