OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

erge branch 'master' into FOM_dev

Browse Source
This commit is contained in:
James McClure 2021-08-30 11:49:38 -04:00
commit f908e6f8ee
11 changed files with 790 additions and 305 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
common/ScaLBL.h
View File

@ -88,12 +88,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros, double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np);

173
cpu/GreyscaleColor.cpp
View File

@ -1341,7 +1341,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl
//CP: capillary penalty
// also turn off recoloring for grey nodes
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta,
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
@ -1375,6 +1375,11 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
/* Corey model parameters */
double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB;
double GreyDiff; // grey diffusion
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
const double mrt_V3=0.04761904761904762;
@ -1394,14 +1399,16 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
nB = Den[Np + n];
porosity = Poros[n];
perm = Perm[n];
GreyDiff = Perm[n];
perm = 1.0;
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
Kn_grey = GreyKn[n];
Kw_grey = GreyKw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
// local density
rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
// local relaxation time
@ -1411,6 +1418,22 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey;
Swn = 0.0;
}
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey);
Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero
// recompute the effective permeability
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5));
mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey;
Swn = 1.0;
}
// Get the 1D index based on regular data layout
ijk = Map[n];
// COMPUTE THE COLOR GRADIENT
@ -2053,21 +2076,28 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
nAB = 1.0/(nA+nB);
Aq[n] = 0.3333333333333333*nA;
Bq[n] = 0.3333333333333333*nB;
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
jA = nA*ux;
jB = nB*ux;
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
// q = 1
//nread = neighborList[n+Np];
@ -2080,17 +2110,24 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//...............................................
// Cq = {0,1,0}
jA = nA*uy;
jB = nB*uy;
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
// q = 3
//nread = neighborList[n+3*Np];
@ -2104,17 +2141,25 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//...............................................
// q = 4
// Cq = {0,0,1}
jA = nA*uz;
jB = nB*uz;
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
// q = 5
//nread = neighborList[n+5*Np];
@ -2131,7 +2176,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map
//CP: capillary penalty
// also turn off recoloring for grey nodes
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
@ -2152,6 +2197,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
/* Corey model parameters */
double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB;
double GreyDiff; // grey diffusion
//double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002)
double porosity;
double perm;//voxel permeability
@ -2180,11 +2230,14 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
nB = Den[Np + n];
porosity = Poros[n];
perm = Perm[n];
GreyDiff = Perm[n];
perm = 1.0;
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
Kn_grey = GreyKn[n];
Kw_grey = GreyKw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -2196,7 +2249,23 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
rlx_setA = 1.f/tau;
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey;
Swn = 0.0;
}
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey);
Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero
// recompute the effective permeability
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5));
mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey;
Swn = 1.0;
}
// Get the 1D index based on regular data layout
ijk = Map[n];
// COMPUTE THE COLOR GRADIENT
@ -2772,21 +2841,28 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
nAB = 1.0/(nA+nB);
Aq[n] = 0.3333333333333333*nA;
Bq[n] = 0.3333333333333333*nB;
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
jA = nA*ux;
jB = nB*ux;
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nx;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
Aq[1*Np+n] = a1;
Bq[1*Np+n] = b1;
@ -2794,38 +2870,53 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
Bq[2*Np+n] = b2;
//...............................................
// q = 2
// Cq = {0,1,0}
jA = nA*uy;
jB = nB*uy;
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*ny;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
Aq[3*Np+n] = a1;
Bq[3*Np+n] = b1;
Aq[4*Np+n] = a2;
Bq[4*Np+n] = b2;
//...............................................
// q = 4
// Cq = {0,0,1}
jA = nA*uz;
jB = nB*uz;
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
//delta = 0.0;
delta = 0.111111111111111*C*W*GreyDiff*nA*nB*nAB*nz;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
Aq[5*Np+n] = a1;
Bq[5*Np+n] = b1;

558
cuda/GreyscaleColor.cu
View File

