From e83e7940218aa20e56e3f41fabe56781f0873fc2 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Mon, 4 May 2020 23:36:27 -0400 Subject: [PATCH] GreyscaleSC model;BGK works but MRT doesnt;save the work --- common/ScaLBL.h | 15 +- gpu/GreyscaleSC.cu | 930 ++++++++++++++++-------------------- models/GreyscaleSCModel.cpp | 240 +++++++--- models/GreyscaleSCModel.h | 4 +- 4 files changed, 579 insertions(+), 610 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 526ab5a4..bf2bb9dc 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -118,32 +118,33 @@ extern "C" void ScaLBL_D3Q19_GreyscaleFE_PressureTensor(int *neighborList, doubl // GREYSCALE SHAN-CHEN MODEL (Two-component) -extern "C" void ScaLBL_D3Q19_GreyscaleSC_Init(double *distA, double *distB, double *Den, int Np); +extern "C" void ScaLBL_D3Q19_GreyscaleSC_Init(int *Map, double *distA, double *distB, double *DenA, double *DenB, int Np); -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(int *NeighborList, double *distA, double *distB, double *Den, int start, int finish, int Np); +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(int *NeighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(double *distA, double *distB, double *Den, int start, int finish, int Np); +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, int *Mpa, double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(int *Map,double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(int *neighborList, double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(int *neighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np); -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(int *Map, double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np); +extern "C" void ScaLBL_D3Q19_GreyscaleSC_Gradient(int *neighborList, int *Map, double *Den, double *DenGrad, int strideY, int strideZ,int start, int finish, int Np); // MRT MODEL extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx, double Fy, double Fz); diff --git a/gpu/GreyscaleSC.cu b/gpu/GreyscaleSC.cu index dea5476f..d0768f8c 100644 --- a/gpu/GreyscaleSC.cu +++ b/gpu/GreyscaleSC.cu @@ -3,11 +3,12 @@ #define NBLOCKS 1024 #define NTHREADS 256 -__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +__global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList,int *Map, double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np){ + int ijk; int n, nread; double vx,vy,vz,v_mag; double ux_A,uy_A,uz_A,ux_B,uy_B,uz_B,u_mag; @@ -16,7 +17,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double double jxA,jyA,jzA; double jxB,jyB,jzB; double rhoA,rhoB; - double rhoA_next,rhoB_next; + double nA,nB; // non-conserved moments double m1A,m2A,m4A,m6A,m8A,m9A,m10A,m11A,m12A,m13A,m14A,m15A,m16A,m17A,m18A; double m1B,m2B,m4B,m6B,m8B,m9B,m10B,m11B,m12B,m13B,m14B,m15B,m16B,m17B,m18B; @@ -32,8 +33,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double double FxA, FyA, FzA;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz) double FxB, FyB, FzB; double rlx_setA,rlx_setB; - double rhoA_gradx,rhoA_grady,rhoA_gradz; - double rhoB_gradx,rhoB_grady,rhoB_gradz; + double nA_gradx,nA_grady,nA_gradz; + double nB_gradx,nB_grady,nB_gradz; double GffA_x,GffA_y,GffA_z; double GfsA_x,GfsA_y,GfsA_z; double GffB_x,GffB_y,GffB_z; @@ -60,18 +61,19 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double if ( n even part of dist) fq = distA[nread]; // reading the f2 data into register fq - rhoA_next += fq; + rhoA += fq; m1A -= 11.0*(fq); m2A -= 4.0*(fq); jxA -= fq; @@ -109,7 +111,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=3 nread = neighborList[n+2*Np]; // neighbor 4 fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A -= 11.0*fq; m2A -= 4.0*fq; jyA = fq; @@ -122,7 +124,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 4 nread = neighborList[n+3*Np]; // neighbor 3 fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A -= 11.0*fq; m2A -= 4.0*fq; jyA -= fq; @@ -135,7 +137,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=5 nread = neighborList[n+4*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A -= 11.0*fq; m2A -= 4.0*fq; jzA = fq; @@ -149,7 +151,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 6 nread = neighborList[n+5*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A -= 11.0*fq; m2A -= 4.0*fq; jzA -= fq; @@ -162,7 +164,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=7 nread = neighborList[n+6*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA += fq; @@ -180,7 +182,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 8 nread = neighborList[n+7*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA -= fq; @@ -198,7 +200,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=9 nread = neighborList[n+8*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA += fq; @@ -216,7 +218,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 10 nread = neighborList[n+9*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA -= fq; @@ -234,7 +236,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=11 nread = neighborList[n+10*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA += fq; @@ -252,7 +254,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=12 nread = neighborList[n+11*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA -= fq; @@ -270,7 +272,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=13 nread = neighborList[n+12*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA += fq; @@ -288,7 +290,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=14 nread = neighborList[n+13*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jxA -= fq; @@ -306,7 +308,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=15 nread = neighborList[n+14*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jyA += fq; @@ -322,7 +324,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=16 nread = neighborList[n+15*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jyA -= fq; @@ -338,7 +340,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=17 nread = neighborList[n+16*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jyA += fq; @@ -354,7 +356,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=18 nread = neighborList[n+17*Np]; fq = distA[nread]; - rhoA_next += fq; + rhoA += fq; m1A += 8.0*fq; m2A += fq; jyA -= fq; @@ -375,14 +377,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double //........................................................................ // q=0 fq = distB[n]; - rhoB_next = fq; + rhoB = fq; m1B = -30.0*fq; m2B = 12.0*fq; // q=1 nread = neighborList[n]; // neighbor 2 fq = distB[nread]; // reading the f1 data into register fq - rhoB_next += fq; + rhoB += fq; m1B -= 11.0*fq; m2B -= 4.0*fq; jxB = fq; @@ -393,7 +395,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=2 nread = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) fq = distB[nread]; // reading the f2 data into register fq - rhoB_next += fq; + rhoB += fq; m1B -= 11.0*(fq); m2B -= 4.0*(fq); jxB -= fq; @@ -404,7 +406,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=3 nread = neighborList[n+2*Np]; // neighbor 4 fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B -= 11.0*fq; m2B -= 4.0*fq; jyB = fq; @@ -417,7 +419,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 4 nread = neighborList[n+3*Np]; // neighbor 3 fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B -= 11.0*fq; m2B -= 4.0*fq; jyB -= fq; @@ -430,7 +432,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=5 nread = neighborList[n+4*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B -= 11.0*fq; m2B -= 4.0*fq; jzB = fq; @@ -444,7 +446,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 6 nread = neighborList[n+5*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B -= 11.0*fq; m2B -= 4.0*fq; jzB -= fq; @@ -457,7 +459,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=7 nread = neighborList[n+6*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB += fq; @@ -475,7 +477,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 8 nread = neighborList[n+7*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB -= fq; @@ -493,7 +495,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=9 nread = neighborList[n+8*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB += fq; @@ -511,7 +513,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q = 10 nread = neighborList[n+9*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB -= fq; @@ -529,7 +531,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=11 nread = neighborList[n+10*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB += fq; @@ -547,7 +549,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=12 nread = neighborList[n+11*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB -= fq; @@ -565,7 +567,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=13 nread = neighborList[n+12*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB += fq; @@ -583,7 +585,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=14 nread = neighborList[n+13*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jxB -= fq; @@ -601,7 +603,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=15 nread = neighborList[n+14*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jyB += fq; @@ -617,7 +619,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=16 nread = neighborList[n+15*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jyB -= fq; @@ -633,7 +635,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=17 nread = neighborList[n+16*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jyB += fq; @@ -649,7 +651,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // q=18 nread = neighborList[n+17*Np]; fq = distB[nread]; - rhoB_next += fq; + rhoB += fq; m1B += 8.0*fq; m2B += fq; jyB -= fq; @@ -665,12 +667,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // Compute SC fluid-fluid interaction force - GffA_x = -Gsc*rhoB_gradx; - GffA_y = -Gsc*rhoB_grady; - GffA_z = -Gsc*rhoB_gradz; - GffB_x = -Gsc*rhoA_gradx; - GffB_y = -Gsc*rhoA_grady; - GffB_z = -Gsc*rhoA_gradz; + GffA_x = -Gsc*nB_gradx; + GffA_y = -Gsc*nB_grady; + GffA_z = -Gsc*nB_gradz; + GffB_x = -Gsc*nA_gradx; + GffB_y = -Gsc*nA_grady; + GffB_z = -Gsc*nA_gradz; // Compute SC fluid-solid force GfsA_x = SolidForceA[n+0*Np]; GfsA_y = SolidForceA[n+1*Np]; @@ -697,13 +699,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double u_mag=sqrt(ux_A*ux_A+uy_A*uy_A+uz_A*uz_A); //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium - FxA = rhoA*(-porosity*muA_eff/permA*ux_A - porosity*GeoFun/sqrt(permA)*u_mag*ux_A + porosity*Gx + GffA_x + GfsA_x); - FyA = rhoA*(-porosity*muA_eff/permA*uy_A - porosity*GeoFun/sqrt(permA)*u_mag*uy_A + porosity*Gy + GffA_y + GfsA_y); - FzA = rhoA*(-porosity*muA_eff/permA*uz_A - porosity*GeoFun/sqrt(permA)*u_mag*uz_A + porosity*Gz + GffA_z + GfsA_z); + FxA = nA*(-porosity*muA_eff/permA*ux_A - porosity*GeoFun/sqrt(permA)*u_mag*ux_A + porosity*Gx + GffA_x + GfsA_x); + FyA = nA*(-porosity*muA_eff/permA*uy_A - porosity*GeoFun/sqrt(permA)*u_mag*uy_A + porosity*Gy + GffA_y + GfsA_y); + FzA = nA*(-porosity*muA_eff/permA*uz_A - porosity*GeoFun/sqrt(permA)*u_mag*uz_A + porosity*Gz + GffA_z + GfsA_z); if (porosity==1.0){ - FxA=rhoA*(Gx + GffA_x + GfsA_x); - FyA=rhoA*(Gy + GffA_y + GfsA_y); - FzA=rhoA*(Gz + GffA_z + GfsA_z); + FxA=nA*(Gx + GffA_x + GfsA_x); + FyA=nA*(Gy + GffA_y + GfsA_y); + FzA=nA*(Gz + GffA_z + GfsA_z); } // ------------------- Fluid Component B -----------------------// // Compute greyscale related parameters @@ -723,62 +725,29 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double u_mag=sqrt(ux_B*ux_B+uy_B*uy_B+uz_B*uz_B); //Update the total force to include linear (Darcy) and nonlinear (Forchheimer) drags due to the porous medium - FxB = rhoB*(-porosity*muB_eff/permB*ux_B - porosity*GeoFun/sqrt(permB)*u_mag*ux_B + porosity*Gx + GffB_x + GfsB_x); - FyB = rhoB*(-porosity*muB_eff/permB*uy_B - porosity*GeoFun/sqrt(permB)*u_mag*uy_B + porosity*Gy + GffB_y + GfsB_y); - FzB = rhoB*(-porosity*muB_eff/permB*uz_B - porosity*GeoFun/sqrt(permB)*u_mag*uz_B + porosity*Gz + GffB_z + GfsB_z); + FxB = nB*(-porosity*muB_eff/permB*ux_B - porosity*GeoFun/sqrt(permB)*u_mag*ux_B + porosity*Gx + GffB_x + GfsB_x); + FyB = nB*(-porosity*muB_eff/permB*uy_B - porosity*GeoFun/sqrt(permB)*u_mag*uy_B + porosity*Gy + GffB_y + GfsB_y); + FzB = nB*(-porosity*muB_eff/permB*uz_B - porosity*GeoFun/sqrt(permB)*u_mag*uz_B + porosity*Gz + GffB_z + GfsB_z); if (porosity==1.0){ - FxB=rhoB*(Gx + GffB_x + GfsB_x); - FyB=rhoB*(Gy + GffB_y + GfsB_y); - FzB=rhoB*(Gz + GffB_z + GfsB_z); + FxB=nB*(Gx + GffB_x + GfsB_x); + FyB=nB*(Gy + GffB_y + GfsB_y); + FzB=nB*(Gz + GffB_z + GfsB_z); } // Calculate barycentric velocity of the fluid mixture - ux = (rhoA*ux_A+rhoB*ux_B)/(rhoA+rhoB); - uy = (rhoA*uy_A+rhoB*uy_B)/(rhoA+rhoB); - uz = (rhoA*uz_A+rhoB*uz_B)/(rhoA+rhoB); + ux = (nA*ux_A+nB*ux_B)/(nA+nB); + uy = (nA*uy_A+nB*uy_B)/(nA+nB); + uz = (nA*uz_A+nB*uz_B)/(nA+nB); -// //..............carry out relaxation process............................................... -// m1 = m1 + rlx_setA*((-30*Den+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*Den - 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*Den) - m4) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx); -// jy = jy + Fy; -// m6 = m6 + rlx_setB*((-0.6666666666666666*uy*Den) - m6) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy); -// jz = jz + Fz; -// m8 = m8 + rlx_setB*((-0.6666666666666666*uz*Den) - m8) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz); -// m9 = m9 + rlx_setA*((Den*(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*Den*((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*((Den*(uy*uy-uz*uz)/porosity) - m11) -// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity; -// m12 = m12 + rlx_setA*(-0.5*(Den*(uy*uy-uz*uz)/porosity)- m12) -// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity; -// m13 = m13 + rlx_setA*((Den*ux*uy/porosity) - m13) -// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity; -// m14 = m14 + rlx_setA*((Den*uy*uz/porosity) - m14) -// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity; -// m15 = m15 + rlx_setA*((Den*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); -// //....................................................................................................... - // ------------------- Fluid Component A -----------------------// rlx_setA = 1.0/tauA; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); //-------------------- MRT collison where body force has NO higher-order terms -------------// //..............carry out relaxation process............................................... //TODO need to incoporate porosity - m1A = m1A + rlx_setA*((19*rhoA*(ux*ux+uy*uy+uz*uz) - 11*rhoA_next) - m1A) + m1A = m1A + rlx_setA*((19*rhoA*(ux*ux+uy*uy+uz*uz) - 11*rhoA) - m1A) + (1-0.5*rlx_setA)*38*(FxA*ux+FyA*uy+FzA*uz); - m2A = m2A + rlx_setA*((3*rhoA_next - 5.5*rhoA*(ux*ux+uy*uy+uz*uz))- m2A) + m2A = m2A + rlx_setA*((3*rhoA - 5.5*rhoA*(ux*ux+uy*uy+uz*uz))- m2A) + (1-0.5*rlx_setA)*11*(-FxA*ux-FyA*uy-FzA*uz); jxA = jxA + FxA; m4A = m4A + rlx_setB*((-0.6666666666666666*ux*rhoA)- m4A) @@ -816,165 +785,109 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // ------------------- Fluid Component A -----------------------// //.................inverse transformation...................................................... // q=0 - fq = mrt_V1*rhoA_next-mrt_V2*m1A+mrt_V3*m2A; - //f0 = mrt_V1*rhoA_next-mrt_V2*m1A+mrt_V3*m2A; + fq = mrt_V1*rhoA-mrt_V2*m1A+mrt_V3*m2A; distA[n] = fq; // q = 1 - fq = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(jxA-m4A)+mrt_V6*(m9A-m10A); - //f1 = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(jxA-m4A)+mrt_V6*(m9A-m10A); + fq = mrt_V1*rhoA-mrt_V4*m1A-mrt_V5*m2A+0.1*(jxA-m4A)+mrt_V6*(m9A-m10A); nread = neighborList[n+Np]; distA[nread] = fq; // q=2 - fq = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(m4A-jxA)+mrt_V6*(m9A-m10A); - //f2 = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(m4A-jxA)+mrt_V6*(m9A-m10A); + fq = mrt_V1*rhoA-mrt_V4*m1A-mrt_V5*m2A+0.1*(m4A-jxA)+mrt_V6*(m9A-m10A); nread = neighborList[n]; distA[nread] = fq; // q = 3 - fq = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(jyA-m6A)+mrt_V7*(m10A-m9A)+mrt_V8*(m11A-m12A); - //f3 = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(jyA-m6A)+mrt_V7*(m10A-m9A)+mrt_V8*(m11A-m12A); + fq = mrt_V1*rhoA-mrt_V4*m1A-mrt_V5*m2A+0.1*(jyA-m6A)+mrt_V7*(m10A-m9A)+mrt_V8*(m11A-m12A); nread = neighborList[n+3*Np]; distA[nread] = fq; // q = 4 - fq = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(m6A-jyA)+mrt_V7*(m10A-m9A)+mrt_V8*(m11A-m12A); - //f4 = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(m6A-jyA)+mrt_V7*(m10A-m9A)+mrt_V8*(m11A-m12A); + fq = mrt_V1*rhoA-mrt_V4*m1A-mrt_V5*m2A+0.1*(m6A-jyA)+mrt_V7*(m10A-m9A)+mrt_V8*(m11A-m12A); nread = neighborList[n+2*Np]; distA[nread] = fq; // q = 5 - fq = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(jzA-m8A)+mrt_V7*(m10A-m9A)+mrt_V8*(m12A-m11A); - //f5 = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(jzA-m8A)+mrt_V7*(m10A-m9A)+mrt_V8*(m12A-m11A); + fq = mrt_V1*rhoA-mrt_V4*m1A-mrt_V5*m2A+0.1*(jzA-m8A)+mrt_V7*(m10A-m9A)+mrt_V8*(m12A-m11A); nread = neighborList[n+5*Np]; distA[nread] = fq; // q = 6 - fq = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(m8A-jzA)+mrt_V7*(m10A-m9A)+mrt_V8*(m12A-m11A); - //f6 = mrt_V1*rhoA_next-mrt_V4*m1A-mrt_V5*m2A+0.1*(m8A-jzA)+mrt_V7*(m10A-m9A)+mrt_V8*(m12A-m11A); + fq = mrt_V1*rhoA-mrt_V4*m1A-mrt_V5*m2A+0.1*(m8A-jzA)+mrt_V7*(m10A-m9A)+mrt_V8*(m12A-m11A); nread = neighborList[n+4*Np]; distA[nread] = fq; // q = 7 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA+jyA)+0.025*(m4A+m6A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A+0.25*m13A+0.125*(m16A-m17A); - //f7 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA+jyA)+0.025*(m4A+m6A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A+0.25*m13A+0.125*(m16A-m17A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA+jyA)+0.025*(m4A+m6A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A+0.25*m13A+0.125*(m16A-m17A); nread = neighborList[n+7*Np]; distA[nread] = fq; // q = 8 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A-0.1*(jxA+jyA)-0.025*(m4A+m6A) +mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A+0.25*m13A+0.125*(m17A-m16A); - //f8 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A-0.1*(jxA+jyA)-0.025*(m4A+m6A) +mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A+0.25*m13A+0.125*(m17A-m16A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A-0.1*(jxA+jyA)-0.025*(m4A+m6A) +mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A+0.25*m13A+0.125*(m17A-m16A); nread = neighborList[n+6*Np]; distA[nread] = fq; // q = 9 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA-jyA)+0.025*(m4A-m6A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A-0.25*m13A+0.125*(m16A+m17A); - //f9 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA-jyA)+0.025*(m4A-m6A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A-0.25*m13A+0.125*(m16A+m17A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA-jyA)+0.025*(m4A-m6A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A-0.25*m13A+0.125*(m16A+m17A); nread = neighborList[n+9*Np]; distA[nread] = fq; // q = 10 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA-jxA)+0.025*(m6A-m4A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A-0.25*m13A-0.125*(m16A+m17A); - //f10 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA-jxA)+0.025*(m6A-m4A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A-0.25*m13A-0.125*(m16A+m17A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA-jxA)+0.025*(m6A-m4A)+mrt_V7*m9A+mrt_V11*m10A+mrt_V8*m11A+mrt_V12*m12A-0.25*m13A-0.125*(m16A+m17A); nread = neighborList[n+8*Np]; distA[nread] = fq; // q = 11 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA+jzA)+0.025*(m4A+m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A+0.25*m15A+0.125*(m18A-m16A); - //f11 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA+jzA)+0.025*(m4A+m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A+0.25*m15A+0.125*(m18A-m16A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA+jzA)+0.025*(m4A+m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A+0.25*m15A+0.125*(m18A-m16A); nread = neighborList[n+11*Np]; distA[nread] = fq; // q = 12 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A-0.1*(jxA+jzA)-0.025*(m4A+m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A+0.25*m15A+0.125*(m16A-m18A); - //f12 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A-0.1*(jxA+jzA)-0.025*(m4A+m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A+0.25*m15A+0.125*(m16A-m18A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A-0.1*(jxA+jzA)-0.025*(m4A+m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A+0.25*m15A+0.125*(m16A-m18A); nread = neighborList[n+10*Np]; distA[nread]= fq; // q = 13 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA-jzA)+0.025*(m4A-m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A-0.25*m15A-0.125*(m16A+m18A); - //f13 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA-jzA)+0.025*(m4A-m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A-0.25*m15A-0.125*(m16A+m18A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jxA-jzA)+0.025*(m4A-m8A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A-0.25*m15A-0.125*(m16A+m18A); nread = neighborList[n+13*Np]; distA[nread] = fq; // q= 14 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jzA-jxA)+0.025*(m8A-m4A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A-0.25*m15A+0.125*(m16A+m18A); - //f14 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jzA-jxA)+0.025*(m8A-m4A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A-0.25*m15A+0.125*(m16A+m18A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jzA-jxA)+0.025*(m8A-m4A)+mrt_V7*m9A+mrt_V11*m10A-mrt_V8*m11A-mrt_V12*m12A-0.25*m15A+0.125*(m16A+m18A); nread = neighborList[n+12*Np]; distA[nread] = fq; // q = 15 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA+jzA)+0.025*(m6A+m8A)-mrt_V6*m9A-mrt_V7*m10A+0.25*m14A+0.125*(m17A-m18A); - //f15 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA+jzA)+0.025*(m6A+m8A)-mrt_V6*m9A-mrt_V7*m10A+0.25*m14A+0.125*(m17A-m18A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA+jzA)+0.025*(m6A+m8A)-mrt_V6*m9A-mrt_V7*m10A+0.25*m14A+0.125*(m17A-m18A); nread = neighborList[n+15*Np]; distA[nread] = fq; // q = 16 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A-0.1*(jyA+jzA)-0.025*(m6A+m8A)-mrt_V6*m9A-mrt_V7*m10A+0.25*m14A+0.125*(m18A-m17A); - //f16 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A-0.1*(jyA+jzA)-0.025*(m6A+m8A)-mrt_V6*m9A-mrt_V7*m10A+0.25*m14A+0.125*(m18A-m17A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A-0.1*(jyA+jzA)-0.025*(m6A+m8A)-mrt_V6*m9A-mrt_V7*m10A+0.25*m14A+0.125*(m18A-m17A); nread = neighborList[n+14*Np]; distA[nread] = fq; // q = 17 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA-jzA)+0.025*(m6A-m8A)-mrt_V6*m9A-mrt_V7*m10A-0.25*m14A+0.125*(m17A+m18A); - //f17 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA-jzA)+0.025*(m6A-m8A)-mrt_V6*m9A-mrt_V7*m10A-0.25*m14A+0.125*(m17A+m18A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jyA-jzA)+0.025*(m6A-m8A)-mrt_V6*m9A-mrt_V7*m10A-0.25*m14A+0.125*(m17A+m18A); nread = neighborList[n+17*Np]; distA[nread] = fq; // q = 18 - fq = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jzA-jyA)+0.025*(m8A-m6A)-mrt_V6*m9A-mrt_V7*m10A-0.25*m14A-0.125*(m17A+m18A); - //f18 = mrt_V1*rhoA_next+mrt_V9*m1A+mrt_V10*m2A+0.1*(jzA-jyA)+0.025*(m8A-m6A)-mrt_V6*m9A-mrt_V7*m10A-0.25*m14A-0.125*(m17A+m18A); + fq = mrt_V1*rhoA+mrt_V9*m1A+mrt_V10*m2A+0.1*(jzA-jyA)+0.025*(m8A-m6A)-mrt_V6*m9A-mrt_V7*m10A-0.25*m14A-0.125*(m17A+m18A); nread = neighborList[n+16*Np]; distA[nread] = fq; //........................................................................ - - //Den[n] = f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; - - -// //..............carry out relaxation process............................................... -// m1 = m1 + rlx_setA*((-30*Den+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*Den - 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*Den) - m4) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fx); -// jy = jy + Fy; -// m6 = m6 + rlx_setB*((-0.6666666666666666*uy*Den) - m6) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fy); -// jz = jz + Fz; -// m8 = m8 + rlx_setB*((-0.6666666666666666*uz*Den) - m8) -// + (1-0.5*rlx_setB)*(-0.6666666666666666*Fz); -// m9 = m9 + rlx_setA*((Den*(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*Den*((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*((Den*(uy*uy-uz*uz)/porosity) - m11) -// + (1-0.5*rlx_setA)*(2*Fy*uy-2*Fz*uz)/porosity; -// m12 = m12 + rlx_setA*(-0.5*(Den*(uy*uy-uz*uz)/porosity)- m12) -// + (1-0.5*rlx_setA)*(-Fy*uy+Fz*uz)/porosity; -// m13 = m13 + rlx_setA*((Den*ux*uy/porosity) - m13) -// + (1-0.5*rlx_setA)*(Fy*ux+Fx*uy)/porosity; -// m14 = m14 + rlx_setA*((Den*uy*uz/porosity) - m14) -// + (1-0.5*rlx_setA)*(Fz*uy+Fy*uz)/porosity; -// m15 = m15 + rlx_setA*((Den*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); -// //....................................................................................................... - // ------------------- Fluid Component B -----------------------// rlx_setA = 1.0/tauB; rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA); //-------------------- MRT collison where body force has NO higher-order terms -------------// //..............carry out relaxation process............................................... //TODO need to incoporate porosity - m1B = m1B + rlx_setA*((19*rhoB*(ux*ux+uy*uy+uz*uz) - 11*rhoB_next) - m1B) + m1B = m1B + rlx_setA*((19*rhoB*(ux*ux+uy*uy+uz*uz) - 11*rhoB) - m1B) + (1-0.5*rlx_setA)*38*(FxB*ux+FyB*uy+FzB*uz); - m2B = m2B + rlx_setA*((3*rhoB_next - 5.5*rhoB*(ux*ux+uy*uy+uz*uz))- m2B) + m2B = m2B + rlx_setA*((3*rhoB - 5.5*rhoB*(ux*ux+uy*uy+uz*uz))- m2B) + (1-0.5*rlx_setA)*11*(-FxB*ux-FyB*uy-FzB*uz); jxB = jxB + FxB; m4B = m4B + rlx_setB*((-0.6666666666666666*ux*rhoB)- m4B) @@ -1012,140 +925,116 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double // ------------------- Fluid Component B -----------------------// //.................inverse transformation...................................................... // q=0 - fq = mrt_V1*rhoB_next-mrt_V2*m1B+mrt_V3*m2B; - //f0 = mrt_V1*rhoB_next-mrt_V2*m1B+mrt_V3*m2B; + fq = mrt_V1*rhoB-mrt_V2*m1B+mrt_V3*m2B; distB[n] = fq; // q = 1 - fq = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(jxB-m4B)+mrt_V6*(m9B-m10B); - //f1 = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(jxB-m4B)+mrt_V6*(m9B-m10B); + fq = mrt_V1*rhoB-mrt_V4*m1B-mrt_V5*m2B+0.1*(jxB-m4B)+mrt_V6*(m9B-m10B); nread = neighborList[n+Np]; distB[nread] = fq; // q=2 - fq = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(m4B-jxB)+mrt_V6*(m9B-m10B); - //f2 = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(m4B-jxB)+mrt_V6*(m9B-m10B); + fq = mrt_V1*rhoB-mrt_V4*m1B-mrt_V5*m2B+0.1*(m4B-jxB)+mrt_V6*(m9B-m10B); nread = neighborList[n]; distB[nread] = fq; // q = 3 - fq = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(jyB-m6B)+mrt_V7*(m10B-m9B)+mrt_V8*(m11B-m12B); - //f3 = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(jyB-m6B)+mrt_V7*(m10B-m9B)+mrt_V8*(m11B-m12B); + fq = mrt_V1*rhoB-mrt_V4*m1B-mrt_V5*m2B+0.1*(jyB-m6B)+mrt_V7*(m10B-m9B)+mrt_V8*(m11B-m12B); nread = neighborList[n+3*Np]; distB[nread] = fq; // q = 4 - fq = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(m6B-jyB)+mrt_V7*(m10B-m9B)+mrt_V8*(m11B-m12B); - //f4 = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(m6B-jyB)+mrt_V7*(m10B-m9B)+mrt_V8*(m11B-m12B); + fq = mrt_V1*rhoB-mrt_V4*m1B-mrt_V5*m2B+0.1*(m6B-jyB)+mrt_V7*(m10B-m9B)+mrt_V8*(m11B-m12B); nread = neighborList[n+2*Np]; distB[nread] = fq; // q = 5 - fq = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(jzB-m8B)+mrt_V7*(m10B-m9B)+mrt_V8*(m12B-m11B); - //f5 = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(jzB-m8B)+mrt_V7*(m10B-m9B)+mrt_V8*(m12B-m11B); + fq = mrt_V1*rhoB-mrt_V4*m1B-mrt_V5*m2B+0.1*(jzB-m8B)+mrt_V7*(m10B-m9B)+mrt_V8*(m12B-m11B); nread = neighborList[n+5*Np]; distB[nread] = fq; // q = 6 - fq = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(m8B-jzB)+mrt_V7*(m10B-m9B)+mrt_V8*(m12B-m11B); - //f6 = mrt_V1*rhoB_next-mrt_V4*m1B-mrt_V5*m2B+0.1*(m8B-jzB)+mrt_V7*(m10B-m9B)+mrt_V8*(m12B-m11B); + fq = mrt_V1*rhoB-mrt_V4*m1B-mrt_V5*m2B+0.1*(m8B-jzB)+mrt_V7*(m10B-m9B)+mrt_V8*(m12B-m11B); nread = neighborList[n+4*Np]; distB[nread] = fq; // q = 7 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB+jyB)+0.025*(m4B+m6B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B+0.25*m13B+0.125*(m16B-m17B); - //f7 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB+jyB)+0.025*(m4B+m6B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B+0.25*m13B+0.125*(m16B-m17B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB+jyB)+0.025*(m4B+m6B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B+0.25*m13B+0.125*(m16B-m17B); nread = neighborList[n+7*Np]; distB[nread] = fq; // q = 8 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B-0.1*(jxB+jyB)-0.025*(m4B+m6B) +mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B+0.25*m13B+0.125*(m17B-m16B); - //f8 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B-0.1*(jxB+jyB)-0.025*(m4B+m6B) +mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B+0.25*m13B+0.125*(m17B-m16B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B-0.1*(jxB+jyB)-0.025*(m4B+m6B) +mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B+0.25*m13B+0.125*(m17B-m16B); nread = neighborList[n+6*Np]; distB[nread] = fq; // q = 9 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB-jyB)+0.025*(m4B-m6B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B-0.25*m13B+0.125*(m16B+m17B); - //f9 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB-jyB)+0.025*(m4B-m6B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B-0.25*m13B+0.125*(m16B+m17B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB-jyB)+0.025*(m4B-m6B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B-0.25*m13B+0.125*(m16B+m17B); nread = neighborList[n+9*Np]; distB[nread] = fq; // q = 10 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB-jxB)+0.025*(m6B-m4B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B-0.25*m13B-0.125*(m16B+m17B); - //f10 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB-jxB)+0.025*(m6B-m4B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B-0.25*m13B-0.125*(m16B+m17B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB-jxB)+0.025*(m6B-m4B)+mrt_V7*m9B+mrt_V11*m10B+mrt_V8*m11B+mrt_V12*m12B-0.25*m13B-0.125*(m16B+m17B); nread = neighborList[n+8*Np]; distB[nread] = fq; // q = 11 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB+jzB)+0.025*(m4B+m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B+0.25*m15B+0.125*(m18B-m16B); - //f11 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB+jzB)+0.025*(m4B+m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B+0.25*m15B+0.125*(m18B-m16B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB+jzB)+0.025*(m4B+m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B+0.25*m15B+0.125*(m18B-m16B); nread = neighborList[n+11*Np]; distB[nread] = fq; // q = 12 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B-0.1*(jxB+jzB)-0.025*(m4B+m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B+0.25*m15B+0.125*(m16B-m18B); - //f12 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B-0.1*(jxB+jzB)-0.025*(m4B+m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B+0.25*m15B+0.125*(m16B-m18B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B-0.1*(jxB+jzB)-0.025*(m4B+m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B+0.25*m15B+0.125*(m16B-m18B); nread = neighborList[n+10*Np]; distB[nread]= fq; // q = 13 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB-jzB)+0.025*(m4B-m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B-0.25*m15B-0.125*(m16B+m18B); - //f13 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB-jzB)+0.025*(m4B-m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B-0.25*m15B-0.125*(m16B+m18B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jxB-jzB)+0.025*(m4B-m8B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B-0.25*m15B-0.125*(m16B+m18B); nread = neighborList[n+13*Np]; distB[nread] = fq; // q= 14 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jzB-jxB)+0.025*(m8B-m4B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B-0.25*m15B+0.125*(m16B+m18B); - //f14 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jzB-jxB)+0.025*(m8B-m4B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B-0.25*m15B+0.125*(m16B+m18B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jzB-jxB)+0.025*(m8B-m4B)+mrt_V7*m9B+mrt_V11*m10B-mrt_V8*m11B-mrt_V12*m12B-0.25*m15B+0.125*(m16B+m18B); nread = neighborList[n+12*Np]; distB[nread] = fq; // q = 15 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB+jzB)+0.025*(m6B+m8B)-mrt_V6*m9B-mrt_V7*m10B+0.25*m14B+0.125*(m17B-m18B); - //f15 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB+jzB)+0.025*(m6B+m8B)-mrt_V6*m9B-mrt_V7*m10B+0.25*m14B+0.125*(m17B-m18B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB+jzB)+0.025*(m6B+m8B)-mrt_V6*m9B-mrt_V7*m10B+0.25*m14B+0.125*(m17B-m18B); nread = neighborList[n+15*Np]; distB[nread] = fq; // q = 16 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B-0.1*(jyB+jzB)-0.025*(m6B+m8B)-mrt_V6*m9B-mrt_V7*m10B+0.25*m14B+0.125*(m18B-m17B); - //f16 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B-0.1*(jyB+jzB)-0.025*(m6B+m8B)-mrt_V6*m9B-mrt_V7*m10B+0.25*m14B+0.125*(m18B-m17B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B-0.1*(jyB+jzB)-0.025*(m6B+m8B)-mrt_V6*m9B-mrt_V7*m10B+0.25*m14B+0.125*(m18B-m17B); nread = neighborList[n+14*Np]; distB[nread] = fq; // q = 17 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB-jzB)+0.025*(m6B-m8B)-mrt_V6*m9B-mrt_V7*m10B-0.25*m14B+0.125*(m17B+m18B); - //f17 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB-jzB)+0.025*(m6B-m8B)-mrt_V6*m9B-mrt_V7*m10B-0.25*m14B+0.125*(m17B+m18B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jyB-jzB)+0.025*(m6B-m8B)-mrt_V6*m9B-mrt_V7*m10B-0.25*m14B+0.125*(m17B+m18B); nread = neighborList[n+17*Np]; distB[nread] = fq; // q = 18 - fq = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jzB-jyB)+0.025*(m8B-m6B)-mrt_V6*m9B-mrt_V7*m10B-0.25*m14B-0.125*(m17B+m18B); - //f18 = mrt_V1*rhoB_next+mrt_V9*m1B+mrt_V10*m2B+0.1*(jzB-jyB)+0.025*(m8B-m6B)-mrt_V6*m9B-mrt_V7*m10B-0.25*m14B-0.125*(m17B+m18B); + fq = mrt_V1*rhoB+mrt_V9*m1B+mrt_V10*m2B+0.1*(jzB-jyB)+0.025*(m8B-m6B)-mrt_V6*m9B-mrt_V7*m10B-0.25*m14B-0.125*(m17B+m18B); nread = neighborList[n+16*Np]; distB[nread] = fq; //........................................................................ - //Den[n+Np] = f0+f2+f1+f4+f3+f6+f5+f8+f7+f10+f9+f12+f11+f14+f13+f16+f15+f18+f17; - //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] = (rhoA+rhoB+Gsc*rhoA*rhoB)/3.0; - //Update density - //Den[n] = rhoA_next; - //Den[n+Np] = rhoB_next; - + Pressure[n] = (nA+nB+Gsc*nA*nB)/3.0; } } } -__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +__global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(int *Map,double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np){ + int ijk; int n; double vx,vy,vz,v_mag; double ux_A,uy_A,uz_A,ux_B,uy_B,uz_B,u_mag; @@ -1154,7 +1043,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(double *distA, double *d double jxA,jyA,jzA; double jxB,jyB,jzB; double rhoA,rhoB; - double rhoA_next,rhoB_next; + double nA,nB; // non-conserved moments double m1A,m2A,m4A,m6A,m8A,m9A,m10A,m11A,m12A,m13A,m14A,m15A,m16A,m17A,m18A; double m1B,m2B,m4B,m6B,m8B,m9B,m10B,m11B,m12B,m13B,m14B,m15B,m16B,m17B,m18B; @@ -1170,8 +1059,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(double *distA, double *d double FxA, FyA, FzA;//The total body force including Brinkman force and user-specified (Gx,Gy,Gz) double FxB, FyB, FzB; double rlx_setA,rlx_setB; - double rhoA_gradx,rhoA_grady,rhoA_gradz; - double rhoB_gradx,rhoB_grady,rhoB_gradz; + double nA_gradx,nA_grady,nA_gradz; + double nB_gradx,nB_grady,nB_gradz; double GffA_x,GffA_y,GffA_z; double GfsA_x,GfsA_y,GfsA_z; double GffB_x,GffB_y,GffB_z; @@ -1198,18 +1087,19 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(double *distA, double *d if ( n>>(distA,distB,Den,Np); +__global__ void dvc_ScaLBL_D3Q19_GreyscaleSC_Gradient(int *neighborList, int *Map, double *Den, double *DenGrad, int strideY, int strideZ,int start, int finish, int Np){ + + int n,nn; + int ijk; + // distributions + double m1,m2,m3,m4,m5,m6,m7,m8,m9; + double m10,m11,m12,m13,m14,m15,m16,m17,m18; + double nx,ny,nz; + + int S = Np/NBLOCKS/NTHREADS + 1; + for (int s=0; s>>(Map,distA,distB,DenA,DenB,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q19_GreyscaleSC_Init: %s \n",cudaGetErrorString(err)); } } -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(int *NeighborList, double *distA, double *distB, double *Den, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(int *NeighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_Density<<>>(NeighborList, distA, distB, Den, start, finish, Np); + dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_Density<<>>(NeighborList, Map, distA, distB, DenA, DenB, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -3358,9 +3233,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(int *NeighborList, double } } -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(double *distA, double *distB, double *Den, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(int *Map, double *distA, double *distB, double *DenA, double *DenB, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_Density<<>>(distA, distB, Den, start, finish, Np); + dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_Density<<>>(Map,distA, distB, DenA, DenB, start, finish, Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -3369,12 +3244,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(double *distA, double *d } -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, int *Map, double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT<<>>(neighborList,distA,distB,Den,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, + dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT<<>>(neighborList,Map,distA,distB,DenA,DenB,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, tauA,tauB,tauA_eff,tauB_eff,Gsc,Gx,Gy,Gz,start,finish,Np); cudaError_t err = cudaGetLastError(); @@ -3383,12 +3258,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(int *neighborList, double *di } } -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(int *Map,double *distA, double *distB, double *DenA,double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT<<>>(distA,distB,Den,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, + dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT<<>>(Map,distA,distB,DenA,DenB,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, tauA,tauB,tauA_eff,tauB_eff,Gsc,Gx,Gy,Gz,start,finish,Np); cudaError_t err = cudaGetLastError(); @@ -3397,12 +3272,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(double *distA, double *distB } } -extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(int *neighborList, double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(int *neighborList, int *Map, double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK<<>>(neighborList,distA,distB,Den,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, + dvc_ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK<<>>(neighborList,Map,distA,distB,DenA,DenB,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, tauA,tauB,tauA_eff,tauB_eff,Gsc,Gx,Gy,Gz,start,finish,Np); cudaError_t err = cudaGetLastError(); @@ -3411,12 +3286,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(int *neighborList, double *di } } -extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(double *distA, double *distB, double *Den, double *DenGradA, double *DenGradB, +extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(int *Map,double *distA, double *distB, double *DenA, double *DenB, double *DenGradA, double *DenGradB, double *SolidForceA, double *SolidForceB, double *Poros,double *Perm, double *Velocity,double *Pressure, double tauA,double tauB,double tauA_eff,double tauB_eff, double Gsc, double Gx, double Gy, double Gz, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK<<>>(distA,distB,Den,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, + dvc_ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK<<>>(Map,distA,distB,DenA,DenB,DenGradA,DenGradB,SolidForceA,SolidForceB,Poros,Perm,Velocity,Pressure, tauA,tauB,tauA_eff,tauB_eff,Gsc,Gx,Gy,Gz,start,finish,Np); cudaError_t err = cudaGetLastError(); @@ -3424,3 +3299,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(double *distA, double *distB printf("CUDA error in ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK: %s \n",cudaGetErrorString(err)); } } + +extern "C" void ScaLBL_D3Q19_GreyscaleSC_Gradient(int *neighborList, int *Map, double *Den, double *DenGrad, int strideY, int strideZ,int start, int finish, int Np){ + + dvc_ScaLBL_D3Q19_GreyscaleSC_Gradient<<>>(neighborList, Map, Den, DenGrad, strideY, strideZ,start, finish, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_D3Q19_GreyscaleSC_Gradient: %s \n",cudaGetErrorString(err)); + } +} diff --git a/models/GreyscaleSCModel.cpp b/models/GreyscaleSCModel.cpp index 7e437b7b..3f50dc34 100644 --- a/models/GreyscaleSCModel.cpp +++ b/models/GreyscaleSCModel.cpp @@ -465,6 +465,81 @@ void ScaLBL_GreyscaleSCModel::AssignGreyscaleAndSolidLabels() delete [] Dst; } +//void ScaLBL_GreyscaleSCModel::Density_Init(){ +// +// size_t NLABELS=0; +// signed char VALUE=0; +// +// vector LabelList{1,2}; +// vector SwList{0.0,1.0}; +// +// if (greyscaleSC_db->keyExists( "GreyNodeLabels" )){ +// LabelList.clear(); +// LabelList = greyscaleSC_db->getVector( "GreyNodeLabels" ); +// } +// if (greyscaleSC_db->keyExists( "GreyNodeSw" )){ +// SwList.clear(); +// SwList = greyscaleSC_db->getVector( "GreyNodeSw" ); +// } +// +// NLABELS=LabelList.size(); +// if (NLABELS != SwList.size()){ +// ERROR("Error: GreyNodeLabels and GreyNodeSw must be the same length! \n"); +// } +// +// 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; +// +// //double *Phi_temp; +// //Phi_temp=new double [Np]; +// //double phi = 0.0; +// +// for (int k=0; kid[n]; +// if (VALUE>0){ +// for (unsigned int idx=0; idx < NLABELS; idx++){ +// if (VALUE == LabelList[idx]){ +// double Sw = SwList[idx]; +// if ((Sw<0.0) || (Sw>1.0)) ERROR("Error: Initial saturation for grey nodes must be between [0.0, 1.0]! \n"); +// nB=Sw; +// nA=1.0-Sw; +// //phi = nA-nB; +// idx = NLABELS; +// } +// } +// if (VALUE==1){//label=1 reserved for NW phase +// //TODO; maybe need rho_major and rho_minor initialization +// nA=rhoA; +// nB=rhoB_minor; +// //phi = nA-nB; +// } +// else if(VALUE==2){//label=2 reserved for W phase +// //TODO; maybe need rho_major and rho_minor initialization +// nA=rhoA_minor; +// nB=rhoB; +// //phi = nA-nB; +// } +// int idx = Map(i,j,k); +// 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(Phi, Phi_temp, 1*Np*sizeof(double)); +// ScaLBL_DeviceBarrier(); +// delete [] Den_temp; +// //delete [] Phi_temp; +//} + void ScaLBL_GreyscaleSCModel::Density_Init(){ size_t NLABELS=0; @@ -487,10 +562,11 @@ void ScaLBL_GreyscaleSCModel::Density_Init(){ ERROR("Error: GreyNodeLabels and GreyNodeSw must be the same length! \n"); } - 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; + double *DenA_temp,*DenB_temp; + DenA_temp=new double [Nx*Ny*Nz]; + DenB_temp=new double [Nx*Ny*Nz]; + double nA=0.0;//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.0; //double *Phi_temp; //Phi_temp=new double [Np]; @@ -513,30 +589,32 @@ void ScaLBL_GreyscaleSCModel::Density_Init(){ } } if (VALUE==1){//label=1 reserved for NW phase - //TODO; maybe need rho_major and rho_minor initialization nA=rhoA; nB=rhoB_minor; //phi = nA-nB; } else if(VALUE==2){//label=2 reserved for W phase - //TODO; maybe need rho_major and rho_minor initialization nA=rhoA_minor; nB=rhoB; //phi = nA-nB; } - int idx = Map(i,j,k); - Den_temp[idx+0*Np] = nA; - Den_temp[idx+1*Np] = nB; - //Phi_temp[idx] = phi; + DenA_temp[n] = nA; + DenB_temp[n] = nB; + } + else{ //for ID<=0, i.e. all sorts of solid minerals, density is zero + DenA_temp[n] = 0.0; + DenB_temp[n] = 0.0; } } } } //copy to device - ScaLBL_CopyToDevice(Den, Den_temp, 2*Np*sizeof(double)); + ScaLBL_CopyToDevice(DenA, DenA_temp, Nx*Ny*Nz*sizeof(double)); + ScaLBL_CopyToDevice(DenB, DenB_temp, Nx*Ny*Nz*sizeof(double)); //ScaLBL_CopyToDevice(Phi, Phi_temp, 1*Np*sizeof(double)); ScaLBL_DeviceBarrier(); - delete [] Den_temp; + delete [] DenA_temp; + delete [] DenB_temp; //delete [] Phi_temp; } @@ -559,6 +637,7 @@ void ScaLBL_GreyscaleSCModel::Create(){ // Create a communicator for the device (will use optimized layout) // ScaLBL_Communicator ScaLBL_Comm(Mask); // original ScaLBL_Comm = std::shared_ptr(new ScaLBL_Communicator(Mask)); + ScaLBL_Comm_Regular = std::shared_ptr(new ScaLBL_Communicator(Mask)); int Npad=(Np/16 + 2)*16; if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); @@ -577,9 +656,11 @@ void ScaLBL_GreyscaleSCModel::Create(){ neighborSize=18*(Np*sizeof(int)); //........................................................................... ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np); ScaLBL_AllocateDeviceMemory((void **) &fqA, 19*dist_mem_size); ScaLBL_AllocateDeviceMemory((void **) &fqB, 19*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &Den, 2*sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &DenA, sizeof(double)*Nx*Ny*Nz); + ScaLBL_AllocateDeviceMemory((void **) &DenB, sizeof(double)*Nx*Ny*Nz); ScaLBL_AllocateDeviceMemory((void **) &Permeability, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Porosity, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Pressure_dvc, sizeof(double)*Np); @@ -592,6 +673,36 @@ void ScaLBL_GreyscaleSCModel::Create(){ // Update GPU data structures if (rank==0) printf ("Setting up device neighbor list \n"); fflush(stdout); + // Copy the Map to device + int *TmpMap; + TmpMap=new int[Np]; + for (int k=1; kLastExterior(); idx++){ + auto n = TmpMap[idx]; + if (n > Nx*Ny*Nz){ + printf("Bad value! idx=%i \n", n); + TmpMap[idx] = Nx*Ny*Nz-1; + } + } + for (int idx=ScaLBL_Comm->FirstInterior(); idxLastInterior(); idx++){ + auto n = TmpMap[idx]; + if ( n > Nx*Ny*Nz ){ + printf("Bad value! idx=%i \n",n); + TmpMap[idx] = Nx*Ny*Nz-1; + } + } + ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np); + ScaLBL_DeviceBarrier(); + delete [] TmpMap; // copy the neighbor list ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); } @@ -631,18 +742,20 @@ void ScaLBL_GreyscaleSCModel::Initialize(){ if (rank==0) printf ("Initializing density field \n"); Density_Init();//initialize density field if (rank==0) printf ("Initializing distributions \n"); - ScaLBL_D3Q19_GreyscaleSC_Init(fqA, fqB, Den, Np); + ScaLBL_D3Q19_GreyscaleSC_Init(dvcMap,fqA, fqB, DenA,DenB, Np); //debug DoubleArray PhaseField(Nx,Ny,Nz); - ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); + //ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); + ScaLBL_CopyToHost(PhaseField.data(), DenA, sizeof(double)*N); FILE *AFILE; sprintf(LocalRankFilename,"A_init.%05i.raw",rank); AFILE = fopen(LocalRankFilename,"wb"); fwrite(PhaseField.data(),8,N,AFILE); fclose(AFILE); - ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField); + //ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField); + ScaLBL_CopyToHost(PhaseField.data(), DenB, sizeof(double)*N); FILE *BFILE; sprintf(LocalRankFilename,"B_init.%05i.raw",rank); BFILE = fopen(LocalRankFilename,"wb"); @@ -707,49 +820,32 @@ void ScaLBL_GreyscaleSCModel::Run(){ // Compute the density field // Read for Aq, Bq happens in this routine (requires communication) ScaLBL_Comm->BiSendD3Q19AA(fqA,fqB); //READ FROM NORMAL - ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(NeighborList, fqA, fqB, Den, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(NeighborList, dvcMap, fqA, fqB, DenA, DenB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->BiRecvD3Q19AA(fqA,fqB); //WRITE INTO OPPOSITE ScaLBL_DeviceBarrier(); - ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(NeighborList, fqA, fqB, Den, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q19_AAodd_GreyscaleSC_Density(NeighborList, dvcMap, fqA, fqB, DenA, DenB, 0, ScaLBL_Comm->LastExterior(), Np); + // Compute density gradient // fluid component A - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[0], DenGradA, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->SendHalo(DenA); + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenA, DenGradA, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(DenA); + ScaLBL_DeviceBarrier(); + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenA, DenGradA, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); // fluid component B - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[Np], DenGradB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - // Compute density gradient - ScaLBL_Comm->SendHalo(&Den[0]); - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[0], DenGradA, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->RecvGrad(&Den[0],DenGradA); + ScaLBL_Comm_Regular->SendHalo(DenB); + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenB, DenGradB, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(DenB); ScaLBL_DeviceBarrier(); - ScaLBL_Comm->SendHalo(&Den[Np]); - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[Np], DenGradB, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->RecvGrad(&Den[Np],DenGradB); - ScaLBL_DeviceBarrier(); - - //debug -// DoubleArray PhaseField(Nx,Ny,Nz); -// ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); -// FILE *AFILE; -// sprintf(LocalRankFilename,"A_beforeCol_time_%i.%05i.raw",timestep,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_beforeCol_time_%i.%05i.raw",timestep,rank); -// BFILE = fopen(LocalRankFilename,"wb"); -// fwrite(PhaseField.data(),8,N,BFILE); -// fclose(BFILE); - + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenB, DenGradB, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); // Collsion - ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(NeighborList, fqA, fqB, Den, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, + ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(NeighborList, dvcMap, fqA, fqB, DenA, DenB, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, tauA, tauB, tauA_eff, tauB_eff, Gsc, Fx, Fy, Fz, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); // Collsion - ScaLBL_D3Q19_AAodd_GreyscaleSC_BGK(NeighborList, fqA, fqB, Den, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, + ScaLBL_D3Q19_AAodd_GreyscaleSC_MRT(NeighborList, dvcMap, fqA, fqB, DenA, DenB, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, tauA, tauB, tauA_eff, tauB_eff, Gsc, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); @@ -760,48 +856,32 @@ void ScaLBL_GreyscaleSCModel::Run(){ // Compute the density field // Read for Aq, Bq happens in this routine (requires communication) ScaLBL_Comm->BiSendD3Q19AA(fqA,fqB); //READ FROM NORMAL - ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(fqA, fqB, Den, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(dvcMap, fqA, fqB, DenA, DenB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->BiRecvD3Q19AA(fqA,fqB); //WRITE INTO OPPOSITE ScaLBL_DeviceBarrier(); - ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(fqA, fqB, Den, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_D3Q19_AAeven_GreyscaleSC_Density(dvcMap, fqA, fqB, DenA, DenB, 0, ScaLBL_Comm->LastExterior(), Np); + // Compute density gradient // fluid component A - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[0], DenGradA, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->SendHalo(DenA); + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenA, DenGradA, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(DenA); + ScaLBL_DeviceBarrier(); + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenA, DenGradA, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); // fluid component B - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[Np], DenGradB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - // Compute density gradient - ScaLBL_Comm->SendHalo(&Den[0]); - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[0], DenGradA, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->RecvGrad(&Den[0],DenGradA); + ScaLBL_Comm_Regular->SendHalo(DenB); + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenB, DenGradB, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(DenB); ScaLBL_DeviceBarrier(); - ScaLBL_Comm->SendHalo(&Den[Np]); - ScaLBL_D3Q19_GreyscaleFE_Gradient(NeighborList, &Den[Np], DenGradB, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->RecvGrad(&Den[Np],DenGradB); - ScaLBL_DeviceBarrier(); - -// //debug -// //DoubleArray PhaseField(Nx,Ny,Nz); -// ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); -// //FILE *AFILE; -// sprintf(LocalRankFilename,"A_beforeCol_time_%i.%05i.raw",timestep,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_beforeCol_time_%i.%05i.raw",timestep,rank); -// BFILE = fopen(LocalRankFilename,"wb"); -// fwrite(PhaseField.data(),8,N,BFILE); -// fclose(BFILE); + ScaLBL_D3Q19_GreyscaleSC_Gradient(NeighborList, dvcMap, DenB, DenGradB, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); // Collsion - ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(fqA, fqB, Den, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, + ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(dvcMap,fqA, fqB, DenA, DenB, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, tauA, tauB, tauA_eff, tauB_eff, Gsc, Fx, Fy, Fz, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); // Collsion - ScaLBL_D3Q19_AAeven_GreyscaleSC_BGK(fqA, fqB, Den, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, + ScaLBL_D3Q19_AAeven_GreyscaleSC_MRT(dvcMap,fqA, fqB, DenA, DenB, DenGradA, DenGradB, SolidForceA, SolidForceB, Porosity,Permeability,Velocity,Pressure_dvc, tauA, tauB, tauA_eff, tauB_eff, Gsc, Fx, Fy, Fz, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); @@ -1114,14 +1194,16 @@ void ScaLBL_GreyscaleSCModel::WriteDebug(){ // fwrite(PhaseField.data(),8,N,OUTFILE); // fclose(OUTFILE); // - ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); + //ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); + ScaLBL_CopyToHost(PhaseField.data(), DenA, sizeof(double)*N); 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); + //ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField); + ScaLBL_CopyToHost(PhaseField.data(), DenB, sizeof(double)*N); FILE *BFILE; sprintf(LocalRankFilename,"B.%05i.raw",rank); BFILE = fopen(LocalRankFilename,"wb"); diff --git a/models/GreyscaleSCModel.h b/models/GreyscaleSCModel.h index 348490c4..bc71dc01 100644 --- a/models/GreyscaleSCModel.h +++ b/models/GreyscaleSCModel.h @@ -54,6 +54,7 @@ public: std::shared_ptr Dm; // this domain is for analysis std::shared_ptr Mask; // this domain is for lbm std::shared_ptr ScaLBL_Comm; + std::shared_ptr ScaLBL_Comm_Regular; // input database std::shared_ptr db; @@ -64,12 +65,13 @@ public: signed char *id; int *NeighborList; + int *dvcMap; double *fqA, *fqB; double *Permeability;//grey voxel permeability double *Porosity; double *Velocity; double *Pressure_dvc; - double *Den; + double *DenA, *DenB; double *DenGradA,*DenGradB; double *SolidForceA,*SolidForceB;