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

add a full form of Guo-type body force scheme, but now using it now

Browse Source
This commit is contained in:
Rex Zhe Li 2020-08-05 17:33:54 -04:00
parent 80c7afc27c
commit 99ee51d8e1
2 changed files with 398 additions and 321 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
gpu/GreyscaleColor.cu
View File

@ -4,7 +4,7 @@
#define NBLOCKS 1024
#define NTHREADS 256
//Model-1
//Model-1 & 4
__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta,
@ -516,6 +516,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Ma
//........................................................................
//..............carry out relaxation process..............................
//..........Toelke, Fruediger et. al. 2006................................
//---------------- 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);
@ -540,6 +541,43 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Ma
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------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;
//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;
//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
// + (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;
//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
// + (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;
////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;
//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;
////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;
//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;
//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;
//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;
//m16 = m16 + rlx_setB*( - m16);
//m17 = m17 + rlx_setB*( - m17);
//m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//.................inverse transformation......................................................
// q=0
@ -727,7 +765,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, int *Ma
}
}
//Model-1
//Model-1 & 4
__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
@ -1183,6 +1221,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
//........................................................................
//..............carry out relaxation process..............................
//..........Toelke, Fruediger et. al. 2006................................
//---------------- 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);
@ -1207,6 +1246,43 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
m16 = m16 + rlx_setB*( - m16);
m17 = m17 + rlx_setB*( - m17);
m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//----------------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;
//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;
//m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0)- m4)
// + (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;
//m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0)- m8)
// + (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;
////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;
//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;
////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;
//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;
//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;
//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;
//m16 = m16 + rlx_setB*( - m16);
//m17 = m17 + rlx_setB*( - m17);
//m18 = m18 + rlx_setB*( - m18);
//----------------------------------------------------------------------//
//.................inverse transformation......................................................
// q=0
@ -2840,7 +2916,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist,
// }
//}
//Model-1
//Model-1 & 4
extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta,
@ -2859,7 +2935,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, doubl
}
//Model-1
//Model-1 & 4
extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den,
double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,
double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta,

635
models/GreyscaleColorModel.cpp
View File