@ -1450,7 +1450,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
//CP: capillary penalty
// also turn off recoloring for grey nodes
__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta,
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
@ -1479,7 +1479,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
double Fcpx,Fcpy,Fcpz;//capillary penalty force
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
/* Corey model parameters */
double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB;
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
const double mrt_V3=0.04761904761904762;
@ -1502,15 +1506,17 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
nA = Den[n];
nB = Den[Np + n];
porosity = Poros[n];
perm = Perm[n];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
porosity = Poros[n];
//perm = Perm[n];
perm = 1.0;
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
Kn_grey = GreyKn[n];
Kw_grey = GreyKw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
// local density
rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA);
// local relaxation time
@ -1518,8 +1524,24 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff);
rlx_setA = 1.f/tau;
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey;
Swn = 0.0;
}
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey);
Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero
// recompute the effective permeability
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5));
mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey;
Swn = 1.0;
}
// Get the 1D index based on regular data layout
ijk = Map[n];
// COMPUTE THE COLOR GRADIENT
@ -1585,35 +1607,35 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Greyscale Potential Gradient.....................
// Fcpx = 0.0;
// Fcpy = 0.0;
// Fcpz = 0.0;
// if (porosity!=1.0){
// //Fcpx = -3.0/18.0*(gp1-gp2+0.5*(gp7-gp8+gp9-gp10+gp11-gp12+gp13-gp14));
// //Fcpy = -3.0/18.0*(gp3-gp4+0.5*(gp7-gp8-gp9+gp10+gp15-gp16+gp17-gp18));
// //Fcpz = -3.0/18.0*(gp5-gp6+0.5*(gp11-gp12-gp13+gp14+gp15-gp16-gp17+gp18));
// Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
// Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
// Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
// Fcpx *= alpha*W/sqrt(perm);
// Fcpy *= alpha*W/sqrt(perm);
// Fcpz *= alpha*W/sqrt(perm);
// //double Fcp_mag_temp = sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
// //double Fcp_mag = Fcp_mag_temp;
// //if (Fcp_mag_temp==0.0) Fcp_mag=1.0;
// //nx = Fcpx/Fcp_mag;
// //ny = Fcpy/Fcp_mag;
// //nz = Fcpz/Fcp_mag;
// }
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
// Fcpx = 0.0;
// Fcpy = 0.0;
// Fcpz = 0.0;
// if (porosity!=1.0){
// //Fcpx = -3.0/18.0*(gp1-gp2+0.5*(gp7-gp8+gp9-gp10+gp11-gp12+gp13-gp14));
// //Fcpy = -3.0/18.0*(gp3-gp4+0.5*(gp7-gp8-gp9+gp10+gp15-gp16+gp17-gp18));
// //Fcpz = -3.0/18.0*(gp5-gp6+0.5*(gp11-gp12-gp13+gp14+gp15-gp16-gp17+gp18));
// Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
// Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
// Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
// Fcpx *= alpha*W/sqrt(perm);
// Fcpy *= alpha*W/sqrt(perm);
// Fcpz *= alpha*W/sqrt(perm);
// //double Fcp_mag_temp = sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
// //double Fcp_mag = Fcp_mag_temp;
// //if (Fcp_mag_temp==0.0) Fcp_mag=1.0;
// //nx = Fcpx/Fcp_mag;
// //ny = Fcpy/Fcp_mag;
// //nz = Fcpz/Fcp_mag;
// }
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0) Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
@ -1944,98 +1966,98 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
m17 -= fq;
m18 -= fq;
// Compute greyscale related parameters
ux = (jx/rho0+0.5*porosity*Gx+0.5*Fcpx/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uy = (jy/rho0+0.5*porosity*Gy+0.5*Fcpy/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uz = (jz/rho0+0.5*porosity*Gz+0.5*Fcpz/rho0)/(1.0+0.5*porosity*mu_eff/perm);
if (porosity==1.0){//i.e. open nodes
ux = (jx/rho0+0.5*porosity*Gx);
uy = (jy/rho0+0.5*porosity*Gy);
uz = (jz/rho0+0.5*porosity*Gz);
}
// Compute greyscale related parameters
ux = (jx/rho0+0.5*porosity*Gx+0.5*Fcpx/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uy = (jy/rho0+0.5*porosity*Gy+0.5*Fcpy/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uz = (jz/rho0+0.5*porosity*Gz+0.5*Fcpz/rho0)/(1.0+0.5*porosity*mu_eff/perm);
if (porosity==1.0){//i.e. open nodes
ux = (jx/rho0+0.5*porosity*Gx);
uy = (jy/rho0+0.5*porosity*Gy);
uz = (jz/rho0+0.5*porosity*Gz);
}
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = rho0*(-porosity*mu_eff/perm*ux + porosity*Gx)+Fcpx;
Fy = rho0*(-porosity*mu_eff/perm*uy + porosity*Gy)+Fcpy;
Fz = rho0*(-porosity*mu_eff/perm*uz + porosity*Gz)+Fcpz;
if (porosity==1.0){
Fx=rho0*(porosity*Gx);
Fy=rho0*(porosity*Gy);
Fz=rho0*(porosity*Gz);
}
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = rho0*(-porosity*mu_eff/perm*ux + porosity*Gx)+Fcpx;
Fy = rho0*(-porosity*mu_eff/perm*uy + porosity*Gy)+Fcpy;
Fz = rho0*(-porosity*mu_eff/perm*uz + porosity*Gz)+Fcpz;
if (porosity==1.0){
Fx=rho0*(porosity*Gx);
Fy=rho0*(porosity*Gy);
Fz=rho0*(porosity*Gz);
}
// write the velocity
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
//Pressure[n] = rho/3.f/porosity;
Pressure[n] = rho/3.f;
//Pressure[n] = rho/3.f/porosity;
Pressure[n] = rho/3.f;
//........................................................................
//..............carry out relaxation process..............................
//..........Toelke, Fruediger et. al. 2006................................
//---------------- NO higher-order force -------------------------------//
//---------------- NO higher-order force -------------------------------//
if (C == 0.0) nx = ny = nz = 0.0;
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2);
jx = jx + Fx;
jx = jx + Fx;
m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
jy = jy + Fy;
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
jy = jy + Fy;
m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
jz = jz + Fz;
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
jz = jz + Fz;
m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9);
m10 = m10 + rlx_setA*( - m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11);
m12 = m12 + rlx_setA*( - m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------------------------------------------------------------//
//----------------With higher-order force ------------------------------//
//----------------With higher-order force ------------------------------//
//if (C == 0.0) nx = ny = nz = 0.0;
//m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1)
// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
//m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2)
// + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
//jx = jx + Fx;
// + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
//jx = jx + Fx;
//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
//jy = jy + Fy;
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
//jy = jy + Fy;
//m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
//jz = jz + Fz;
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
//jz = jz + Fz;
//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
//m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9)
// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
////m10 = m10 + rlx_setA*( - m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
// + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
// + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
//m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11)
// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
////m12 = m12 + rlx_setA*( - m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
//m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
//m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
//m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
// + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
// + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
//m16 = m16 + rlx_setB*( - m16);
//m17 = m17 + rlx_setB*( - m17);
//m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------------------------------------------------------------//
//.................inverse transformation......................................................
// q=0
@ -2162,21 +2184,27 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
nAB = 1.0/(nA+nB);
Aq[n] = 0.3333333333333333*nA;
Bq[n] = 0.3333333333333333*nB;
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
jA = nA*ux;
jB = nB*ux;
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
delta = 0.0;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
// q = 1
//nread = neighborList[n+Np];
@ -2189,17 +2217,23 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
//...............................................
// Cq = {0,1,0}
jA = nA*uy;
jB = nB*uy;
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
delta = 0.0;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
// q = 3
//nread = neighborList[n+3*Np];
@ -2213,17 +2247,24 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
//...............................................
// q = 4
// Cq = {0,0,1}
jA = nA*uz;
jB = nB*uz;
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
delta = 0.0;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
// q = 5
//nread = neighborList[n+5*Np];
@ -2241,7 +2282,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int
//CP: capillary penalty
// also turn off recoloring for grey nodes
__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm, double *Velocity, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
int ijk,nn,n;
@ -2265,6 +2306,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
double Fcpx,Fcpy,Fcpz;//capillary penalty force
double W;//greyscale wetting strength
double Sn_grey,Sw_grey;
/* Corey model parameters */
double Kn_grey,Kw_grey;
double Swn,Krn_grey,Krw_grey,mobility_ratio,jA,jB;
const double mrt_V1=0.05263157894736842;
const double mrt_V2=0.012531328320802;
@ -2284,17 +2329,18 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
//........Get 1-D index for this thread....................
n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if (n<finish) {
// read the component number densities
nA = Den[n];
nB = Den[Np + n];
porosity = Poros[n];
perm = Perm[n];
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
porosity = Poros[n];
//perm = Perm[n];
perm = 1.0;
W = GreySolidW[n];
Sn_grey = GreySn[n];
Sw_grey = GreySw[n];
Kn_grey = GreyKn[n];
Kw_grey = GreyKw[n];
// compute phase indicator field
phi=(nA-nB)/(nA+nB);
@ -2305,8 +2351,24 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff);
rlx_setA = 1.f/tau;
rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
mu_eff = (tau_eff-0.5)/3.0;//kinematic viscosity
if (nA/(nA+nB)<Sn_grey && porosity !=1.0){
perm = Kw_grey;
Swn = 0.0;
}
else if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
Swn = (nA/(nA+nB) - Sn_grey) /(Sw_grey - Sn_grey);
Krn_grey = Kn_grey*Swn*Swn; // Corey model with exponent = 2, make sure that W cannot shift to zero
Krw_grey = Kw_grey*(1.0-Swn)*(1.0-Swn); // Corey model with exponent = 2, make sure that W cannot shift to zero
// recompute the effective permeability
perm = mu_eff*(Krn_grey*3.0/(tauA-0.5) + Krw_grey*3.0/(tauA-0.5));
mobility_ratio =(nA*Krn_grey*3.0/(tauA-0.5) - nB*Krw_grey*3.0/(tauB-0.5))/(nA*Krn_grey*3.0/(tauA-0.5) + nB*Krw_grey*3.0/(tauB-0.5));
}
else if (nA/(nA+nB)>Sw_grey && porosity !=1.0){
perm = Kn_grey;
Swn = 1.0;
}
// Get the 1D index based on regular data layout
ijk = Map[n];
// COMPUTE THE COLOR GRADIENT
@ -2372,35 +2434,35 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
//............Compute the Greyscale Potential Gradient.....................
// Fcpx = 0.0;
// Fcpy = 0.0;
// Fcpz = 0.0;
// if (porosity!=1.0){
// //Fcpx = -3.0/18.0*(gp1-gp2+0.5*(gp7-gp8+gp9-gp10+gp11-gp12+gp13-gp14));
// //Fcpy = -3.0/18.0*(gp3-gp4+0.5*(gp7-gp8-gp9+gp10+gp15-gp16+gp17-gp18));
// //Fcpz = -3.0/18.0*(gp5-gp6+0.5*(gp11-gp12-gp13+gp14+gp15-gp16-gp17+gp18));
// Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
// Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
// Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
// Fcpx *= alpha*W/sqrt(perm);
// Fcpy *= alpha*W/sqrt(perm);
// Fcpz *= alpha*W/sqrt(perm);
// double Fcp_mag_temp = sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
// double Fcp_mag = Fcp_mag_temp;
// if (Fcp_mag_temp==0.0) Fcp_mag=1.0;
// nx = Fcpx/Fcp_mag;
// ny = Fcpy/Fcp_mag;
// nz = Fcpz/Fcp_mag;
// }
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
// Fcpx = 0.0;
// Fcpy = 0.0;
// Fcpz = 0.0;
// if (porosity!=1.0){
// //Fcpx = -3.0/18.0*(gp1-gp2+0.5*(gp7-gp8+gp9-gp10+gp11-gp12+gp13-gp14));
// //Fcpy = -3.0/18.0*(gp3-gp4+0.5*(gp7-gp8-gp9+gp10+gp15-gp16+gp17-gp18));
// //Fcpz = -3.0/18.0*(gp5-gp6+0.5*(gp11-gp12-gp13+gp14+gp15-gp16-gp17+gp18));
// Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14));
// Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18));
// Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18));
// Fcpx *= alpha*W/sqrt(perm);
// Fcpy *= alpha*W/sqrt(perm);
// Fcpz *= alpha*W/sqrt(perm);
// double Fcp_mag_temp = sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
// double Fcp_mag = Fcp_mag_temp;
// if (Fcp_mag_temp==0.0) Fcp_mag=1.0;
// nx = Fcpx/Fcp_mag;
// ny = Fcpy/Fcp_mag;
// nz = Fcpz/Fcp_mag;
// }
Fcpx = nx;
Fcpy = ny;
Fcpz = nz;
double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz);
if (Fcp_mag==0.0) Fcpx=Fcpy=Fcpz=0.0;
//NOTE for open node (porosity=1.0),Fcp=0.0
Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm);
Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm);
//...........Normalize the Color Gradient.................................
C = sqrt(nx*nx+ny*ny+nz*nz);
@ -2680,98 +2742,98 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
m17 -= fq;
m18 -= fq;
// Compute greyscale related parameters
ux = (jx/rho0+0.5*porosity*Gx+0.5*Fcpx/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uy = (jy/rho0+0.5*porosity*Gy+0.5*Fcpy/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uz = (jz/rho0+0.5*porosity*Gz+0.5*Fcpz/rho0)/(1.0+0.5*porosity*mu_eff/perm);
if (porosity==1.0){//i.e. open nodes
ux = (jx/rho0+0.5*porosity*Gx);
uy = (jy/rho0+0.5*porosity*Gy);
uz = (jz/rho0+0.5*porosity*Gz);
}
// Compute greyscale related parameters
ux = (jx/rho0+0.5*porosity*Gx+0.5*Fcpx/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uy = (jy/rho0+0.5*porosity*Gy+0.5*Fcpy/rho0)/(1.0+0.5*porosity*mu_eff/perm);
uz = (jz/rho0+0.5*porosity*Gz+0.5*Fcpz/rho0)/(1.0+0.5*porosity*mu_eff/perm);
if (porosity==1.0){//i.e. open nodes
ux = (jx/rho0+0.5*porosity*Gx);
uy = (jy/rho0+0.5*porosity*Gy);
uz = (jz/rho0+0.5*porosity*Gz);
}
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = rho0*(-porosity*mu_eff/perm*ux + porosity*Gx)+Fcpx;
Fy = rho0*(-porosity*mu_eff/perm*uy + porosity*Gy)+Fcpy;
Fz = rho0*(-porosity*mu_eff/perm*uz + porosity*Gz)+Fcpz;
if (porosity==1.0){
Fx=rho0*(porosity*Gx);
Fy=rho0*(porosity*Gy);
Fz=rho0*(porosity*Gz);
}
//Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium
Fx = rho0*(-porosity*mu_eff/perm*ux + porosity*Gx)+Fcpx;
Fy = rho0*(-porosity*mu_eff/perm*uy + porosity*Gy)+Fcpy;
Fz = rho0*(-porosity*mu_eff/perm*uz + porosity*Gz)+Fcpz;
if (porosity==1.0){
Fx=rho0*(porosity*Gx);
Fy=rho0*(porosity*Gy);
Fz=rho0*(porosity*Gz);
}
// write the velocity
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
//Pressure[n] = rho/3.f/porosity;
Pressure[n] = rho/3.f;
//Pressure[n] = rho/3.f/porosity;
Pressure[n] = rho/3.f;
//........................................................................
//..............carry out relaxation process..............................
//..........Toelke, Fruediger et. al. 2006................................
//---------------- NO higher-order force -------------------------------//
//---------------- NO higher-order force -------------------------------//
if (C == 0.0) nx = ny = nz = 0.0;
m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1);
m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2);
jx = jx + Fx;
jx = jx + Fx;
m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
jy = jy + Fy;
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
jy = jy + Fy;
m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
jz = jz + Fz;
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
jz = jz + Fz;
m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
+ (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9);
m10 = m10 + rlx_setA*( - m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10);
m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11);
m12 = m12 + rlx_setA*( - m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12);
m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------------------------------------------------------------//
//----------------With higher-order force ------------------------------//
//----------------With higher-order force ------------------------------//
//if (C == 0.0) nx = ny = nz = 0.0;
//m1 = m1 + rlx_setA*((19*(ux*ux+uy*uy+uz*uz)*rho0/porosity - 11*rho) -19*alpha*C - m1)
// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity;
//m2 = m2 + rlx_setA*((3*rho - 5.5*(ux*ux+uy*uy+uz*uz)*rho0/porosity)- m2)
// + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
//jx = jx + Fx;
// + (1-0.5*rlx_setA)*11*(-Fx*ux-Fy*uy-Fz*uz)/porosity;
//jx = jx + Fx;
//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
//jy = jy + Fy;
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx);
//jy = jy + Fy;
//m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0)- m6)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
//jz = jz + Fz;
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy);
//jz = jz + Fz;
//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz);
//m9 = m9 + rlx_setA*(((2*ux*ux-uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(2*nx*nx-ny*ny-nz*nz) - m9)
// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*Fz*uz)/porosity;
////m10 = m10 + rlx_setA*( - m10);
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
// + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
//m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10)
// + (1-0.5*rlx_setA)*(-2*Fx*ux+Fy*uy+Fz*uz)/porosity;
//m11 = m11 + rlx_setA*(((uy*uy-uz*uz)*rho0/porosity) + 0.5*alpha*C*(ny*ny-nz*nz)- m11)
// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity;
////m12 = m12 + rlx_setA*( - m12);
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
//m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12)
// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity;
//m13 = m13 + rlx_setA*( (ux*uy*rho0/porosity) + 0.5*alpha*C*nx*ny - m13);
// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity;
//m14 = m14 + rlx_setA*( (uy*uz*rho0/porosity) + 0.5*alpha*C*ny*nz - m14);
// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity;
//m15 = m15 + rlx_setA*( (ux*uz*rho0/porosity) + 0.5*alpha*C*nx*nz - m15);
// + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
// + (1-0.5*rlx_setA)*(Fz*ux+Fx*uz)/porosity;
//m16 = m16 + rlx_setB*( - m16);
//m17 = m17 + rlx_setB*( - m17);
//m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------------------------------------------------------------//
//.................inverse transformation......................................................
// q=0
@ -2882,21 +2944,27 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
nAB = 1.0/(nA+nB);
Aq[n] = 0.3333333333333333*nA;
Bq[n] = 0.3333333333333333*nB;
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
jA = nA*ux;
jB = nB*ux;
delta = beta*nA*nB*nAB*0.1111111111111111*nx;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*ux))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*ux))+delta;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
delta = 0.0;
jA = 0.5*ux*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*ux*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
Aq[1*Np+n] = a1;
Bq[1*Np+n] = b1;
@ -2904,45 +2972,57 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis
Bq[2*Np+n] = b2;
//...............................................
// q = 2
// Cq = {0,1,0}
jA = nA*uy;
jB = nB*uy;
delta = beta*nA*nB*nAB*0.1111111111111111*ny;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uy))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uy))+delta;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
delta = 0.0;
jA = 0.5*uy*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uy*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
Aq[3*Np+n] = a1;
Bq[3*Np+n] = b1;
Aq[4*Np+n] = a2;
Bq[4*Np+n] = b2;
//...............................................
// q = 4
// Cq = {0,0,1}
jA = nA*uz;
jB = nB*uz;
delta = beta*nA*nB*nAB*0.1111111111111111*nz;
if (!(nA*nB*nAB>0)) delta=0;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0;
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta;
b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta;
a2 = nA*(0.1111111111111111*(1-4.5*uz))-delta;
b2 = nB*(0.1111111111111111*(1-4.5*uz))+delta;
//----------------newly added for better control of recoloring---------------//
if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0){
delta = 0.0;
jA = 0.5*uz*(nA+nB)*(1.0+mobility_ratio);
jB = 0.5*uz*(nA+nB)*(1.0-mobility_ratio);
}
if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta;
//---------------------------------------------------------------------------//
if (RecoloringOff==true && porosity !=1.0) delta=0;
a1 = (0.1111111111111111*(nA+4.5*jA))+delta;
b1 = (0.1111111111111111*(nB+4.5*jB))-delta;
a2 = (0.1111111111111111*(nA-4.5*jA))-delta;
b2 = (0.1111111111111111*(nB-4.5*jB))+delta;
Aq[5*Np+n] = a1;
Bq[5*Np+n] = b1;
Aq[6*Np+n] = a2;
Bq[6*Np+n] = b2;
//...............................................
}
}
}
@ -4530,11 +4610,11 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl
//Model-1 & 4 with capillary pressure penalty
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel, double *Pressure,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure,
dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, GreyKn, GreyKw, Poros, Perm, Vel, Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
cudaError_t err = cudaGetLastError();
@ -4546,11 +4626,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do
//Model-1 & 4 with capillary pressure penalty
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure,
double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *GreyKn, double *GreyKw, double *Poros,double *Perm,double *Vel,double *Pressure,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,
double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure,
dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<<NBLOCKS,NTHREADS >>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, GreyKn, GreyKw, Poros, Perm,Vel,Pressure,
rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np);
cudaError_t err = cudaGetLastError();

