diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index e0f3879c..0feb5558 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1151,41 +1151,39 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq){ //...Packing for x face(2,8,10,12,14)................................ ScaLBL_D3Q19_Pack(2,dvcSendList_x,0,sendCount_x,sendbuf_x,Aq,N); - req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, 2*sendCount_x,rank_x,sendtag); - req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, 2*recvCount_X,rank_X,recvtag); + req1[0] = MPI_COMM_SCALBL.Isend(sendbuf_x, sendCount_x,rank_x,sendtag); + req2[0] = MPI_COMM_SCALBL.Irecv(recvbuf_X, recvCount_X,rank_X,recvtag); //...Packing for X face(1,7,9,11,13)................................ ScaLBL_D3Q19_Pack(1,dvcSendList_X,0,sendCount_X,sendbuf_X,Aq,N); - req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, 2*sendCount_X,rank_X,sendtag); - req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, 2*recvCount_x,rank_x,recvtag); + req1[1] = MPI_COMM_SCALBL.Isend(sendbuf_X, sendCount_X,rank_X,sendtag); + req2[1] = MPI_COMM_SCALBL.Irecv(recvbuf_x, recvCount_x,rank_x,recvtag); //...Packing for y face(4,8,9,16,18)................................. ScaLBL_D3Q19_Pack(4,dvcSendList_y,0,sendCount_y,sendbuf_y,Aq,N); - req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, 2*sendCount_y,rank_y,sendtag); - req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, 2*recvCount_Y,rank_Y,recvtag); + req1[2] = MPI_COMM_SCALBL.Isend(sendbuf_y, sendCount_y,rank_y,sendtag); + req2[2] = MPI_COMM_SCALBL.Irecv(recvbuf_Y, recvCount_Y,rank_Y,recvtag); //...Packing for Y face(3,7,10,15,17)................................. ScaLBL_D3Q19_Pack(3,dvcSendList_Y,0,sendCount_Y,sendbuf_Y,Aq,N); - req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, 2*sendCount_Y,rank_Y,sendtag); - req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, 2*recvCount_y,rank_y,recvtag); + req1[3] = MPI_COMM_SCALBL.Isend(sendbuf_Y, sendCount_Y,rank_Y,sendtag); + req2[3] = MPI_COMM_SCALBL.Irecv(recvbuf_y, recvCount_y,rank_y,recvtag); //...Packing for z face(6,12,13,16,17)................................ ScaLBL_D3Q19_Pack(6,dvcSendList_z,0,sendCount_z,sendbuf_z,Aq,N); - req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, 2*sendCount_z,rank_z,sendtag); - req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, 2*recvCount_Z,rank_Z,recvtag); + req1[4] = MPI_COMM_SCALBL.Isend(sendbuf_z, sendCount_z,rank_z,sendtag); + req2[4] = MPI_COMM_SCALBL.Irecv(recvbuf_Z, recvCount_Z,rank_Z,recvtag); //...Packing for Z face(5,11,14,15,18)................................ ScaLBL_D3Q19_Pack(5,dvcSendList_Z,0,sendCount_Z,sendbuf_Z,Aq,N); + req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, sendCount_Z,rank_Z,sendtag); + req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, recvCount_z,rank_z,recvtag); //................................................................................... - // Send all the distributions - req1[5] = MPI_COMM_SCALBL.Isend(sendbuf_Z, 2*sendCount_Z,rank_Z,sendtag); - req2[5] = MPI_COMM_SCALBL.Irecv(recvbuf_z, 2*recvCount_z,rank_z,recvtag); - } void ScaLBL_Communicator::RecvD3Q7AA(double *Aq){ @@ -1194,6 +1192,7 @@ void ScaLBL_Communicator::RecvD3Q7AA(double *Aq){ //................................................................................... // Wait for completion of D3Q19 communication MPI_COMM_SCALBL.waitAll(6,req1); + MPI_COMM_SCALBL.waitAll(6,req2); ScaLBL_DeviceBarrier(); //................................................................................... diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 4557adb3..849b8252 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -84,27 +84,16 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, double *dis double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double Gsc, double Gx, double Gy, double Gz, double *Poros,double *Perm, double *Velocity,double *Pressure); -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Aq, double *Bq, double *Den,double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, -// double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap); -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Aq, double *Bq, double *Den,double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, -// double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap); - -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Cq, double *Phi, double *Den,double *SolidForce, int start, int finish, int Np, +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, double Gx, double Gy, double Gz, double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap); -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Cq, double *Phi, double *Den,double *SolidForce, int start, int finish, int Np, +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, double Gx, double Gy, double Gz, double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap); -//extern "C" void ScaLBL_D3Q7_GreyColorIMRT_Init(double *Den, double *Aq, double *Bq, double *Phi, int start, int finish, int Np); extern "C" void ScaLBL_D3Q7_GreyColorIMRT_Init(double *Den, double *Cq, double *PhiLap, double gamma, double kappaA, double kappaB, double lambdaA, double lambdaB, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_GreyColorIMRT_Init(double *dist, double *Den, double rhoA, double rhoB, int Np); @@ -113,9 +102,9 @@ extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(int *NeighborList, doubl extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(double *Aq, double *Bq, double *Den, double *Phi, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(int *NeighborList, double *Cq, double *Den, double *Phi, int start, int finish, int Np); +extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(int *NeighborList, double *Cq, double *Phi, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(double *Cq, double *Den, double *Phi, int start, int finish, int Np); +extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(double *Cq, double *Phi, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_GreyscaleColor_Gradient(int *neighborList, double *Den, double *DenGrad, int start, int finish, int Np); @@ -124,7 +113,7 @@ extern "C" void ScaLBL_D3Q19_GreyscaleColor_Laplacian(int *neighborList, double extern "C" void ScaLBL_D3Q19_GreyscaleColor_Pressure(double *dist, double *Den, double *Porosity,double *Velocity, double *Pressure, double rhoA,double rhoB, int Np); -extern "C" void ScaLBL_D3Q19_GreyscaleColor_PressureTensor(int *neighborList, double *Phi, double *PressTensor, double *PhiLap, +extern "C" void ScaLBL_D3Q19_GreyscaleColor_PressureTensor(int *neighborList, double *Phi,double *Pressure, double *PressTensor, double *PhiLap, double kappaA,double kappaB,double lambdaA,double lambdaB, int start, int finish, int Np); // MRT MODEL diff --git a/gpu/GreyscaleColor.cu b/gpu/GreyscaleColor.cu index b9519cd7..4841e212 100644 --- a/gpu/GreyscaleColor.cu +++ b/gpu/GreyscaleColor.cu @@ -1329,1334 +1329,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor(int *neighborList, double } } -//__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Aq, double *Bq, double *Den,double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, -// double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap){ -// int n; -// double vx,vy,vz,v_mag; -// double ux,uy,uz,u_mag; -// double pressure;//defined for this incompressible model -// // conserved momemnts -// double jx,jy,jz; -// // non-conserved moments -// double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; -// double fq; -// // currently disable 'GeoFun' -// double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002) -// double porosity; -// double perm;//voxel permeability -// double c0, c1; //Guo's model parameters -// double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz) -// double tau,tau_eff,rlx_setA,rlx_setB; -// double mu_eff;//effective kinematic viscosity for Darcy term -// double rho0; -// double phi; -// double phi_lap;//laplacian of phase field -// double nA,nB; -// double a1,b1,a2,b2; -// double Gfs_x,Gfs_y,Gfs_z; -// double Gff_x,Gff_y,Gff_z; -// double chem_a,chem_b; -// double rlx_massA,rlx_massB; -// // *---------------------------------Pressure Tensor Gradient------------------------------------*// -// double Pxx_x,Pyy_y,Pzz_z; -// double Pxy_x,Pxy_y; -// double Pyz_y,Pyz_z; -// double Pxz_x,Pxz_z; -// double px,py,pz; //pressure gradient -// -// -// const double mrt_V1=0.05263157894736842; -// const double mrt_V2=0.012531328320802; -// const double mrt_V3=0.04761904761904762; -// const double mrt_V4=0.004594820384294068; -// const double mrt_V5=0.01587301587301587; -// const double mrt_V6=0.0555555555555555555555555; -// const double mrt_V7=0.02777777777777778; -// const double mrt_V8=0.08333333333333333; -// const double mrt_V9=0.003341687552213868; -// const double mrt_V10=0.003968253968253968; -// const double mrt_V11=0.01388888888888889; -// const double mrt_V12=0.04166666666666666; -// -// -// int S = Np/NBLOCKS/NTHREADS + 1; -// for (int s=0; s0.0)+Perm[n]*nB/(nA+nB)*int(phi<0.0); -// -// //Load pressure gradient -// px=PressureGrad[0*Np+n]; -// py=PressureGrad[1*Np+n]; -// pz=PressureGrad[2*Np+n]; -// -// //Load pressure tensor gradient -// //For reference full list of PressTensorGrad -// //PressTensorGrad[n+0*Np] = Pxx_x -// //PressTensorGrad[n+1*Np] = Pxx_y -// //PressTensorGrad[n+2*Np] = Pxx_z -// //PressTensorGrad[n+3*Np] = Pyy_x -// //PressTensorGrad[n+4*Np] = Pyy_y -// //PressTensorGrad[n+5*Np] = Pyy_z -// //PressTensorGrad[n+6*Np] = Pzz_x -// //PressTensorGrad[n+7*Np] = Pzz_y -// //PressTensorGrad[n+8*Np] = Pzz_z -// //PressTensorGrad[n+9*Np] = Pxy_x -// //PressTensorGrad[n+10*Np] = Pxy_y -// //PressTensorGrad[n+11*Np] = Pxy_z -// //PressTensorGrad[n+12*Np] = Pyz_x -// //PressTensorGrad[n+13*Np] = Pyz_y -// //PressTensorGrad[n+14*Np] = Pyz_z -// //PressTensorGrad[n+15*Np] = Pxz_x -// //PressTensorGrad[n+16*Np] = Pxz_y -// //PressTensorGrad[n+17*Np] = Pxz_z -// Pxx_x = PressTensorGrad[0*Np+n]; -// Pyy_y = PressTensorGrad[4*Np+n]; -// Pzz_z = PressTensorGrad[8*Np+n]; -// Pxy_x = PressTensorGrad[9*Np+n]; -// Pxz_x = PressTensorGrad[15*Np+n]; -// Pxy_y = PressTensorGrad[10*Np+n]; -// Pyz_y = PressTensorGrad[13*Np+n]; -// Pyz_z = PressTensorGrad[14*Np+n]; -// Pxz_z = PressTensorGrad[17*Np+n]; -// //............Compute the fluid-fluid force (gfx,gfy,gfz)................................... -// //TODO double check if you need porosity as a fre-factor -// Gff_x = porosity*px-(Pxx_x+Pxy_y+Pxz_z); -// Gff_y = porosity*py-(Pxy_x+Pyy_y+Pyz_z); -// Gff_z = porosity*pz-(Pxz_x+Pyz_y+Pzz_z); -// // fluid-solid force -// Gfs_x = (nA-nB)*SolidForce[n+0*Np]; -// Gfs_y = (nA-nB)*SolidForce[n+1*Np]; -// Gfs_z = (nA-nB)*SolidForce[n+2*Np]; -// -// // local density -// rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); -// // local relaxation time -// tau=tauA + 0.5*(1.0-phi)*(tauB-tauA); -// rlx_setA = 1.f/tau; -// rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); -// tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff); -// mu_eff = (tau_eff-0.5)/3.f;//kinematic viscosity -// -// -// //........................................................................ -// // READ THE DISTRIBUTIONS -// // (read from opposite array due to previous swap operation) -// //........................................................................ -// // q=0 -// fq = dist[n]; -// m1 = -30.0*fq; -// m2 = 12.0*fq; -// -// // q=1 -// fq = dist[2*Np+n]; -// pressure = fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jx = fq; -// m4 = -4.0*fq; -// m9 = 2.0*fq; -// m10 = -4.0*fq; -// -// // f2 = dist[10*Np+n]; -// fq = dist[1*Np+n]; -// pressure += fq; -// m1 -= 11.0*(fq); -// m2 -= 4.0*(fq); -// jx -= fq; -// m4 += 4.0*(fq); -// m9 += 2.0*(fq); -// m10 -= 4.0*(fq); -// -// // q=3 -// fq = dist[4*Np+n]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy = fq; -// m6 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 = fq; -// m12 = -2.0*fq; -// -// // q = 4 -// fq = dist[3*Np+n]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy -= fq; -// m6 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 += fq; -// m12 -= 2.0*fq; -// -// // q=5 -// fq = dist[6*Np+n]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz = fq; -// m8 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// // q = 6 -// fq = dist[5*Np+n]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz -= fq; -// m8 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// // q=7 -// fq = dist[8*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 = fq; -// m16 = fq; -// m17 = -fq; -// -// // q = 8 -// fq = dist[7*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 += fq; -// m16 -= fq; -// m17 += fq; -// -// // q=9 -// fq = dist[10*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 += fq; -// m17 += fq; -// -// // q = 10 -// fq = dist[9*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 -= fq; -// m17 -= fq; -// -// // q=11 -// fq = dist[12*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 = fq; -// m16 -= fq; -// m18 = fq; -// -// // q=12 -// fq = dist[11*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 += fq; -// m16 += fq; -// m18 -= fq; -// -// // q=13 -// fq = dist[14*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 -= fq; -// m18 -= fq; -// -// // q=14 -// fq = dist[13*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 += fq; -// m18 += fq; -// -// // q=15 -// fq = dist[16*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 = fq; -// m17 += fq; -// m18 -= fq; -// -// // q=16 -// fq = dist[15*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 += fq; -// m17 -= fq; -// m18 += fq; -// -// // q=17 -// fq = dist[18*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 += fq; -// m18 += fq; -// -// // q=18 -// fq = dist[17*Np+n]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 -= fq; -// m18 -= fq; -// //---------------------------------------------------------------------// -// -// c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm); -// if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes -// //GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity); -// c1 = porosity*0.5*GeoFun/sqrt(perm); -// if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes -// -// vx = jx/rho0+0.5*(porosity*Gx+Gff_x+Gfs_x); -// vy = jy/rho0+0.5*(porosity*Gy+Gff_y+Gfs_y); -// vz = jz/rho0+0.5*(porosity*Gz+Gff_z+Gfs_z); -// v_mag=sqrt(vx*vx+vy*vy+vz*vz); -// ux = vx/(c0+sqrt(c0*c0+c1*v_mag)); -// uy = vy/(c0+sqrt(c0*c0+c1*v_mag)); -// uz = vz/(c0+sqrt(c0*c0+c1*v_mag)); -// u_mag=sqrt(ux*ux+uy*uy+uz*uz); -// -// //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*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx + Gff_x + Gfs_x); -// Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy + Gff_y + Gfs_y); -// Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz + Gff_z + Gfs_z); -// if (porosity==1.0){ -// Fx=rho0*(Gx + Gff_x + Gfs_x); -// Fy=rho0*(Gy + Gff_y + Gfs_y); -// Fz=rho0*(Gz + Gff_z + Gfs_z); -// } -// -// //Calculate pressure for Incompressible-MRT model -// pressure=0.5/porosity*(pressure-0.5*rho0*u_mag*u_mag/porosity); -// -//// //..............carry out relaxation process............................................... -//// m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1) -//// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity; -//// m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9) -//// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*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*((rho0*(uy*uy-uz*uz)/porosity) - m11) -//// + (1-0.5*rlx_setA)*(2*Fy*uy-2*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*((rho0*ux*uy/porosity) - m13) -//// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity; -//// m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14) -//// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity; -//// m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - 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); -//// //....................................................................................................... -// -// //-------------------- IMRT collison where body force has NO higher-order terms -------------// -// //..............carry out relaxation process............................................... -// m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1); -// m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*porosity) - m2); -// 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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9); -// m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10); -// m11 = m11 + rlx_setA*((rho0*(uy*uy-uz*uz)/porosity) - m11); -// m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12); -// m13 = m13 + rlx_setA*((rho0*ux*uy/porosity) - m13); -// m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14); -// m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - m15); -// m16 = m16 + rlx_setB*( - m16); -// m17 = m17 + rlx_setB*( - m17); -// m18 = m18 + rlx_setB*( - m18); -// //....................................................................................................... -// -// //.................inverse transformation...................................................... -// // q=0 -// fq = mrt_V1*rho0-mrt_V2*m1+mrt_V3*m2; -// dist[n] = fq; -// -// // q = 1 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10); -// dist[1*Np+n] = fq; -// -// // q=2 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10); -// dist[2*Np+n] = fq; -// -// // q = 3 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// dist[3*Np+n] = fq; -// -// // q = 4 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// dist[4*Np+n] = fq; -// -// // q = 5 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// dist[5*Np+n] = fq; -// -// // q = 6 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// dist[6*Np+n] = fq; -// -// // q = 7 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17); -// dist[7*Np+n] = fq; -// -// // q = 8 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m17-m16); -// dist[8*Np+n] = fq; -// -// // q = 9 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17); -// dist[9*Np+n] = fq; -// -// // q = 10 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17); -// dist[10*Np+n] = fq; -// -// // q = 11 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m18-m16); -// dist[11*Np+n] = fq; -// -// // q = 12 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18); -// dist[12*Np+n] = fq; -// -// // q = 13 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15-0.125*(m16+m18); -// dist[13*Np+n] = fq; -// -// // q= 14 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15+0.125*(m16+m18); -// dist[14*Np+n] = fq; -// -// // q = 15 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18); -// dist[15*Np+n] = fq; -// -// // q = 16 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17); -// dist[16*Np+n] = fq; -// -// // q = 17 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18); -// dist[17*Np+n] = fq; -// -// // q = 18 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18); -// dist[18*Np+n] = fq; -// //........................................................................ -// -// //Update velocity on device -// Velocity[0*Np+n] = ux; -// Velocity[1*Np+n] = uy; -// Velocity[2*Np+n] = uz; -// //Update pressure on device -// Pressure[n] = pressure; -// -// //-----------------------Mass transport------------------------// -// // calcuale chemical potential -// chem_a = lambdaA*(nA*nA*nA-1.5*nA*nA+0.5*nA)-0.25*kappaA*phi_lap; -// chem_b = -lambdaB*(nB*nB*nB-1.5*nB*nB+0.5*nB)-0.25*kappaB*phi_lap; -// rlx_massA = 3.f-sqrt(3.f); -// rlx_massB = 3.f-sqrt(3.f); -// -// //............................................... -// // q = 0,2,4 -// // Cq = {1,0,0}, {0,1,0}, {0,0,1} -// a1 = Aq[1*Np+n]; -// b1 = Bq[1*Np+n]; -// a2 = Aq[2*Np+n]; -// b2 = Bq[2*Np+n]; -// a1 = (1.0-rlx_massA)*a1+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a+nA*ux)); -// b1 = (1.0-rlx_massB)*b1+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b+nB*ux)); -// a2 = (1.0-rlx_massA)*a2+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a-nA*ux)); -// b2 = (1.0-rlx_massB)*b2+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b-nB*ux)); -// -// Aq[1*Np+n] = a1; -// Bq[1*Np+n] = b1; -// Aq[2*Np+n] = a2; -// Bq[2*Np+n] = b2; -// -// //............................................... -// // q = 2 -// // Cq = {0,1,0} -// a1 = Aq[3*Np+n]; -// b1 = Bq[3*Np+n]; -// a2 = Aq[4*Np+n]; -// b2 = Bq[4*Np+n]; -// a1 = (1.0-rlx_massA)*a1+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a+nA*uy)); -// b1 = (1.0-rlx_massB)*b1+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b+nB*uy)); -// a2 = (1.0-rlx_massA)*a2+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a-nA*uy)); -// b2 = (1.0-rlx_massB)*b2+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b-nB*uy)); -// -// 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} -// a1 = Aq[5*Np+n]; -// b1 = Bq[5*Np+n]; -// a2 = Aq[6*Np+n]; -// b2 = Bq[6*Np+n]; -// a1 = (1.0-rlx_massA)*a1+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a+nA*uz)); -// b1 = (1.0-rlx_massB)*b1+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b+nB*uz)); -// a2 = (1.0-rlx_massA)*a2+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a-nA*uz)); -// b2 = (1.0-rlx_massB)*b2+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b-nB*uz)); -// -// Aq[5*Np+n] = a1; -// Bq[5*Np+n] = b1; -// Aq[6*Np+n] = a2; -// Bq[6*Np+n] = b2; -// //............................................... -// -// // Instantiate mass transport distributions -// // Stationary value - distribution 0 -// a1=Aq[n]; -// b1=Bq[n]; -// Aq[n] = (1.0-rlx_massA)*a1+rlx_massA*(nA-3.0*gamma*chem_a); -// Bq[n] = (1.0-rlx_massB)*b1+rlx_massB*(nB-3.0*gamma*chem_b); -// -// -// } -// } -//} - -//__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Aq, double *Bq, double *Den,double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, -// double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap){ -// -// int n, nread, nr1,nr2,nr3,nr4,nr5,nr6; -// double vx,vy,vz,v_mag; -// double ux,uy,uz,u_mag; -// double pressure;//defined for this incompressible model -// // conserved momemnts -// double jx,jy,jz; -// // non-conserved moments -// double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; -// double fq; -// // currently disable 'GeoFun' -// double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002) -// double porosity; -// double perm;//voxel permeability -// double c0, c1; //Guo's model parameters -// double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz) -// double tau,tau_eff,rlx_setA,rlx_setB; -// double mu_eff;//effective kinematic viscosity for Darcy term -// double rho0; -// double phi; -// double phi_lap;//laplacian of phase field -// double nA,nB; -// double a1,b1,a2,b2; -// double Gfs_x,Gfs_y,Gfs_z; -// double Gff_x,Gff_y,Gff_z; -// double chem_a,chem_b; -// double rlx_massA,rlx_massB; -// // *---------------------------------Pressure Tensor Gradient------------------------------------*// -// double Pxx_x,Pyy_y,Pzz_z; -// double Pxy_x,Pxy_y; -// double Pyz_y,Pyz_z; -// double Pxz_x,Pxz_z; -// double px,py,pz; //pressure gradient -// -// const double mrt_V1=0.05263157894736842; -// const double mrt_V2=0.012531328320802; -// const double mrt_V3=0.04761904761904762; -// const double mrt_V4=0.004594820384294068; -// const double mrt_V5=0.01587301587301587; -// const double mrt_V6=0.0555555555555555555555555; -// const double mrt_V7=0.02777777777777778; -// const double mrt_V8=0.08333333333333333; -// const double mrt_V9=0.003341687552213868; -// const double mrt_V10=0.003968253968253968; -// const double mrt_V11=0.01388888888888889; -// const double mrt_V12=0.04166666666666666; -// -// int S = Np/NBLOCKS/NTHREADS + 1; -// for (int s=0; s0.0)+Perm[n]*nB/(nA+nB)*int(phi<0.0); -// -// //Load pressure gradient -// px=PressureGrad[0*Np+n]; -// py=PressureGrad[1*Np+n]; -// pz=PressureGrad[2*Np+n]; -// -// //Load pressure tensor gradient -// //For reference full list of PressTensorGrad -// //PressTensorGrad[n+0*Np] = Pxx_x -// //PressTensorGrad[n+1*Np] = Pxx_y -// //PressTensorGrad[n+2*Np] = Pxx_z -// //PressTensorGrad[n+3*Np] = Pyy_x -// //PressTensorGrad[n+4*Np] = Pyy_y -// //PressTensorGrad[n+5*Np] = Pyy_z -// //PressTensorGrad[n+6*Np] = Pzz_x -// //PressTensorGrad[n+7*Np] = Pzz_y -// //PressTensorGrad[n+8*Np] = Pzz_z -// //PressTensorGrad[n+9*Np] = Pxy_x -// //PressTensorGrad[n+10*Np] = Pxy_y -// //PressTensorGrad[n+11*Np] = Pxy_z -// //PressTensorGrad[n+12*Np] = Pyz_x -// //PressTensorGrad[n+13*Np] = Pyz_y -// //PressTensorGrad[n+14*Np] = Pyz_z -// //PressTensorGrad[n+15*Np] = Pxz_x -// //PressTensorGrad[n+16*Np] = Pxz_y -// //PressTensorGrad[n+17*Np] = Pxz_z -// Pxx_x = PressTensorGrad[0*Np+n]; -// Pyy_y = PressTensorGrad[4*Np+n]; -// Pzz_z = PressTensorGrad[8*Np+n]; -// Pxy_x = PressTensorGrad[9*Np+n]; -// Pxz_x = PressTensorGrad[15*Np+n]; -// Pxy_y = PressTensorGrad[10*Np+n]; -// Pyz_y = PressTensorGrad[13*Np+n]; -// Pyz_z = PressTensorGrad[14*Np+n]; -// Pxz_z = PressTensorGrad[17*Np+n]; -// //............Compute the fluid-fluid force (gfx,gfy,gfz)................................... -// //TODO double check if you need porosity as a fre-factor -// Gff_x = porosity*px-(Pxx_x+Pxy_y+Pxz_z); -// Gff_y = porosity*py-(Pxy_x+Pyy_y+Pyz_z); -// Gff_z = porosity*pz-(Pxz_x+Pyz_y+Pzz_z); -// // fluid-solid force -// Gfs_x = (nA-nB)*SolidForce[n+0*Np]; -// Gfs_y = (nA-nB)*SolidForce[n+1*Np]; -// Gfs_z = (nA-nB)*SolidForce[n+2*Np]; -// -// // local density -// rho0=rhoA + 0.5*(1.0-phi)*(rhoB-rhoA); -// // local relaxation time -// tau=tauA + 0.5*(1.0-phi)*(tauB-tauA); -// rlx_setA = 1.f/tau; -// rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); -// tau_eff=tauA_eff + 0.5*(1.0-phi)*(tauB_eff-tauA_eff); -// mu_eff = (tau_eff-0.5)/3.f;//kinematic viscosity -// -// //........................................................................ -// // READ THE DISTRIBUTIONS -// // (read from opposite array due to previous swap operation) -// //........................................................................ -// // q=0 -// fq = dist[n]; -// m1 = -30.0*fq; -// m2 = 12.0*fq; -// -// // q=1 -// nr1 = neighborList[n]; // neighbor 2 ( > 10Np => odd part of dist) -// fq = dist[nr1]; // reading the f1 data into register fq -// pressure = fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jx = fq; -// m4 = -4.0*fq; -// m9 = 2.0*fq; -// m10 = -4.0*fq; -// -// // q=2 -// nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) -// fq = dist[nr2]; // reading the f2 data into register fq -// pressure += fq; -// m1 -= 11.0*(fq); -// m2 -= 4.0*(fq); -// jx -= fq; -// m4 += 4.0*(fq); -// m9 += 2.0*(fq); -// m10 -= 4.0*(fq); -// -// // q=3 -// nr3 = neighborList[n+2*Np]; // neighbor 4 -// fq = dist[nr3]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy = fq; -// m6 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 = fq; -// m12 = -2.0*fq; -// -// // q = 4 -// nr4 = neighborList[n+3*Np]; // neighbor 3 -// fq = dist[nr4]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jy -= fq; -// m6 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 += fq; -// m12 -= 2.0*fq; -// -// // q=5 -// nr5 = neighborList[n+4*Np]; -// fq = dist[nr5]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz = fq; -// m8 = -4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// // q = 6 -// nr6 = neighborList[n+5*Np]; -// fq = dist[nr6]; -// pressure += fq; -// m1 -= 11.0*fq; -// m2 -= 4.0*fq; -// jz -= fq; -// m8 += 4.0*fq; -// m9 -= fq; -// m10 += 2.0*fq; -// m11 -= fq; -// m12 += 2.0*fq; -// -// // q=7 -// nread = neighborList[n+6*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 = fq; -// m16 = fq; -// m17 = -fq; -// -// // q = 8 -// nread = neighborList[n+7*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 += fq; -// m16 -= fq; -// m17 += fq; -// -// // q=9 -// nread = neighborList[n+8*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jy -= fq; -// m6 -= fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 += fq; -// m17 += fq; -// -// // q = 10 -// nread = neighborList[n+9*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jy += fq; -// m6 += fq; -// m9 += fq; -// m10 += fq; -// m11 += fq; -// m12 += fq; -// m13 -= fq; -// m16 -= fq; -// m17 -= fq; -// -// // q=11 -// nread = neighborList[n+10*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 = fq; -// m16 -= fq; -// m18 = fq; -// -// // q=12 -// nread = neighborList[n+11*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 += fq; -// m16 += fq; -// m18 -= fq; -// -// // q=13 -// nread = neighborList[n+12*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx += fq; -// m4 += fq; -// jz -= fq; -// m8 -= fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 -= fq; -// m18 -= fq; -// -// // q=14 -// nread = neighborList[n+13*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jx -= fq; -// m4 -= fq; -// jz += fq; -// m8 += fq; -// m9 += fq; -// m10 += fq; -// m11 -= fq; -// m12 -= fq; -// m15 -= fq; -// m16 += fq; -// m18 += fq; -// -// // q=15 -// nread = neighborList[n+14*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 = fq; -// m17 += fq; -// m18 -= fq; -// -// // q=16 -// nread = neighborList[n+15*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 += fq; -// m17 -= fq; -// m18 += fq; -// -// // q=17 -// nread = neighborList[n+16*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy += fq; -// m6 += fq; -// jz -= fq; -// m8 -= fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 += fq; -// m18 += fq; -// -// // q=18 -// nread = neighborList[n+17*Np]; -// fq = dist[nread]; -// pressure += fq; -// m1 += 8.0*fq; -// m2 += fq; -// jy -= fq; -// m6 -= fq; -// jz += fq; -// m8 += fq; -// m9 -= 2.0*fq; -// m10 -= 2.0*fq; -// m14 -= fq; -// m17 -= fq; -// m18 -= fq; -// //---------------------------------------------------------------------// -// -// c0 = 0.5*(1.0+porosity*0.5*mu_eff/perm); -// if (porosity==1.0) c0 = 0.5;//i.e. apparent pore nodes -// //GeoFun = 1.75/sqrt(150.0*porosity*porosity*porosity); -// c1 = porosity*0.5*GeoFun/sqrt(perm); -// if (porosity==1.0) c1 = 0.0;//i.e. apparent pore nodes -// -// vx = jx/rho0+0.5*(porosity*Gx+Gff_x+Gfs_x); -// vy = jy/rho0+0.5*(porosity*Gy+Gff_y+Gfs_y); -// vz = jz/rho0+0.5*(porosity*Gz+Gff_z+Gfs_z); -// v_mag=sqrt(vx*vx+vy*vy+vz*vz); -// ux = vx/(c0+sqrt(c0*c0+c1*v_mag)); -// uy = vy/(c0+sqrt(c0*c0+c1*v_mag)); -// uz = vz/(c0+sqrt(c0*c0+c1*v_mag)); -// u_mag=sqrt(ux*ux+uy*uy+uz*uz); -// -// //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*GeoFun/sqrt(perm)*u_mag*ux + porosity*Gx + Gff_x + Gfs_x); -// Fy = rho0*(-porosity*mu_eff/perm*uy - porosity*GeoFun/sqrt(perm)*u_mag*uy + porosity*Gy + Gff_y + Gfs_y); -// Fz = rho0*(-porosity*mu_eff/perm*uz - porosity*GeoFun/sqrt(perm)*u_mag*uz + porosity*Gz + Gff_z + Gfs_z); -// if (porosity==1.0){ -// Fx=rho0*(Gx + Gff_x + Gfs_x); -// Fy=rho0*(Gy + Gff_y + Gfs_y); -// Fz=rho0*(Gz + Gff_z + Gfs_z); -// } -// -// //Calculate pressure for Incompressible-MRT model -// pressure=0.5/porosity*(pressure-0.5*rho0*u_mag*u_mag/porosity); -// -//// //..............carry out relaxation process............................................... -//// m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1) -//// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity; -//// m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9) -//// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*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*((rho0*(uy*uy-uz*uz)/porosity) - m11) -//// + (1-0.5*rlx_setA)*(2*Fy*uy-2*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*((rho0*ux*uy/porosity) - m13) -//// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity; -//// m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14) -//// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity; -//// m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - 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); -//// //....................................................................................................... -// -// //-------------------- IMRT collison where body force has NO higher-order terms -------------// -// //..............carry out relaxation process............................................... -// m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1); -// m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*porosity) - m2); -// 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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9); -// m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10); -// m11 = m11 + rlx_setA*((rho0*(uy*uy-uz*uz)/porosity) - m11); -// m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12); -// m13 = m13 + rlx_setA*((rho0*ux*uy/porosity) - m13); -// m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14); -// m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - m15); -// m16 = m16 + rlx_setB*( - m16); -// m17 = m17 + rlx_setB*( - m17); -// m18 = m18 + rlx_setB*( - m18); -// //....................................................................................................... -// -// -// //.................inverse transformation...................................................... -// // q=0 -// fq = mrt_V1*rho0-mrt_V2*m1+mrt_V3*m2; -// dist[n] = fq; -// -// // q = 1 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10); -// //nread = neighborList[n+Np]; -// dist[nr2] = fq; -// -// // q=2 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10); -// //nread = neighborList[n]; -// dist[nr1] = fq; -// -// // q = 3 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// //nread = neighborList[n+3*Np]; -// dist[nr4] = fq; -// -// // q = 4 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); -// //nread = neighborList[n+2*Np]; -// dist[nr3] = fq; -// -// // q = 5 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// //nread = neighborList[n+5*Np]; -// dist[nr6] = fq; -// -// // q = 6 -// fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); -// //nread = neighborList[n+4*Np]; -// dist[nr5] = fq; -// -// // q = 7 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17); -// nread = neighborList[n+7*Np]; -// dist[nread] = fq; -// -// // q = 8 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m17-m16); -// nread = neighborList[n+6*Np]; -// dist[nread] = fq; -// -// // q = 9 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17); -// nread = neighborList[n+9*Np]; -// dist[nread] = fq; -// -// // q = 10 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17); -// nread = neighborList[n+8*Np]; -// dist[nread] = fq; -// -// // q = 11 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m18-m16); -// nread = neighborList[n+11*Np]; -// dist[nread] = fq; -// -// // q = 12 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18); -// nread = neighborList[n+10*Np]; -// dist[nread]= fq; -// -// // q = 13 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15-0.125*(m16+m18); -// nread = neighborList[n+13*Np]; -// dist[nread] = fq; -// -// // q= 14 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15+0.125*(m16+m18); -// nread = neighborList[n+12*Np]; -// dist[nread] = fq; -// -// // q = 15 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18); -// nread = neighborList[n+15*Np]; -// dist[nread] = fq; -// -// // q = 16 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17); -// nread = neighborList[n+14*Np]; -// dist[nread] = fq; -// -// // q = 17 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18); -// nread = neighborList[n+17*Np]; -// dist[nread] = fq; -// -// // q = 18 -// fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18); -// nread = neighborList[n+16*Np]; -// dist[nread] = fq; -// //........................................................................ -// -// //Update velocity on device -// Velocity[0*Np+n] = ux; -// Velocity[1*Np+n] = uy; -// Velocity[2*Np+n] = uz; -// //Update pressure on device -// Pressure[n] = pressure; -// -// //-----------------------Mass transport------------------------// -// // calcuale chemical potential -// chem_a = lambdaA*(nA*nA*nA-1.5*nA*nA+0.5*nA)-0.25*kappaA*phi_lap; -// chem_b = -lambdaB*(nB*nB*nB-1.5*nB*nB+0.5*nB)-0.25*kappaB*phi_lap; -// rlx_massA = 3.f-sqrt(3.f); -// rlx_massB = 3.f-sqrt(3.f); -// -// //............................................... -// // q = 0,2,4 -// // Cq = {1,0,0}, {0,1,0}, {0,0,1} -// a1 = Aq[nr2]; -// b1 = Bq[nr2]; -// a2 = Aq[nr1]; -// b2 = Bq[nr1]; -// a1 = (1.0-rlx_massA)*a1+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a+nA*ux)); -// b1 = (1.0-rlx_massB)*b1+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b+nB*ux)); -// a2 = (1.0-rlx_massA)*a2+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a-nA*ux)); -// b2 = (1.0-rlx_massB)*b2+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b-nB*ux)); -// -// // q = 1 -// //nread = neighborList[n+Np]; -// Aq[nr2] = a1; -// Bq[nr2] = b1; -// // q=2 -// //nread = neighborList[n]; -// Aq[nr1] = a2; -// Bq[nr1] = b2; -// -// //............................................... -// // Cq = {0,1,0} -// a1 = Aq[nr4]; -// b1 = Bq[nr4]; -// a2 = Aq[nr3]; -// b2 = Bq[nr3]; -// a1 = (1.0-rlx_massA)*a1+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a+nA*uy)); -// b1 = (1.0-rlx_massB)*b1+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b+nB*uy)); -// a2 = (1.0-rlx_massA)*a2+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a-nA*uy)); -// b2 = (1.0-rlx_massB)*b2+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b-nB*uy)); -// -// // q = 3 -// //nread = neighborList[n+3*Np]; -// Aq[nr4] = a1; -// Bq[nr4] = b1; -// // q = 4 -// //nread = neighborList[n+2*Np]; -// Aq[nr3] = a2; -// Bq[nr3] = b2; -// -// //............................................... -// // q = 4 -// // Cq = {0,0,1} -// a1 = Aq[nr6]; -// b1 = Bq[nr6]; -// a2 = Aq[nr5]; -// b2 = Bq[nr5]; -// a1 = (1.0-rlx_massA)*a1+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a+nA*uz)); -// b1 = (1.0-rlx_massB)*b1+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b+nB*uz)); -// a2 = (1.0-rlx_massA)*a2+rlx_massA*(0.1111111111111111*4.5*(gamma*chem_a-nA*uz)); -// b2 = (1.0-rlx_massB)*b2+rlx_massB*(0.1111111111111111*4.5*(gamma*chem_b-nB*uz)); -// -// // q = 5 -// //nread = neighborList[n+5*Np]; -// Aq[nr6] = a1; -// Bq[nr6] = b1; -// // q = 6 -// //nread = neighborList[n+4*Np]; -// Aq[nr5] = a2; -// Bq[nr5] = b2; -// //............................................... -// -// // Instantiate mass transport distributions -// // Stationary value - distribution 0 -// a1=Aq[n]; -// b1=Bq[n]; -// Aq[n] = (1.0-rlx_massA)*a1+rlx_massA*(nA-3.0*gamma*chem_a); -// Bq[n] = (1.0-rlx_massB)*b1+rlx_massB*(nB-3.0*gamma*chem_b); -// -// -// } -// } -//} - -__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Cq, double *Phi, double *Den,double *SolidForce, int start, int finish, int Np, +__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, double Gx, double Gy, double Gz, double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap){ @@ -2678,15 +1351,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz) double tau,tau_eff,rlx_setA,rlx_setB; double mu_eff;//effective kinematic viscosity for Darcy term - double rho0; + double rho,rho0; double phi; double phi_lap;//laplacian of phase field double nA,nB; - //double a1,b1,a2,b2; double Gfs_x,Gfs_y,Gfs_z; double Gff_x,Gff_y,Gff_z; double chem; - //double rlx_massA,rlx_massB; double rlx_phi; double a1,a2;//PDF of phase field // *---------------------------------Pressure Tensor Gradient------------------------------------*// @@ -2716,11 +1387,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou if ( n 10Np => odd part of dist) fq = dist[nr1]; // reading the f1 data into register fq - pressure = fq; + rho += fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jx = fq; @@ -2804,7 +1475,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=2 nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) fq = dist[nr2]; // reading the f2 data into register fq - pressure += fq; + rho += fq; m1 -= 11.0*(fq); m2 -= 4.0*(fq); jx -= fq; @@ -2815,7 +1486,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=3 nr3 = neighborList[n+2*Np]; // neighbor 4 fq = dist[nr3]; - pressure += fq; + rho += fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jy = fq; @@ -2828,7 +1499,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q = 4 nr4 = neighborList[n+3*Np]; // neighbor 3 fq = dist[nr4]; - pressure += fq; + rho += fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jy -= fq; @@ -2841,7 +1512,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=5 nr5 = neighborList[n+4*Np]; fq = dist[nr5]; - pressure += fq; + rho += fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jz = fq; @@ -2854,7 +1525,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q = 6 nr6 = neighborList[n+5*Np]; fq = dist[nr6]; - pressure += fq; + rho += fq; m1 -= 11.0*fq; m2 -= 4.0*fq; jz -= fq; @@ -2867,7 +1538,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=7 nread = neighborList[n+6*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; @@ -2885,7 +1556,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q = 8 nread = neighborList[n+7*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; @@ -2903,7 +1574,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=9 nread = neighborList[n+8*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; @@ -2921,7 +1592,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q = 10 nread = neighborList[n+9*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; @@ -2939,7 +1610,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=11 nread = neighborList[n+10*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; @@ -2957,7 +1628,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=12 nread = neighborList[n+11*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; @@ -2975,7 +1646,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=13 nread = neighborList[n+12*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx += fq; @@ -2993,7 +1664,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=14 nread = neighborList[n+13*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jx -= fq; @@ -3011,7 +1682,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=15 nread = neighborList[n+14*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jy += fq; @@ -3027,7 +1698,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=16 nread = neighborList[n+15*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jy -= fq; @@ -3043,7 +1714,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=17 nread = neighborList[n+16*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jy += fq; @@ -3059,7 +1730,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // q=18 nread = neighborList[n+17*Np]; fq = dist[nread]; - pressure += fq; + rho += fq; m1 += 8.0*fq; m2 += fq; jy -= fq; @@ -3099,159 +1770,127 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou } //Calculate pressure for Incompressible-MRT model - pressure=0.5/porosity*(pressure-0.5*rho0*u_mag*u_mag/porosity); + //pressure=0.5/porosity*(pressure-0.5*rho0*u_mag*u_mag/porosity); + pressure=rho/3.0; -// //..............carry out relaxation process............................................... -// m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1) -// + (1-0.5*rlx_setA)*38*(Fx*ux+Fy*uy+Fz*uz)/porosity; -// m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*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*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9) -// + (1-0.5*rlx_setA)*(4*Fx*ux-2*Fy*uy-2*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*((rho0*(uy*uy-uz*uz)/porosity) - m11) -// + (1-0.5*rlx_setA)*(2*Fy*uy-2*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*((rho0*ux*uy/porosity) - m13) -// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity; -// m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14) -// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity; -// m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - 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); -// //....................................................................................................... - //-------------------- IMRT collison where body force has NO higher-order terms -------------// //..............carry out relaxation process............................................... - m1 = m1 + rlx_setA*((-30*rho0+19*(ux*ux+uy*uy+uz*uz)/porosity + 57*pressure*porosity) - m1); - m2 = m2 + rlx_setA*((12*rho0 - 5.5*(ux*ux+uy*uy+uz*uz)/porosity-27*pressure*porosity) - m2); + m1 = m1 + rlx_setA*((19*(jx*jx+jy*jy+jz*jz)/rho0 - 11*rho) - m1); + m2 = m2 + rlx_setA*((3*rho - 5.5*(jx*jx+jy*jy+jz*jz)/rho0)- m2); jx = jx + Fx; - m4 = m4 + rlx_setB*((-0.6666666666666666*ux*rho0) - m4) + m4 = m4 + rlx_setB*((-0.6666666666666666*jx)- m4) + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx); jy = jy + Fy; - m6 = m6 + rlx_setB*((-0.6666666666666666*uy*rho0) - m6) + m6 = m6 + rlx_setB*((-0.6666666666666666*jy)- m6) + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy); jz = jz + Fz; - m8 = m8 + rlx_setB*((-0.6666666666666666*uz*rho0) - m8) + m8 = m8 + rlx_setB*((-0.6666666666666666*jz)- m8) + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz); - m9 = m9 + rlx_setA*((rho0*(2*ux*ux-uy*uy-uz*uz)/porosity) - m9); - m10 = m10 + rlx_setA*(-0.5*rho0*((2*ux*ux-uy*uy-uz*uz)/porosity)- m10); - m11 = m11 + rlx_setA*((rho0*(uy*uy-uz*uz)/porosity) - m11); - m12 = m12 + rlx_setA*(-0.5*(rho0*(uy*uy-uz*uz)/porosity)- m12); - m13 = m13 + rlx_setA*((rho0*ux*uy/porosity) - m13); - m14 = m14 + rlx_setA*((rho0*uy*uz/porosity) - m14); - m15 = m15 + rlx_setA*((rho0*ux*uz/porosity) - m15); - m16 = m16 + rlx_setB*( - m16); - m17 = m17 + rlx_setB*( - m17); - m18 = m18 + rlx_setB*( - m18); + m9 = m9 + rlx_setA*(((2*jx*jx-jy*jy-jz*jz)/rho0) - m9); + m10 = m10 + rlx_setA*( - m10); + m11 = m11 + rlx_setA*(((jy*jy-jz*jz)/rho0) - m11); + m12 = m12 + rlx_setA*( - m12); + m13 = m13 + rlx_setA*( (jx*jy/rho0) - m13); + m14 = m14 + rlx_setA*( (jy*jz/rho0) - m14); + m15 = m15 + rlx_setA*( (jx*jz/rho0) - m15); + m16 = m16 + rlx_setB*( - m16); + m17 = m17 + rlx_setB*( - m17); + m18 = m18 + rlx_setB*( - m18); //....................................................................................................... //.................inverse transformation...................................................... // q=0 - fq = mrt_V1*rho0-mrt_V2*m1+mrt_V3*m2; + fq = mrt_V1*rho-mrt_V2*m1+mrt_V3*m2; dist[n] = fq; // q = 1 - fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10); + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jx-m4)+mrt_V6*(m9-m10); //nread = neighborList[n+Np]; dist[nr2] = fq; // q=2 - fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10); + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m4-jx)+mrt_V6*(m9-m10); //nread = neighborList[n]; dist[nr1] = fq; // q = 3 - fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jy-m6)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); //nread = neighborList[n+3*Np]; dist[nr4] = fq; // q = 4 - fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m6-jy)+mrt_V7*(m10-m9)+mrt_V8*(m11-m12); //nread = neighborList[n+2*Np]; dist[nr3] = fq; // q = 5 - fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(jz-m8)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); //nread = neighborList[n+5*Np]; dist[nr6] = fq; // q = 6 - fq = mrt_V1*rho0-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); + fq = mrt_V1*rho-mrt_V4*m1-mrt_V5*m2+0.1*(m8-jz)+mrt_V7*(m10-m9)+mrt_V8*(m12-m11); //nread = neighborList[n+4*Np]; dist[nr5] = fq; // q = 7 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jy)+0.025*(m4+m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m16-m17); nread = neighborList[n+7*Np]; dist[nread] = fq; // q = 8 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m17-m16); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jy)-0.025*(m4+m6) +mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12+0.25*m13+0.125*(m17-m16); nread = neighborList[n+6*Np]; dist[nread] = fq; // q = 9 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jy)+0.025*(m4-m6)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13+0.125*(m16+m17); nread = neighborList[n+9*Np]; dist[nread] = fq; // q = 10 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jx)+0.025*(m6-m4)+mrt_V7*m9+mrt_V11*m10+mrt_V8*m11+mrt_V12*m12-0.25*m13-0.125*(m16+m17); nread = neighborList[n+8*Np]; dist[nread] = fq; // q = 11 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m18-m16); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx+jz)+0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m18-m16); nread = neighborList[n+11*Np]; dist[nread] = fq; // q = 12 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jx+jz)-0.025*(m4+m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12+0.25*m15+0.125*(m16-m18); nread = neighborList[n+10*Np]; dist[nread]= fq; // q = 13 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15-0.125*(m16+m18); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jx-jz)+0.025*(m4-m8)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15-0.125*(m16+m18); nread = neighborList[n+13*Np]; dist[nread] = fq; // q= 14 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15+0.125*(m16+m18); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jx)+0.025*(m8-m4)+mrt_V7*m9+mrt_V11*m10-mrt_V8*m11-mrt_V12*m12-0.25*m15+0.125*(m16+m18); nread = neighborList[n+12*Np]; dist[nread] = fq; // q = 15 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy+jz)+0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m17-m18); nread = neighborList[n+15*Np]; dist[nread] = fq; // q = 16 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2-0.1*(jy+jz)-0.025*(m6+m8)-mrt_V6*m9-mrt_V7*m10+0.25*m14+0.125*(m18-m17); nread = neighborList[n+14*Np]; dist[nread] = fq; // q = 17 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jy-jz)+0.025*(m6-m8)-mrt_V6*m9-mrt_V7*m10-0.25*m14+0.125*(m17+m18); nread = neighborList[n+17*Np]; dist[nread] = fq; // q = 18 - fq = mrt_V1*rho0+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18); + fq = mrt_V1*rho+mrt_V9*m1+mrt_V10*m2+0.1*(jz-jy)+0.025*(m8-m6)-mrt_V6*m9-mrt_V7*m10-0.25*m14-0.125*(m17+m18); nread = neighborList[n+16*Np]; dist[nread] = fq; //........................................................................ @@ -3265,17 +1904,19 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou //-----------------------Mass transport------------------------// // calcuale chemical potential - chem = lambdaA*(nA*nA*nA-1.5*nA*nA+0.5*nA)-lambdaB*(nB*nB*nB-1.5*nB*nB+0.5*nB)-0.25*(kappaA+kappaB)*phi_lap; + chem = 0.125*(lambdaA+lambdaB)*(-phi+phi*phi*phi)-0.25*(kappaA+kappaB)*phi_lap; //rlx_phi = 3.f-sqrt(3.f); rlx_phi = 1.0; //............................................... // q = 0,2,4 // Cq = {1,0,0}, {0,1,0}, {0,0,1} - a1 = Cq[nr2]; - a2 = Cq[nr1]; - a1 = (1.0-rlx_phi)*a1+rlx_phi*(0.1111111111111111*4.5*(gamma*chem+phi*ux)); - a2 = (1.0-rlx_phi)*a2+rlx_phi*(0.1111111111111111*4.5*(gamma*chem-phi*ux)); + //a1 = Cq[nr2]; + //a2 = Cq[nr1]; + //a1 = (1.0-rlx_phi)*a1+rlx_phi*(0.1111111111111111*4.5*(gamma*chem+phi*ux)); + //a2 = (1.0-rlx_phi)*a2+rlx_phi*(0.1111111111111111*4.5*(gamma*chem-phi*ux)); + a1 = 0.1111111111111111*4.5*(gamma*chem+phi*ux); + a2 = 0.1111111111111111*4.5*(gamma*chem-phi*ux); // q = 1 //nread = neighborList[n+Np]; @@ -3286,10 +1927,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou //............................................... // Cq = {0,1,0} - a1 = Cq[nr4]; - a2 = Cq[nr3]; - a1 = (1.0-rlx_phi)*a1+rlx_phi*(0.1111111111111111*4.5*(gamma*chem+phi*uy)); - a2 = (1.0-rlx_phi)*a2+rlx_phi*(0.1111111111111111*4.5*(gamma*chem-phi*uy)); + //a1 = Cq[nr4]; + //a2 = Cq[nr3]; + //a1 = (1.0-rlx_phi)*a1+rlx_phi*(0.1111111111111111*4.5*(gamma*chem+phi*uy)); + //a2 = (1.0-rlx_phi)*a2+rlx_phi*(0.1111111111111111*4.5*(gamma*chem-phi*uy)); + a1 = 0.1111111111111111*4.5*(gamma*chem+phi*uy); + a2 = 0.1111111111111111*4.5*(gamma*chem-phi*uy); // q = 3 //nread = neighborList[n+3*Np]; @@ -3301,10 +1944,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou //............................................... // q = 4 // Cq = {0,0,1} - a1 = Cq[nr6]; - a2 = Cq[nr5]; - a1 = (1.0-rlx_phi)*a1+rlx_phi*(0.1111111111111111*4.5*(gamma*chem+phi*uz)); - a2 = (1.0-rlx_phi)*a2+rlx_phi*(0.1111111111111111*4.5*(gamma*chem-phi*uz)); + //a1 = Cq[nr6]; + //a2 = Cq[nr5]; + //a1 = (1.0-rlx_phi)*a1+rlx_phi*(0.1111111111111111*4.5*(gamma*chem+phi*uz)); + //a2 = (1.0-rlx_phi)*a2+rlx_phi*(0.1111111111111111*4.5*(gamma*chem-phi*uz)); + a1 = 0.1111111111111111*4.5*(gamma*chem+phi*uz); + a2 = 0.1111111111111111*4.5*(gamma*chem-phi*uz); // q = 5 //nread = neighborList[n+5*Np]; @@ -3316,14 +1961,15 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, dou // Instantiate mass transport distributions // Stationary value - distribution 0 - a1=Cq[n]; - Cq[n] = (1.0-rlx_phi)*a1+rlx_phi*(a1-3.0*gamma*chem); + //a1=Cq[n]; + //Cq[n] = (1.0-rlx_phi)*a1+rlx_phi*(phi-3.0*gamma*chem); + Cq[n] = phi-3.0*gamma*chem; } } } -__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Cq, double *Phi, double *Den,double *SolidForce, int start, int finish, int Np, +__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, double Gx, double Gy, double Gz, double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap){ @@ -3344,15 +1990,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double double Fx, Fy, Fz;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz) double tau,tau_eff,rlx_setA,rlx_setB; double mu_eff;//effective kinematic viscosity for Darcy term - double rho0; + double rho,rho0; double phi; double phi_lap;//laplacian of phase field double nA,nB; - //double a1,b1,a2,b2; double Gfs_x,Gfs_y,Gfs_z; double Gff_x,Gff_y,Gff_z; double chem; - //double rlx_massA,rlx_massB; double rlx_phi; double a1,a2;//PDF of phase field // *---------------------------------Pressure Tensor Gradient------------------------------------*// @@ -3384,12 +2028,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double if ( n>>(dist, Aq, Bq, Den, SolidForce, start, finish, Np, -// tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Gx, Gy, Gz, Poros, Perm, Velocity, Pressure,PressureGrad,PressTensorGrad,PhiLap); -// -// cudaError_t err = cudaGetLastError(); -// if (cudaSuccess != err){ -// printf("CUDA error in ScaLBL_D3Q19_AAeven_GreyscaleColorChem: %s \n",cudaGetErrorString(err)); -// } -//} -// -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Aq, double *Bq, double *Den,double *SolidForce, int start, int finish, int Np, -// double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, -// double Gx, double Gy, double Gz, -// double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap){ -// -// dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem<<>>(neighborList, dist, Aq, Bq, Den, SolidForce, start, finish, Np, -// tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Gx, Gy, Gz, -// Poros, Perm, Velocity, Pressure,PressureGrad,PressTensorGrad,PhiLap); -// -// cudaError_t err = cudaGetLastError(); -// if (cudaSuccess != err){ -// printf("CUDA error in ScaLBL_D3Q19_AAodd_GreyscaleColorChem: %s \n",cudaGetErrorString(err)); -// } -//} - -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Cq, double *Phi, double *Den,double *SolidForce, int start, int finish, int Np, +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, double Gx, double Gy, double Gz, double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleColorChem<<>>(dist, Cq, Phi, Den, SolidForce, start, finish, Np, + dvc_ScaLBL_D3Q19_AAeven_GreyscaleColorChem<<>>(dist, Cq, Phi, SolidForce, start, finish, Np, tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Gx, Gy, Gz, Poros, Perm, Velocity, Pressure,PressureGrad,PressTensorGrad,PhiLap); cudaError_t err = cudaGetLastError(); @@ -4682,12 +3232,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColorChem(double *dist, double *Cq, } } -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Cq, double *Phi, double *Den,double *SolidForce, int start, int finish, int Np, +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double *dist, double *Cq, double *Phi, double *SolidForce, int start, int finish, int Np, double tauA,double tauB,double tauA_eff,double tauB_eff,double rhoA,double rhoB,double gamma,double kappaA,double kappaB,double lambdaA,double lambdaB, double Gx, double Gy, double Gz, double *Poros,double *Perm, double *Velocity,double *Pressure,double *PressureGrad,double *PressTensorGrad,double *PhiLap){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem<<>>(neighborList, dist, Cq, Phi, Den, SolidForce, start, finish, Np, + dvc_ScaLBL_D3Q19_AAodd_GreyscaleColorChem<<>>(neighborList, dist, Cq, Phi, SolidForce, start, finish, Np, tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Gx, Gy, Gz, Poros, Perm, Velocity, Pressure,PressureGrad,PressTensorGrad,PhiLap); @@ -4697,16 +3247,8 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColorChem(int *neighborList, double } } -//extern "C" void ScaLBL_D3Q7_GreyColorIMRT_Init(double *Den, double *Aq, double *Bq, double *Phi, int start, int finish, int Np){ -// dvc_ScaLBL_D3Q7_GreyColorIMRT_Init<<>>(Den, Aq, Bq, Phi, start, finish, Np); -// cudaError_t err = cudaGetLastError(); -// if (cudaSuccess != err){ -// printf("CUDA error in ScaLBL_D3Q7_GreyColorIMRT_Init: %s \n",cudaGetErrorString(err)); -// } -//} - -extern "C" void ScaLBL_D3Q7_GreyColorIMRT_Init(double *Den, double *Cq, double *PhiLap, double gamma, double kappaA, double kappaB, double lambdaA, double lambdaB, int start, int finish, int Np){ - dvc_ScaLBL_D3Q7_GreyColorIMRT_Init<<>>(Den, Cq, PhiLap,gamma,kappaA,kappaB,lambdaA,lambdaB, start, finish, Np); +extern "C" void ScaLBL_D3Q7_GreyColorIMRT_Init(double *Phi, double *Cq, double *PhiLap, double gamma, double kappaA, double kappaB, double lambdaA, double lambdaB, int start, int finish, int Np){ + dvc_ScaLBL_D3Q7_GreyColorIMRT_Init<<>>(Phi, Cq, PhiLap,gamma,kappaA,kappaB,lambdaA,lambdaB, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q7_GreyColorIMRT_Init: %s \n",cudaGetErrorString(err)); @@ -4740,9 +3282,9 @@ extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(double *Aq, double *Bq, } } -extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(int *NeighborList, double *Cq, double *Den, double *Phi, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(int *NeighborList, double *Cq, double *Phi, int start, int finish, int Np){ - dvc_ScaLBL_D3Q7_AAodd_GreyscaleColorPhi<<>>(NeighborList, Cq, Den, Phi, start, finish, Np); + dvc_ScaLBL_D3Q7_AAodd_GreyscaleColorPhi<<>>(NeighborList, Cq, Phi, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -4750,9 +3292,9 @@ extern "C" void ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(int *NeighborList, double *C } } -extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(double *Cq, double *Den, double *Phi, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(double *Cq, double *Phi, int start, int finish, int Np){ - dvc_ScaLBL_D3Q7_AAeven_GreyscaleColorPhi<<>>(Cq, Den, Phi, start, finish, Np); + dvc_ScaLBL_D3Q7_AAeven_GreyscaleColorPhi<<>>(Cq, Phi, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q7_AAeven_GreyscaleColorPhi: %s \n",cudaGetErrorString(err)); @@ -4786,9 +3328,9 @@ extern "C" void ScaLBL_D3Q19_GreyscaleColor_Pressure(double *dist, double *Den, } } -extern "C" void ScaLBL_D3Q19_GreyscaleColor_PressureTensor(int *neighborList, double *Phi, double *PressTensor, double *PhiLap, +extern "C" void ScaLBL_D3Q19_GreyscaleColor_PressureTensor(int *neighborList, double *Phi,double *Pressure, double *PressTensor, double *PhiLap, double kappaA,double kappaB,double lambdaA,double lambdaB, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_GreyscaleColor_PressureTensor<<>>(neighborList,Phi,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,start,finish,Np); + dvc_ScaLBL_D3Q19_GreyscaleColor_PressureTensor<<>>(neighborList,Phi,Pressure,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,start,finish,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q19_GreyscaleColor_PressureTensor: %s \n",cudaGetErrorString(err)); diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp index 1680939b..1827b961 100644 --- a/models/GreyscaleColorModel.cpp +++ b/models/GreyscaleColorModel.cpp @@ -74,9 +74,6 @@ void ScaLBL_GreyscaleColorModel::ReadParams(string filename){ if (greyscaleColor_db->keyExists( "rhoB" )){ rhoB = greyscaleColor_db->getScalar( "rhoB" ); } -// if (greyscaleColor_db->keyExists( "Gsc" )){ -// Gsc = greyscaleColor_db->getScalar( "Gsc" ); -// } if (greyscaleColor_db->keyExists( "gamma" )){ gamma = greyscaleColor_db->getScalar( "gamma" ); } @@ -477,8 +474,8 @@ void ScaLBL_GreyscaleColorModel::Density_and_Phase_Init(){ ERROR("Error: GreyNodeLabels and GreyNodeSw must be the same length! \n"); } - double *Den_temp; - Den_temp=new double [2*Np]; +// double *Den_temp; +// Den_temp=new double [2*Np]; double nA=0.5;//to prevent use may forget to specify all greynodes, then must initialize something to start with, givning just zeros is too risky. double nB=0.5; @@ -513,18 +510,18 @@ void ScaLBL_GreyscaleColorModel::Density_and_Phase_Init(){ phi = nA-nB; } int idx = Map(i,j,k); - Den_temp[idx+0*Np] = nA; - Den_temp[idx+1*Np] = nB; + //Den_temp[idx+0*Np] = nA; + //Den_temp[idx+1*Np] = nB; Phi_temp[idx] = phi; } } } } //copy to device - ScaLBL_CopyToDevice(Den, Den_temp, 2*Np*sizeof(double)); + //ScaLBL_CopyToDevice(Den, Den_temp, 2*Np*sizeof(double)); ScaLBL_CopyToDevice(Phi, Phi_temp, 1*Np*sizeof(double)); ScaLBL_DeviceBarrier(); - delete [] Den_temp; + //delete [] Den_temp; delete [] Phi_temp; } @@ -566,20 +563,20 @@ void ScaLBL_GreyscaleColorModel::Create(){ //........................................................................... ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Cq, 7*sizeof(double)*Np);//phase field distribution ScaLBL_AllocateDeviceMemory((void **) &Permeability, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Porosity, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Pressure_dvc, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &PressureGrad, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); - //ScaLBL_AllocateDeviceMemory((void **) &Aq, 7*sizeof(double)*Np); - //ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*sizeof(double)*Np); - ScaLBL_AllocateDeviceMemory((void **) &Cq, 7*sizeof(double)*Np); - ScaLBL_AllocateDeviceMemory((void **) &Den, 2*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &PhiLap, sizeof(double)*Np);//laplacian of phase field ScaLBL_AllocateDeviceMemory((void **) &SolidForce, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &PressTensor, 6*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &PressTensorGrad, 18*sizeof(double)*Np); + //ScaLBL_AllocateDeviceMemory((void **) &Den, 2*sizeof(double)*Np); + //ScaLBL_AllocateDeviceMemory((void **) &Aq, 7*sizeof(double)*Np); + //ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*sizeof(double)*Np); //ScaLBL_AllocateDeviceMemory((void **) &DenGradA, 3*sizeof(double)*Np); //ScaLBL_AllocateDeviceMemory((void **) &DenGradB, 3*sizeof(double)*Np); //ScaLBL_AllocateDeviceMemory((void **) &DenLapA, sizeof(double)*Np); @@ -613,6 +610,7 @@ void ScaLBL_GreyscaleColorModel::Create(){ void ScaLBL_GreyscaleColorModel::Initialize(){ if (Restart == true){ + //TODO: Restart funtion is currently not working; need updates if (rank==0){ printf("Initializing density field and distributions from Restart! \n"); } @@ -635,10 +633,10 @@ void ScaLBL_GreyscaleColorModel::Initialize(){ //ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Aq, Bq, Phi, 0, ScaLBL_Comm->LastExterior(), Np);//initialize D3Q7 density components //ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Aq, Bq, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q19_GreyscaleColor_Laplacian(NeighborList, Phi, PhiLap, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_D3Q19_GreyscaleColor_Laplacian(NeighborList, Phi, PhiLap, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, 0, ScaLBL_Comm->LastExterior(), Np);//initialize D3Q7 density components - ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + //ScaLBL_D3Q19_GreyscaleColor_Laplacian(NeighborList, Phi, PhiLap, 0, ScaLBL_Comm->LastExterior(), Np); + //ScaLBL_D3Q19_GreyscaleColor_Laplacian(NeighborList, Phi, PhiLap, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + //ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, 0, ScaLBL_Comm->LastExterior(), Np);//initialize D3Q7 density components + //ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); //TODO need to initialize velocity field ! //this is required for calculating the pressure_dvc @@ -651,20 +649,21 @@ void ScaLBL_GreyscaleColorModel::Initialize(){ //ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Aq, Bq, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q19_GreyscaleColor_Laplacian(NeighborList, Phi, PhiLap, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_D3Q19_GreyscaleColor_Laplacian(NeighborList, Phi, PhiLap, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, 0, ScaLBL_Comm->LastExterior(), Np);//initialize D3Q7 density components - ScaLBL_D3Q7_GreyColorIMRT_Init(Den, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_GreyColorIMRT_Init(Phi, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, 0, ScaLBL_Comm->LastExterior(), Np);//initialize D3Q7 density components + ScaLBL_D3Q7_GreyColorIMRT_Init(Phi, Cq, PhiLap, gamma,kappaA,kappaB,lambdaA,lambdaB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); if (rank==0) printf ("Initializing distributions \n"); - ScaLBL_D3Q19_GreyColorIMRT_Init(fq, Den, rhoA, rhoB, Np); + //ScaLBL_D3Q19_GreyColorIMRT_Init(fq, Den, rhoA, rhoB, Np); + ScaLBL_D3Q19_Init(fq, Np); - //Velocity also needs initialization - if (rank==0) printf ("Initializing velocity field \n"); - double *vel_init; - vel_init = new double [3*Np]; - for (int i=0;i<3*Np;i++) vel_init[i]=0.0; - ScaLBL_CopyToDevice(Velocity,vel_init,3*Np*sizeof(double)); - ScaLBL_DeviceBarrier(); - delete [] vel_init; + //Velocity also needs initialization (for old incompressible momentum transport) + //if (rank==0) printf ("Initializing velocity field \n"); + //double *vel_init; + //vel_init = new double [3*Np]; + //for (int i=0;i<3*Np;i++) vel_init[i]=0.0; + //ScaLBL_CopyToDevice(Velocity,vel_init,3*Np*sizeof(double)); + //ScaLBL_DeviceBarrier(); + //delete [] vel_init; } } @@ -714,18 +713,15 @@ void ScaLBL_GreyscaleColorModel::Run(){ timestep++; // Compute the density field // Read for Aq, Bq happens in this routine (requires communication) - //ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL ScaLBL_Comm->SendD3Q7AA(Cq); //READ FROM NORMAL - //ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(NeighborList, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(NeighborList, Cq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - //ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE + ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(NeighborList, Cq, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q7AA(Cq); //WRITE INTO OPPOSITE ScaLBL_DeviceBarrier(); - //ScaLBL_D3Q7_AAodd_GreyscaleColorDensity(NeighborList, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(NeighborList, Cq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q7_AAodd_GreyscaleColorPhi(NeighborList, Cq, Phi, 0, ScaLBL_Comm->LastExterior(), Np); // Update local pressure - ScaLBL_D3Q19_GreyscaleColor_Pressure(fq, Den, Porosity, Velocity, Pressure_dvc, rhoA, rhoB, Np); + //ScaLBL_D3Q19_GreyscaleColor_Pressure(fq, Den, Porosity, Velocity, Pressure_dvc, rhoA, rhoB, Np); + ScaLBL_D3Q19_Pressure(fq, Pressure_dvc, Np); // Compute pressure gradient ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, Pressure_dvc, PressureGrad, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->SendHalo(Pressure_dvc); @@ -733,11 +729,13 @@ void ScaLBL_GreyscaleColorModel::Run(){ ScaLBL_Comm->RecvGrad(Pressure_dvc,PressureGrad); ScaLBL_DeviceBarrier(); // Compute Pressure Tensor - ScaLBL_Comm->SendHalo(Phi); - ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(),Np); - ScaLBL_Comm->RecvHalo(Phi); - ScaLBL_DeviceBarrier(); - ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,0,ScaLBL_Comm->LastExterior(),Np); + //NOTE send and recv halo causes problems - it errorneously changes Phi + //ScaLBL_Comm->SendHalo(Phi); + ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,Pressure_dvc,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(),Np); + //ScaLBL_Comm->RecvHalo(Phi); + //ScaLBL_DeviceBarrier(); + ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,Pressure_dvc,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,0,ScaLBL_Comm->LastExterior(),Np); + /* Compute gradient of the pressure tensor */ // call the recv Grad function once per tensor element // 1st tensor element @@ -771,22 +769,9 @@ void ScaLBL_GreyscaleColorModel::Run(){ ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &PressTensor[5*Np], &PressTensorGrad[15*Np], 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->RecvGrad(&PressTensor[5*Np],&PressTensorGrad[15*Np]); -// //compute Den gradients - component A -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); -// ScaLBL_Comm->SendHalo(&Den[0*Np]); -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, 0, ScaLBL_Comm->LastExterior(), Np); -// ScaLBL_Comm->RecvGrad(&Den[0*Np],DenGradA); -// //compute Den gradients - component B -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); -// ScaLBL_Comm->SendHalo(&Den[1*Np]); -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, 0, ScaLBL_Comm->LastExterior(), Np); -// ScaLBL_Comm->RecvGrad(&Den[1*Np],DenGradB); - ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL -// ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, -// tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc); - ScaLBL_D3Q19_AAodd_GreyscaleColorChem(NeighborList, fq, Cq, Phi, Den, SolidForce, + ScaLBL_D3Q19_AAodd_GreyscaleColorChem(NeighborList, fq, Cq, Phi, SolidForce, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Fx, Fy, Fz, Porosity, Permeability, Velocity, Pressure_dvc,PressureGrad,PressTensorGrad,PhiLap); @@ -798,31 +783,27 @@ void ScaLBL_GreyscaleColorModel::Run(){ // ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); // ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); // } -// ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, 0, ScaLBL_Comm->LastExterior(), Np, -// tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc); - ScaLBL_D3Q19_AAodd_GreyscaleColorChem(NeighborList, fq, Cq, Phi, Den, SolidForce, + ScaLBL_D3Q19_AAodd_GreyscaleColorChem(NeighborList, fq, Cq, Phi, SolidForce, 0, ScaLBL_Comm->LastExterior(), Np, tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Fx, Fy, Fz, Porosity, Permeability, Velocity, Pressure_dvc,PressureGrad,PressTensorGrad,PhiLap); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + // *************EVEN TIMESTEP*************// timestep++; // Compute the density field // Read for Aq, Bq happens in this routine (requires communication) - //ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL ScaLBL_Comm->SendD3Q7AA(Cq); //READ FROM NORMAL - //ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(Cq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - //ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE + ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(Cq, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q7AA(Cq); //WRITE INTO OPPOSITE ScaLBL_DeviceBarrier(); - //ScaLBL_D3Q7_AAeven_GreyscaleColorDensity(Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(Cq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q7_AAeven_GreyscaleColorPhi(Cq, Phi, 0, ScaLBL_Comm->LastExterior(), Np); // Update local pressure - ScaLBL_D3Q19_GreyscaleColor_Pressure(fq, Den, Porosity, Velocity, Pressure_dvc, rhoA, rhoB, Np); + //ScaLBL_D3Q19_GreyscaleColor_Pressure(fq, Den, Porosity, Velocity, Pressure_dvc, rhoA, rhoB, Np); + ScaLBL_D3Q19_Pressure(fq, Pressure_dvc, Np); // Compute pressure gradient ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, Pressure_dvc, PressureGrad, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->SendHalo(Pressure_dvc); @@ -830,11 +811,12 @@ void ScaLBL_GreyscaleColorModel::Run(){ ScaLBL_Comm->RecvGrad(Pressure_dvc,PressureGrad); ScaLBL_DeviceBarrier(); // Compute Pressure Tensor - ScaLBL_Comm->SendHalo(Phi); - ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(),Np); - ScaLBL_Comm->RecvHalo(Phi); - ScaLBL_DeviceBarrier(); - ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,0,ScaLBL_Comm->LastExterior(),Np); + //ScaLBL_Comm->SendHalo(Phi); + //NOTE send and recv halo causes problems - it errorneously changes Phi + ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,Pressure_dvc,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(),Np); + //ScaLBL_Comm->RecvHalo(Phi); + //ScaLBL_DeviceBarrier(); + ScaLBL_D3Q19_GreyscaleColor_PressureTensor(NeighborList,Phi,Pressure_dvc,PressTensor,PhiLap,kappaA,kappaB,lambdaA,lambdaB,0,ScaLBL_Comm->LastExterior(),Np); /* Compute gradient of the pressure tensor */ // call the recv Grad function once per tensor element // 1st tensor element @@ -868,21 +850,9 @@ void ScaLBL_GreyscaleColorModel::Run(){ ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &PressTensor[5*Np], &PressTensorGrad[15*Np], 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->RecvGrad(&PressTensor[5*Np],&PressTensorGrad[15*Np]); -// //compute Den gradients - component A -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); -// ScaLBL_Comm->SendHalo(&Den[0*Np]); -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[0*Np], DenGradA, 0, ScaLBL_Comm->LastExterior(), Np); -// ScaLBL_Comm->RecvGrad(&Den[0*Np],DenGradA); -// //compute Den gradients - component B -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); -// ScaLBL_Comm->SendHalo(&Den[1*Np]); -// ScaLBL_D3Q19_GreyscaleColor_Gradient(NeighborList, &Den[1*Np], DenGradB, 0, ScaLBL_Comm->LastExterior(), Np); -// ScaLBL_Comm->RecvGrad(&Den[1*Np],DenGradB); ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL -// ScaLBL_D3Q19_AAeven_GreyscaleColor(fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, -// tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc); - ScaLBL_D3Q19_AAeven_GreyscaleColorChem(fq, Cq, Phi, Den, SolidForce, + ScaLBL_D3Q19_AAeven_GreyscaleColorChem(fq, Cq, Phi, SolidForce, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Fx, Fy, Fz, Porosity, Permeability, Velocity, Pressure_dvc,PressureGrad,PressTensorGrad,PhiLap); @@ -894,9 +864,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ // ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); // ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); // } -// ScaLBL_D3Q19_AAeven_GreyscaleColor(fq, Aq, Bq, Den, DenGradA, DenGradB, SolidForce, 0, ScaLBL_Comm->LastExterior(), Np, -// tauA,tauB,tauA_eff,tauB_eff,rhoA,rhoB,Gsc,Fx,Fy,Fz,Porosity,Permeability,Velocity,Pressure_dvc); - ScaLBL_D3Q19_AAeven_GreyscaleColorChem(fq, Cq, Phi, Den, SolidForce, + ScaLBL_D3Q19_AAeven_GreyscaleColorChem(fq, Cq, Phi, SolidForce, 0, ScaLBL_Comm->LastExterior(), Np, tauA, tauB, tauA_eff, tauB_eff, rhoA, rhoB, gamma,kappaA,kappaB,lambdaA,lambdaB, Fx, Fy, Fz, Porosity, Permeability, Velocity, Pressure_dvc,PressureGrad,PressTensorGrad,PhiLap); @@ -1219,19 +1187,20 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){ // fwrite(PhaseField.data(),8,N,OUTFILE); // fclose(OUTFILE); // - ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); - FILE *AFILE; - sprintf(LocalRankFilename,"A.%05i.raw",rank); - AFILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,AFILE); - fclose(AFILE); - ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField); - FILE *BFILE; - sprintf(LocalRankFilename,"B.%05i.raw",rank); - BFILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,BFILE); - fclose(BFILE); +// ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); +// FILE *AFILE; +// sprintf(LocalRankFilename,"A.%05i.raw",rank); +// AFILE = fopen(LocalRankFilename,"wb"); +// fwrite(PhaseField.data(),8,N,AFILE); +// fclose(AFILE); +// +// ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField); +// FILE *BFILE; +// sprintf(LocalRankFilename,"B.%05i.raw",rank); +// BFILE = fopen(LocalRankFilename,"wb"); +// fwrite(PhaseField.data(),8,N,BFILE); +// fclose(BFILE); ScaLBL_Comm->RegularLayout(Map,Pressure_dvc,PhaseField); FILE *PFILE; @@ -1261,6 +1230,14 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){ fwrite(PhaseField.data(),8,N,VELZ_FILE); fclose(VELZ_FILE); + + ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); + FILE *PhiFILE; + sprintf(LocalRankFilename,"Phase.%05i.raw",rank); + PhiFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,PhiFILE); + fclose(PhiFILE); + // ScaLBL_Comm->RegularLayout(Map,&Porosity[0],PhaseField); // FILE *POROS_FILE; // sprintf(LocalRankFilename,"Porosity.%05i.raw",rank);