OpenCAEPoro  0.2.0 Sep/22/2022
A simulator for multicomponent porous media flow
Public Member Functions | List of all members
BulkConn Class Reference

Properties and operations on connections between bulks (active grids). More...

#include <BulkConn.hpp>

Public Member Functions

 BulkConn ()=default
 Default constructor.
 
void Setup (const Grid &myGrid, const Bulk &myBulk)
 Setup active connections and calculate necessary properties using Grid and Bulk. More...
 
void SetupWellBulk_K (Bulk &myBulk) const
 Setup k-neighbor for bulks.
 
void AllocateMat (LinearSystem &myLS) const
 Allocate memory for the coefficient matrix. More...
 
void SetupMatSparsity (LinearSystem &myLS) const
 Setup sparsity pattern of the coefficient matrix.
 
void UpdateLastStep ()
 Update physcial values of the previous step.
 
void Reset ()
 Reset physcial values of the current step with the previous step.
 
void CheckDiff () const
 Check differences between the current and previous steps.
 
OCP_USI GetBulkNum () const
 Return number of bulks.
 
void PrintConnectionInfo (const Grid &myGrid) const
 Print information of connections on screen.
 
void PrintConnectionInfoCoor (const Grid &myGrid) const
 
void AllocateAuxIMPEC (const USI &np)
 Allocate memory for auxiliary variables used by the IMPEC method.
 
void AssembleMatIMPEC (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 Assmeble coefficient matrix for IMPEC, terms related to bulks only.
 
void CalCFL (const Bulk &myBulk, const OCP_DBL &dt) const
 Calculate the CFL number for flow between bulks???
 
void CalFluxIMPEC (const Bulk &myBulk)
 Calculate flux information about flow between bulks for IMPEC.
 
void MassConserveIMPEC (Bulk &myBulk, const OCP_DBL &dt) const
 Update mole composition of each bulk according to mass conservation for IMPEC.
 
void AllocateAuxFIM (const USI &np)
 Allocate memory for auxiliary variables used by the FIM method.
 
void AssembleMat_FIM (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 Assmeble coefficient matrix for FIM, terms related to bulks only.
 
void CalFluxFIM (const Bulk &myBulk)
 Calculate flux for FIM, considering upwinding.
 
void CalResFIM (vector< OCP_DBL > &res, const Bulk &myBulk, const OCP_DBL &dt)
 Calculate resiual for the Newton iteration in FIM.
 
void CalFluxFIMS (const Bulk &myBulk)
 rho = (S1*rho1 + S2*rho2)/(S1+S2)
 
void CalResFIMS (vector< OCP_DBL > &res, const Bulk &myBulk, const OCP_DBL &dt)
 
void AssembleMat_FIM_new (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 
void AssembleMat_FIM_newS (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 OCP_NEW_FIM rho = (S1*rho1 + S2*rho2)/(S1+S2)
 
void AssembleMat_FIM_new_n (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 OCP_NEW_FIMn.
 
void SetupFIMBulk (Bulk &myBulk, const bool &NRflag=false) const
 
void AddFIMBulk (Bulk &myBulk)
 
void SetupFIMBulkBoundAIMs (Bulk &myBulk)
 
void AllocateAuxAIMt ()
 Allocate memory for auxiliary variables used by the AIMt method.
 
void SetupMatSparsityAIMt (LinearSystem &myLS, const Bulk &myBulk) const
 Setup sparsity pattern of the coefficient matrix for AIMt.
 
void AssembleMat_AIMt (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 Assmeble coefficient matrix for FIM, terms related to bulks only.
 
void CalResAIMt (vector< OCP_DBL > &res, const Bulk &myBulk, const OCP_DBL &dt)
 Calculate resiual for the Newton iteration in local FIM.
 
void CalResAIMs (vector< OCP_DBL > &res, const Bulk &myBulk, const OCP_DBL &dt)
 
void AssembleMat_AIMs (LinearSystem &myLS, vector< OCP_DBL > &res, const Bulk &myBulk, const OCP_DBL &dt) const
 
void AllocateAuxAIMc (const USI &np)
 Allocate memory for auxiliary variables used by the AIMc method.
 
void AssembleMat_AIMc (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 
void AssembleMat_AIMc01 (LinearSystem &myLS, const Bulk &myBulk, const OCP_DBL &dt) const
 
void CalResAIMc (vector< OCP_DBL > &res, const Bulk &myBulk, const OCP_DBL &dt)
 Calculate resiual for the Newton iteration in FIM.
 

Detailed Description

Properties and operations on connections between bulks (active grids).

Definition at line 57 of file BulkConn.hpp.

Member Function Documentation

◆ AllocateMat()

void BulkConn::AllocateMat ( LinearSystem &  myLS) const

Allocate memory for the coefficient matrix.

This method should be called only once at the beginning.

Definition at line 107 of file BulkConn.cpp.

108 {
109  OCP_FUNCNAME;
110 
111  for (OCP_USI n = 0; n < numBulk; n++) {
112  MySolver.rowCapacity[n] += neighborNum[n];
113  }
114 }
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:24
OCP_FUNCNAME
#define OCP_FUNCNAME
Print Function Name.
Definition: UtilError.hpp:73

References OCP_FUNCNAME.

◆ AssembleMat_AIMs()

void BulkConn::AssembleMat_AIMs ( LinearSystem &  myLS,
vector< OCP_DBL > &  res,
const Bulk &  myBulk,
const OCP_DBL &  dt 
) const

Assmeble coefficient matrix for AIMs, terms related to bulks only parts related to FIM A, IMPEC A, and IMPEC b

Definition at line 2573 of file BulkConn.cpp.

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 }
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
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_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:25
CONV2
const OCP_DBL CONV2
Darcy constant in field unit.
Definition: OCPConst.hpp:60

References CONV1, CONV2, DaABpbC(), Dscalar(), and GRAVITY_FACTOR.

◆ AssembleMat_FIM_new()

void BulkConn::AssembleMat_FIM_new ( LinearSystem &  myLS,
const Bulk &  myBulk,
const OCP_DBL &  dt 
) const

Assmeble coefficient matrix for FIM, terms related to bulks only. OCP_NEW_FIM

Definition at line 1002 of file BulkConn.cpp.

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 }
CheckNan
bool CheckNan(const int &N, const T *x)
check NaN
Definition: DenseMat.hpp:132
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47

References CheckNan(), CONV1, CONV2, DaABpbC(), Dscalar(), GRAVITY_FACTOR, OCP_ABORT, and OCP_FUNCNAME.

◆ CalResAIMs()

void BulkConn::CalResAIMs ( vector< OCP_DBL > &  res,
const Bulk &  myBulk,
const OCP_DBL &  dt 
)

Calculate resiual for the Newton iteration in AIMs. Only parts using local FIM are considered.

Definition at line 2477 of file BulkConn.cpp.

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 }

References CONV1, CONV2, GRAVITY_FACTOR, and OCP_FUNCNAME.

◆ Setup()

void BulkConn::Setup ( const Grid &  myGrid,
const Bulk &  myBulk 
)

Setup active connections and calculate necessary properties using Grid and Bulk.

It should be called after Grid and Bulk Setup.

Definition at line 25 of file BulkConn.cpp.

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 }
BulkPair
Connection between two bulks (BId, EId); usually, indices BId > EId.
Definition: BulkConn.hpp:30
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

References GB_Pair::GetId(), GB_Pair::IsAct(), and OCP_FUNCNAME.

The documentation for this class was generated from the following files:
Generated on Sat Oct 15 2022 00:08:52 for OpenCAEPoro by doxygen 1.9.1