38
docs/source/userGuide/models/color/analysis/basic.rst Normal file
View File

@ -0,0 +1,38 @@
======================================
Color model -- Basic Analysis
======================================
The basic analysis routine for the LBPM color model logs a time series
of averaged information to the space-delimited CSV file ``timelog.csv``.
The ``analysis_interval`` key should be specified to control the interval at
to perform basic analysis, which corresponds to the number of simulation timesteps
to perform between analyses. The basic analysis routine is designed to
be lightweight so that it can be performed frequently, and it is almost always
useful to enable it. Results will be immediately logged
to ``timelog.csv`` which can be used to assess whether or not the simulation is
behaving as expected. Furthermore, the analysis framework will check
simulation measures to verfiy that no ``NaN`` are detected for the fluid
pressure and flow rate.
The list of measures logged to ``timelog.csv`` are defined as follows.
The region of space occupied by the wetting fluid is determined from the
phase indicator field, :math:`\Omega_w:\phi<0`
* ``sw`` -- water saturation (fluid component 2)
* ``krw`` -- water effective permeability
* ``krw`` -- effective permeability for fluid w
* ``krn`` -- effective permeability for fluid n
* ``krwf`` -- effective permeability for fluid w (with film uncertainty estimate)
* ``krnf`` -- effective permeability for fluid n (with film uncertainty estimate)
* ``vw`` -- speed for the non-wetting fluid
* ``vn`` -- speed for the wetting fluid
* ``force`` -- magnitude for effective driving force
* ``pw`` -- average pressure for fluid w
* ``pn`` -- average pressure for fluid n
* ``wet`` -- total solid wetting energy
More comprehensive analysis is performed in the ``subphase`` analysis module.

