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