@ -332,158 +332,168 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
delete [] phase;
}
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-2 & Model-3
//{
// // ONLY initialize grey nodes
// // Key input parameters:
// // 1. GreySolidLabels
// // labels for grey nodes
// // 2. GreySolidAffinity
// // affinity ranges [-1,1]
// // oil-wet > 0
// // water-wet < 0
// // neutral = 0
//
// double *GreySolidPhi_host = new double [Nx*Ny*Nz];
// //initialize grey solid phase field
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// GreySolidPhi_host[n]=0.f;
// }
// }
// }
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// size_t NLABELS=0;
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// double *Dst;
// Dst = new double [3*3*3];
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
// int index = kk*9+jj*3+ii;
// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
// }
// }
// }
// double w_face = 1.f;
// double w_edge = 1.f;
// double w_corner = 0.f;
// //local
// Dst[13] = 0.f;
// //faces
// Dst[4] = w_face;
// Dst[10] = w_face;
// Dst[12] = w_face;
// Dst[14] = w_face;
// Dst[16] = w_face;
// Dst[22] = w_face;
// // corners
// Dst[0] = w_corner;
// Dst[2] = w_corner;
// Dst[6] = w_corner;
// Dst[8] = w_corner;
// Dst[18] = w_corner;
// Dst[20] = w_corner;
// Dst[24] = w_corner;
// Dst[26] = w_corner;
// // edges
// Dst[1] = w_edge;
// Dst[3] = w_edge;
// Dst[5] = w_edge;
// Dst[7] = w_edge;
// Dst[9] = w_edge;
// Dst[11] = w_edge;
// Dst[15] = w_edge;
// Dst[17] = w_edge;
// Dst[19] = w_edge;
// Dst[21] = w_edge;
// Dst[23] = w_edge;
// Dst[25] = w_edge;
//
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
//
// int n = k*Nx*Ny+j*Nx+i;
// signed char VALUE=Mask->id[n];
// double AFFINITY=0.f;
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
//
// if (VALUE>2){//i.e. a grey node
// double neighbor_counter = 0;
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
//
// int index = kk*9+jj*3+ii;
// double weight= Dst[index];
//
// int idi=i+ii-1;
// int idj=j+jj-1;
// int idk=k+kk-1;
//
// if (idi < 0) idi=0;
// if (idj < 0) idj=0;
// if (idk < 0) idk=0;
// if (!(idi < Nx)) idi=Nx-1;
// if (!(idj < Ny)) idj=Ny-1;
// if (!(idk < Nz)) idk=Nz-1;
//
// int nn = idk*Nx*Ny + idj*Nx + idi;
// //if (Mask->id[nn] != VALUE){//Model-2:i.e. open nodes, impermeable solid nodes or any other type of greynodes
// if (Mask->id[nn] <=0){//Model-3:i.e. only impermeable solid nodes or any other type of greynodes
// neighbor_counter +=weight;
// }
// }
// }
// }
// if (neighbor_counter>0){
// GreySolidPhi_host[n] = AFFINITY;
// }
// }
// }
// }
// }
//
// if (rank==0){
// printf("Number of grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// signed char VALUE=LabelList[idx];
// double AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_DeviceBarrier();
//
// //debug
// //FILE *OUTFILE;
// //sprintf(LocalRankFilename,"GreySolidInit.%05i.raw",rank);
// //OUTFILE = fopen(LocalRankFilename,"wb");
// //fwrite(GreySolidPhi_host,8,N,OUTFILE);
// //fclose(OUTFILE);
//
// delete [] GreySolidPhi_host;
// delete [] Dst;
//}
void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
{
// ONLY initialize grey nodes
// Key input parameters:
// 1. GreySolidLabels
// labels for grey nodes
// 2. GreySolidAffinity
// affinity ranges [-1,1]
// oil-wet > 0
// water-wet < 0
// neutral = 0
double *SolidPotential_host = new double [Nx*Ny*Nz];
double *GreySolidGrad_host = new double [3*Np];
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
NLABELS=LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
}
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=id[n];
AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component
}
}
SolidPotential_host[n] = AFFINITY;
}
}
}
// Calculate grey-solid color-gradient
double *Dst;
Dst = new double [3*3*3];
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
}
}
}
double w_face = 1.f;
double w_edge = 0.5;
double w_corner = 0.f;
//local
Dst[13] = 0.f;
//faces
Dst[4] = w_face;
Dst[10] = w_face;
Dst[12] = w_face;
Dst[14] = w_face;
Dst[16] = w_face;
Dst[22] = w_face;
// corners
Dst[0] = w_corner;
Dst[2] = w_corner;
Dst[6] = w_corner;
Dst[8] = w_corner;
Dst[18] = w_corner;
Dst[20] = w_corner;
Dst[24] = w_corner;
Dst[26] = w_corner;
// edges
Dst[1] = w_edge;
Dst[3] = w_edge;
Dst[5] = w_edge;
Dst[7] = w_edge;
Dst[9] = w_edge;
Dst[11] = w_edge;
Dst[15] = w_edge;
Dst[17] = w_edge;
Dst[19] = w_edge;
Dst[21] = w_edge;
Dst[23] = w_edge;
Dst[25] = w_edge;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx < 0)){
double phi_x = 0.f;
double phi_y = 0.f;
double phi_z = 0.f;
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
double weight= Dst[index];
int idi=i+ii-1;
int idj=j+jj-1;
int idk=k+kk-1;
if (idi < 0) idi=0;
if (idj < 0) idj=0;
if (idk < 0) idk=0;
if (!(idi < Nx)) idi=Nx-1;
if (!(idj < Ny)) idj=Ny-1;
if (!(idk < Nz)) idk=Nz-1;
int nn = idk*Nx*Ny + idj*Nx + idi;
double vec_x = double(ii-1);
double vec_y = double(jj-1);
double vec_z = double(kk-1);
double GWNS=SolidPotential_host[nn];
phi_x += GWNS*weight*vec_x;
phi_y += GWNS*weight*vec_y;
phi_z += GWNS*weight*vec_z;
}
}
}
if (Averages->SDs(i,j,k)<2.0){
GreySolidGrad_host[idx+0*Np] = phi_x;
GreySolidGrad_host[idx+1*Np] = phi_y;
GreySolidGrad_host[idx+2*Np] = phi_z;
}
else{
GreySolidGrad_host[idx+0*Np] = 0.0;
GreySolidGrad_host[idx+1*Np] = 0.0;
GreySolidGrad_host[idx+2*Np] = 0.0;
}
}
}
}
}
if (rank==0){
printf("Number of Grey-solid labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
}
}
ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
ScaLBL_DeviceBarrier();
delete [] SolidPotential_host;
delete [] GreySolidGrad_host;
delete [] Dst;
}
////----------------------------------------------------------------------------------------------------------//
void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
{
@ -1007,7 +1017,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
// Halo exchange for phase field
ScaLBL_Comm_Regular->SendHalo(Phi);
if (greyMode==true){
//Model-1
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@ -1033,7 +1043,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
if (greyMode==true){
//Model-1
//Model-1&4
ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
@ -1067,7 +1077,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
ScaLBL_Comm_Regular->SendHalo(Phi);
if (greyMode==true){
//Model-1
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
@ -1093,7 +1103,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
}
if (greyMode==true){
//Model-1
//Model-1&4
ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidGrad,Porosity_dvc,Permeability_dvc,Velocity,
rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
@ -2207,168 +2217,159 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
//}
////----------------------------------------------------------------------------------------------------------//
void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
{
// ONLY initialize grey nodes
// Key input parameters:
// 1. GreySolidLabels
// labels for grey nodes
// 2. GreySolidAffinity
// affinity ranges [-1,1]
// oil-wet > 0
// water-wet < 0
// neutral = 0
double *SolidPotential_host = new double [Nx*Ny*Nz];
double *GreySolidGrad_host = new double [3*Np];
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
NLABELS=LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
}
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=id[n];
AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component
}
}
SolidPotential_host[n] = AFFINITY;
}
}
}
// Calculate grey-solid color-gradient
double *Dst;
Dst = new double [3*3*3];
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
}
}
}
double w_face = 1.f;
double w_edge = 0.5;
double w_corner = 0.f;
//local
Dst[13] = 0.f;
//faces
Dst[4] = w_face;
Dst[10] = w_face;
Dst[12] = w_face;
Dst[14] = w_face;
Dst[16] = w_face;
Dst[22] = w_face;
// corners
Dst[0] = w_corner;
Dst[2] = w_corner;
Dst[6] = w_corner;
Dst[8] = w_corner;
Dst[18] = w_corner;
Dst[20] = w_corner;
Dst[24] = w_corner;
Dst[26] = w_corner;
// edges
Dst[1] = w_edge;
Dst[3] = w_edge;
Dst[5] = w_edge;
Dst[7] = w_edge;
Dst[9] = w_edge;
Dst[11] = w_edge;
Dst[15] = w_edge;
Dst[17] = w_edge;
Dst[19] = w_edge;
Dst[21] = w_edge;
Dst[23] = w_edge;
Dst[25] = w_edge;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
if (!(idx < 0)){
double phi_x = 0.f;
double phi_y = 0.f;
double phi_z = 0.f;
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
int index = kk*9+jj*3+ii;
double weight= Dst[index];
int idi=i+ii-1;
int idj=j+jj-1;
int idk=k+kk-1;
if (idi < 0) idi=0;
if (idj < 0) idj=0;
if (idk < 0) idk=0;
if (!(idi < Nx)) idi=Nx-1;
if (!(idj < Ny)) idj=Ny-1;
if (!(idk < Nz)) idk=Nz-1;
int nn = idk*Nx*Ny + idj*Nx + idi;
double vec_x = double(ii-1);
double vec_y = double(jj-1);
double vec_z = double(kk-1);
double GWNS=SolidPotential_host[nn];
phi_x += GWNS*weight*vec_x;
phi_y += GWNS*weight*vec_y;
phi_z += GWNS*weight*vec_z;
}
}
}
if (Averages->SDs(i,j,k)<2.0){
GreySolidGrad_host[idx+0*Np] = phi_x;
GreySolidGrad_host[idx+1*Np] = phi_y;
GreySolidGrad_host[idx+2*Np] = phi_z;
}
else{
GreySolidGrad_host[idx+0*Np] = 0.0;
GreySolidGrad_host[idx+1*Np] = 0.0;
GreySolidGrad_host[idx+2*Np] = 0.0;
}
}
}
}
}
if (rank==0){
printf("Number of Grey-solid labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
}
}
ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
ScaLBL_DeviceBarrier();
delete [] SolidPotential_host;
delete [] GreySolidGrad_host;
delete [] Dst;
}
////----------------------------------------------------------------------------------------------------------//
//void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-2 & Model-3
//{
// // ONLY initialize grey nodes
// // Key input parameters:
// // 1. GreySolidLabels
// // labels for grey nodes
// // 2. GreySolidAffinity
// // affinity ranges [-1,1]
// // oil-wet > 0
// // water-wet < 0
// // neutral = 0
//
// double *GreySolidPhi_host = new double [Nx*Ny*Nz];
// //initialize grey solid phase field
// for (int k=0;k<Nz;k++){
// for (int j=0;j<Ny;j++){
// for (int i=0;i<Nx;i++){
// int n = k*Nx*Ny+j*Nx+i;
// GreySolidPhi_host[n]=0.f;
// }
// }
// }
//
// auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
// auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
//
// size_t NLABELS=0;
// NLABELS=LabelList.size();
// if (NLABELS != AffinityList.size()){
// ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
// }
//
// double *Dst;
// Dst = new double [3*3*3];
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
// int index = kk*9+jj*3+ii;
// Dst[index] = sqrt(double(ii-1)*double(ii-1) + double(jj-1)*double(jj-1)+ double(kk-1)*double(kk-1));
// }
// }
// }
// double w_face = 1.f;
// double w_edge = 1.f;
// double w_corner = 0.f;
// //local
// Dst[13] = 0.f;
// //faces
// Dst[4] = w_face;
// Dst[10] = w_face;
// Dst[12] = w_face;
// Dst[14] = w_face;
// Dst[16] = w_face;
// Dst[22] = w_face;
// // corners
// Dst[0] = w_corner;
// Dst[2] = w_corner;
// Dst[6] = w_corner;
// Dst[8] = w_corner;
// Dst[18] = w_corner;
// Dst[20] = w_corner;
// Dst[24] = w_corner;
// Dst[26] = w_corner;
// // edges
// Dst[1] = w_edge;
// Dst[3] = w_edge;
// Dst[5] = w_edge;
// Dst[7] = w_edge;
// Dst[9] = w_edge;
// Dst[11] = w_edge;
// Dst[15] = w_edge;
// Dst[17] = w_edge;
// Dst[19] = w_edge;
// Dst[21] = w_edge;
// Dst[23] = w_edge;
// Dst[25] = w_edge;
//
// for (int k=1; k<Nz-1; k++){
// for (int j=1; j<Ny-1; j++){
// for (int i=1; i<Nx-1; i++){
//
// int n = k*Nx*Ny+j*Nx+i;
// signed char VALUE=Mask->id[n];
// double AFFINITY=0.f;
// // Assign the affinity from the paired list
// for (unsigned int idx=0; idx < NLABELS; idx++){
// //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
// if (VALUE == LabelList[idx]){
// AFFINITY=AffinityList[idx];
// idx = NLABELS;
// //Mask->id[n] = 0; // set mask to zero since this is an immobile component
// }
// }
//
// if (VALUE>2){//i.e. a grey node
// double neighbor_counter = 0;
// for (int kk=0; kk<3; kk++){
// for (int jj=0; jj<3; jj++){
// for (int ii=0; ii<3; ii++){
//
// int index = kk*9+jj*3+ii;
// double weight= Dst[index];
//
// int idi=i+ii-1;
// int idj=j+jj-1;
// int idk=k+kk-1;
//
// if (idi < 0) idi=0;
// if (idj < 0) idj=0;
// if (idk < 0) idk=0;
// if (!(idi < Nx)) idi=Nx-1;
// if (!(idj < Ny)) idj=Ny-1;
// if (!(idk < Nz)) idk=Nz-1;
//
// int nn = idk*Nx*Ny + idj*Nx + idi;
// //if (Mask->id[nn] != VALUE){//Model-2:i.e. open nodes, impermeable solid nodes or any other type of greynodes
// if (Mask->id[nn] <=0){//Model-3:i.e. only impermeable solid nodes or any other type of greynodes
// neighbor_counter +=weight;
// }
// }
// }
// }
// if (neighbor_counter>0){
// GreySolidPhi_host[n] = AFFINITY;
// }
// }
// }
// }
// }
//
// if (rank==0){
// printf("Number of grey-solid labels: %lu \n",NLABELS);
// for (unsigned int idx=0; idx<NLABELS; idx++){
// signed char VALUE=LabelList[idx];
// double AFFINITY=AffinityList[idx];
// printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
// }
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_DeviceBarrier();
//
// //debug
// //FILE *OUTFILE;
// //sprintf(LocalRankFilename,"GreySolidInit.%05i.raw",rank);
// //OUTFILE = fopen(LocalRankFilename,"wb");
// //fwrite(GreySolidPhi_host,8,N,OUTFILE);
// //fclose(OUTFILE);
//
// delete [] GreySolidPhi_host;
// delete [] Dst;
//}
//--------- This is another old version of calculating greyscale-solid color-gradient modification-------//
// **not working effectively, to be deprecated