157
docs/source/userGuide/models/color/analysis/subphase.rst Normal file
View File

@ -0,0 +1,157 @@
======================================
Color model -- Subphase Analysis
======================================
The subphase analysis routine for the LBPM color model logs a time series
of averaged information to the space-delimited CSV file ``subphase.csv``.
The ``subphase_analysis_interval`` key should be specified to control the interval at
to perform basic analysis, which corresponds to the number of simulation timesteps
to perform between analyses. The subphase analys routine performs the following functions
* analyzes the connectivity of fluid phases using a connected components algorithm
* constructs iso-surfaces to represent interfaces within the system
* computes averages of physical quantities based on the defined entities
Since it is more computationally expensive to carry out these operations compared to the
basic analysis, it may be useful to choose ``subphase_analysis_interval`` to be larger than
``analysis_interval``. Nevertheless, since analysis is performed in memory, it is orders of
magnitude faster than to analyze data *in situ* rather than writing fields to disc. Monitoring
the performance in MLUPS provides an easy way to tune the analysis interval so that the
overall simulation performance is not subject to significant penalties.
There are three main entities that are constructed for subphase analysis
* :math:`\Omega_w:\phi<0` : region of space filled by fluid w
* :math:`\Omega_n:\phi<0` : region of space filled by fluid n
* :math:`\Omega_i: |\nabla \phi|<0 > \epsilon` : region of space for the interface region
The phase regions are defined identical to what is reported in ``timelog.csv``.
The interface region is defined explicitly as the portion of space where
significant composition gradients are present. For each entity, sub-entities are
constructed by running a connected components algorithm on the set. This is done to
separate the connected part of the phase from the disconnected part. This subset operation
is performed by identifying the largest connected component (based on volume)
and denoting this as the connected part of that region. The remaining portion of the
phase is lumped into a disconnected sub-entity. Subphase analysis is therefore performed
for the following six entities
* ``wc`` -- connected part of phase w
* ``wd`` -- disconnected part of phase w
* ``nc`` -- connected part of phase n
* ``nd`` -- disconnected part of phase n
* ``ic`` -- connected part of interface region
* ``id`` -- disconnected part of interface region
For each entity :math:`\Omega_k` with :math:`k\in\{wc,wd,nc,nd,ic,id\}`
an isosurface is constructed to approximate the boundary of the region,
:math:`\Gamma_k`. Once each region is identified, the following measures are determined
**Geometric invariants**
* Volume -- :math:`V_k=\int_{\Omega_k} dV`
* Surface area -- :math:`A_k=\int_{\Gamma_k} dS`
* Integral mean curvature -- :math:`H_k=\int_{\Gamma_k} \frac 12 (\kappa_1 + \kappa_2) dS`
* Euler characteristic -- :math:`\chi_k= \frac{1}{4\pi} \int_{\Gamma_k} \kappa_1 \kappa_2 dS`
**Conserved quantities**
* Total mass within the region :math:`M_k=\int_{\Omega_k} \rho dV`
* Total momentum within the region :math:`\mathbf{P}_k=\int_{\Omega_k} \rho \mathbf{u} dV`
* Total kinetic energy within the region :math:`K_k=\frac 12 \int_{\Omega_k} \rho \mathbf{u} \cdot \mathbf{u} dV`
**Thermodynamic quantities**
* Pressure -- :math:`p_k=\frac{1}{V_k}\int_{\Omega_k} p dV`
* Solid wetting energy -- :math:`\gamma_s=\int_{\Gamma_s}\gamma dS`
* Viscous dissipation -- :math:`\Phi_k=\int_{\Omega_k}\bm{\varepsilon} : \nabla \mathbf{u} dV`
The total solid wetting energy is determined by integrating the interfacial stresses in the
immediate vicinity of the solid surface :math:`\Gamma_s`. The integral of the
dissipation function is determined based on the viscous stress tensor, denoted by :math:`\bm{\varepsilon}`.
The full list of measures are associated with the labels in ``subphase.csv``
* ``time`` -- timestep
* ``rn`` -- density for phase n (input parameter)
* ``rw`` -- density for phase w (input parameter)
* ``nun`` -- kinematic viscosity for phase n (input parameter)
* ``nuw`` -- kinematic viscosity for phase w (input parameter)
* ``Fx`` -- external force in x direction (input parameter)
* ``Fy`` -- external force in y direction (input parameter)
* ``Fz`` -- external force in z direction (input parameter)
* ``iftwn`` -- interfacial tension (input parameter)
* ``wet`` -- total solid wetting energy
* ``pwc`` -- average pressure for connected part of phase w
* ``pwd`` -- average pressure for disconnected part of phase w
* ``pnc`` -- average pressure for connected part of phase n
* ``pnd`` -- average pressure for disconnected part of phase n
* ``Mwc`` -- mass for connected part of phase w
* ``Mwd`` --mass for disconnected part of phase w
* ``Mwi`` -- mass for phase within diffuse interface region
* ``Mnc`` -- mass for connected part of phase n
* ``Mnd`` -- mass for disconnected part of phase n
* ``Mni`` -- mass for phase n within diffuse interface region
* ``Msw`` -- mass for component w within 2 voxels of solid
* ``Msn`` -- mass for component n within 2 voxels of solid
* ``Pwc_x`` -- x- momentum for connected part of phase w
* ``Pwd_x`` -- x- momentum for disconnected part of phase w
* ``Pwi_x`` -- x- momentum for phase w within diffuse interface
* ``Pnc_x`` -- x- momentum for connected part of phase n
* ``Pnd_x`` -- x- momentum for disconnected part of phase n
* ``Pni_x`` -- x- momentum for phase n within diffuse interface
* ``Psw_x`` -- x- momentum for phase w within 2 voxels of solid
* ``Psn_x`` -- x- momentum for phase n within 2 voxels of solid
* ``Pwc_y`` -- y- momentum for connected part of phase w
* ``Pwd_y`` -- y- momentum for disconnected part of phase w
* ``Pwi_y`` -- y- momentum for phase w within diffuse interface
* ``Pnc_y`` -- y- momentum for connected part of phase n
* ``Pnd_y`` -- y- momentum for connected part of phase n
* ``Pni_y`` -- y- momentum for phase n within diffuse interface
* ``Psw_y`` -- y- momentum for phase w within 2 voxels of solid
* ``Psn_y`` -- y- momentum for phase n within 2 voxels of solid
* ``Pwc_z`` -- z- momentum for connected part of phase w
* ``Pwd_z`` -- z- momentum for disconnected part of phase w
* ``Pwi_z`` -- z- momentum for phase w within diffuse interface
* ``Pnc_z`` -- z- momentum for connected part of phase n
* ``Pnd_z`` -- z- momentum for disconnected part of phase n
* ``Pni_z`` -- z- momentum for phase n within diffuse interface
* ``Psw_z`` -- z- momentum for phase w within 2 voxels of solid
* ``Psn_z`` -- z- momentum for phase n within 2 voxels of solid
* ``Kwc`` -- Kinetic energy for transport within connected part of phase w
* ``Kwd`` -- Kinetic energy for transport within disconnected part of phase w
* ``Kwi`` -- Kinetic energy for transport of phase w within diffuse interface region
* ``Knc`` -- Kinetic energy for transport in connected part of phase n
* ``Knd`` -- Kinetic energy for transport within disconnected part of phase n
* ``Kni`` -- Kinetic energy for transport of phase n within diffuse interface region
* ``Dwc`` -- Viscous dissipation for conneced pathway for phase w
* ``Dwd`` -- Viscous dissipation for disconnected part of phase w
* ``Dnc`` -- Viscous dissipation for connected pathway for phase n
* ``Dnd`` -- Viscous dissipation for disconnected part of phase n
* ``Vwc`` -- Volume for connected pathway for phase w
* ``Awc`` -- Surface area for connected pathway for phase w
* ``Hwc`` -- Integral mean curvature for connected pathway for phase w
* ``Xwc`` -- Euler characteristic for connected pathway for phase w
* ``Vwd`` -- Volume for disconnected phase w
* ``Awd`` -- Surface area for disconnected phase w
* ``Hwd`` -- Integral mean curvature for disconnected phase w
* ``Xwd`` -- Euler characteristic for disconnected phase w
* ``Nwd`` -- Number of connected components in disconnected phase w
* ``Vnc`` -- Volume for connected pathway for phase n
* ``Anc`` -- Surface area for connected pathway for phase n
* ``Hnc`` -- Integral mean curvature for connected pathway for phase n
* ``Xnc`` -- Euler characteristic for connected pathway for phase n
* ``Vnd`` -- Volume for disconnected phase n
* ``And`` -- Surface area for disconnected phase n
* ``Hnd`` -- Integral mean curvature for disconnected phase n
* ``Xnd`` -- Euler characteristic for disconnected phase n
* ``Nnd`` -- number of connected components within disconnected phase n
* ``Vi`` -- volume for diffuse interface region
* ``Ai`` -- surface area for boundary of diffuse interface region
* ``Hi`` -- integral mean curvature for boundary of diffuse interface region
* ``Xi`` -- Euler characteristic for diffuse interface region
* ``Vic`` -- volume for connected interface region
* ``Aic`` -- surface area for boundary of connected interface region
* ``Hic`` -- Integral mean curvature for connected interface region
* ``Xic`` -- Euler characteristic for connected interface region
* ``Nic`` -- number of connected components in connected interface region

20
docs/source/userGuide/models/color/index.rst
View File

@ -37,6 +37,16 @@ that can be over-ridden to develop customized simulations.
protocols/*
****************************
Analysis capabilities
****************************
.. toctree::
:glob:
:maxdepth: 2
analysis/*
***************************
Model parameters
@ -44,10 +54,12 @@ Model parameters
The essential model parameters for the color model are
- :math:`\alpha` -- control the interfacial tension between fluids with key ``alpha``
- :math:`\beta` -- control the width of the interface with key ``beta``
- :math:`\tau_A` -- control the viscosity of fluid A with key ``tauA``
- :math:`\tau_B` -- control the viscosity of fluid B with key ``tauB``
- ``alpha`` -- control the interfacial tension between fluids -- :math:`0 < \alpha < 0.01`
- ``beta`` -- control the width of the interface -- :math:`\beta < 1`
- ``tauA`` -- control the viscosity of fluid A -- :math:`0.7 < \tau_A < 1.5`
- ``tauB`` -- control the viscosity of fluid B -- :math:`0.7 < \tau_B < 1.5`
- ``rhoA`` -- control the viscosity of fluid A -- :math:`0.05 < \rho_A < 1.0`
- ``rhoB`` -- control the viscosity of fluid B -- :math:`0.05 < \rho_B < 1.0`
****************************
Model Formulation

35
example/DropletCoalescence/Droplets.py Normal file
View File

@ -0,0 +1,35 @@
import numpy
import math
nx=400
ny=200
nz=200
N=nx*ny*nz
Radius=64
mesh=(nx,ny,nz)
data=numpy.ones(mesh,dtype=numpy.int8)
#print(data)
print("Create two droplets")
print("Mesh size: "+repr(mesh))
print("Droplet radius: "+repr(Radius))
gap = 6
c1x = nx/2 - gap/2 - Radius
c2x = nx/2 + gap/2 + Radius
# assign a bubble on each side
for x in range(0,200):
for y in range(0,ny):
for z in range(0,nz):
if math.sqrt((x-c1x)*(x-c1x)+(y-ny/2)*(y-ny/2)+(z-nz/2)*(z-nz/2) ) < Radius:
data[x,y,z]=2
for x in range(200,nx):
for y in range(0,ny):
for z in range(0,nz):
if math.sqrt((x-c2x)*(x-c2x)+(y-ny/2)*(y-ny/2)+(z-nz/2)*(z-nz/2) ) < Radius:
data[x,y,z]=2
data.tofile("Droplets.raw")

51
example/DropletCoalescence/input.db Normal file
View File

@ -0,0 +1,51 @@
MRT {
timestepMax = 10000
tau = 0.7
F = 1e-05, 0, 0
Restart = false
din = 1.0
dout = 1.0
flux = 0.0
}
Color {
tauA = 0.7;
tauB = 1.0;
rhoA = 1.0;
rhoB = 1.0;
alpha = 5e-3;
beta = 0.95;
F = 0, 0, 0
Restart = false
timestepMax = 40000
ComponentLabels = -2
ComponentAffinity = -0.5
}
Domain {
Filename = "Droplets.raw"
nproc = 1, 1, 2 // Number of processors (Npx,Npy,Npz)
n = 200, 200, 200 // Size of local domain (Nx,Ny,Nz)
N = 200, 200, 400 // size of the input image
voxel_length = 1.0
BC = 0 // Boundary condition type
Sw = 0.15
ReadType = "8bit"
ReadValues = -2, 1, 2
WriteValues = -2, 1, 2
ComponentLabels = -2
HistoryLabels = -2
}
Analysis {
analysis_interval = 1000 // Frequency to perform analysis
subphase_analysis_interval = 5000 // Frequency to perform analysis
restart_interval = 60000 // Frequency to write restart data
visualization_interval = 100000 // Frequency to write visualization data
restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
Visualization {
}

36
models/GreyscaleColorModel.cpp
View File

@ -311,23 +311,31 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
double *GreySolidW_host = new double [Np];
double *GreySn_host = new double [Np];
double *GreySw_host = new double [Np];
double *GreyKn_host = new double [Np];
double *GreyKw_host = new double [Np];
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
double Sn,Sw;//end-point saturation of greynodes set by users
double Kn,Kw; // endpoint effective permeability
auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
auto SnList = greyscaleColor_db->getVector<double>( "GreySnList" );
auto SwList = greyscaleColor_db->getVector<double>( "GreySwList" );
auto SnList = greyscaleColor_db->getVector<double>( "grey_endpoint_A" );
auto SwList = greyscaleColor_db->getVector<double>( "grey_endpoint_B" );
auto KnList = greyscaleColor_db->getVector<double>( "grey_endpoint_permeability_A" );
auto KwList = greyscaleColor_db->getVector<double>( "grey_endpoint_permeability_B" );
NLABELS=LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
}
if (NLABELS != SnList.size() || NLABELS != SwList.size()){
ERROR("Error: GreySolidLabels, GreySnList, and GreySwList must be the same length! \n");
ERROR("Error: GreySolidLabels, grey_endpoint_A, and grey_endpoint_B must be the same length! \n");
}
if (NLABELS != KnList.size() || NLABELS != KwList.size()){
ERROR("Error: GreySolidLabels, grey_endpoint_permeability_A, and grey_endpoint_permeability_B must be the same length! \n");
}
for (int k=0;k<Nz;k++){
@ -344,6 +352,8 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
AFFINITY=AffinityList[idx];
Sn = SnList[idx];
Sw = SwList[idx];
Kn = SnList[idx];
Kw = SwList[idx];
idx = NLABELS;
}
}
@ -352,6 +362,8 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
GreySolidW_host[idx] = AFFINITY;
GreySn_host[idx] = Sn;
GreySw_host[idx] = Sw;
GreyKn_host[idx] = Kn;
GreyKw_host[idx] = Kw;
}
}
}
@ -374,6 +386,8 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt
ScaLBL_CopyToDevice(GreySolidW, GreySolidW_host, Np*sizeof(double));
ScaLBL_CopyToDevice(GreySn, GreySn_host, Np*sizeof(double));
ScaLBL_CopyToDevice(GreySw, GreySw_host, Np*sizeof(double));
ScaLBL_CopyToDevice(GreyKn, GreySn_host, Np*sizeof(double));
ScaLBL_CopyToDevice(GreyKw, GreySw_host, Np*sizeof(double));
ScaLBL_Comm->Barrier();
delete [] GreySolidW_host;
delete [] GreySn_host;
@ -620,7 +634,9 @@ void ScaLBL_GreyscaleColorModel::Create(){
//ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySn, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySw, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreySw, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreyKn, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &GreyKw, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Porosity_dvc, sizeof(double)*Np);
ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np);
//...........................................................................
@ -944,7 +960,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//Model-1&4
@ -973,7 +989,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
//Model-1&4
@ -1006,7 +1022,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
ScaLBL_Comm_Regular->SendHalo(Phi);
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
//Model-1&4
@ -1035,7 +1051,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
//Model-1&4 with capillary pressure penalty for grey nodes
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,GreyKn,GreyKw,Porosity_dvc,Permeability_dvc,Velocity,Pressure,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
//Model-1&4

2
models/GreyscaleColorModel.h
View File

@ -70,6 +70,8 @@ public:
double *GreySolidW;
double *GreySn;
double *GreySw;
double *GreyKn;
double *GreyKw;
//double *ColorGrad;
double *Velocity;
double *Pressure;

21
tests/lbpm_color_simulator.cpp
View File

@ -34,15 +34,15 @@ int main( int argc, char **argv )
// Load the input database
auto db = std::make_shared<Database>( argv[1] );
if (argc > 2) {
SimulationMode = "development";
SimulationMode = "legacy";
}
if ( rank == 0 ) {
printf( "********************************************************\n" );
printf( "Running Color LBM \n" );
printf( "********************************************************\n" );
if (SimulationMode == "development")
printf("**** DEVELOPMENT MODE ENABLED *************\n");
if (SimulationMode == "legacy")
printf("**** LEGACY MODE ENABLED *************\n");
}
// Initialize compute device
int device = ScaLBL_SetDevice( rank );
@ -66,13 +66,16 @@ int main( int argc, char **argv )
// structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
if (SimulationMode == "development"){
if (SimulationMode == "legacy"){
ColorModel.Run();
}
else {
double MLUPS=0.0;
int timestep = 0;
bool ContinueSimulation = true;
/* Variables for simulation protocols */
auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "none" );
auto PROTOCOL = ColorModel.color_db->getWithDefault<std::string>( "protocol", "default" );
/* image sequence protocol */
int IMAGE_INDEX = 0;
int IMAGE_COUNT = 0;
@ -123,6 +126,7 @@ int main( int argc, char **argv )
else{
if (rank==0) printf("Finished simulating image sequence \n");
ColorModel.timestep = ColorModel.timestepMax;
ContinueSimulation = false;
}
}
/*********************************************************/
@ -144,7 +148,7 @@ int main( int argc, char **argv )
double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
/* stop simulation if previous point was sufficiently close to the endpoint*/
if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false;
if (ContinueSimulation){
if (ContinueSimulation && SKIP_TIMESTEPS > 0 ){
while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){
timestep += ANALYSIS_INTERVAL;
if (PROTOCOL == "fractional flow") {
@ -173,9 +177,8 @@ int main( int argc, char **argv )
/*********************************************************/
}
}
else
ColorModel.Run();
PROFILE_STOP( "Main" );
auto file = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
auto level = db->getWithDefault<int>( "TimerLevel", 1 );