diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 73b89f1d..42c51525 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -180,7 +180,9 @@ extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *NeighborList, double *Phi, double // FREE ENERGY LEE MODEL -extern "C" void ScaLBL_D3Q19_FreeLeeModel_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np); +extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np); + +extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np); extern "C" void ScaLBL_FreeLeeModel_PhaseField_Init(int *Map, double *Phi, double *Den, double *hq, double *ColorGrad, double rhonA, double rhoB, double tauM, double W, int start, int finish, int Np); @@ -199,6 +201,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel(int *Map, double *dist, double double rhoA, double rhoB, double tauA, double tauB, double tauM, double kappa, double beta, double W, double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np); +extern "C" void ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(int *neighborList, double *dist, double *Vel, double *Pressure, + double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np); + +extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, + double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np); + // BOUNDARY CONDITION ROUTINES diff --git a/cpu/FreeLee.cpp b/cpu/FreeLee.cpp index f28af185..32a7b568 100644 --- a/cpu/FreeLee.cpp +++ b/cpu/FreeLee.cpp @@ -2,7 +2,7 @@ #define STOKES -extern "C" void ScaLBL_D3Q19_FreeLeeModel_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np) +extern "C" void ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(double *gqbar, double *mu_phi, double *ColorGrad, double Fx, double Fy, double Fz, int Np) { int n; double p = 1.0;//NOTE: take initial pressure p=1.0 @@ -15,7 +15,7 @@ extern "C" void ScaLBL_D3Q19_FreeLeeModel_Init(double *gqbar, double *mu_phi, do cg_y = ColorGrad[1*Np+n]; cg_z = ColorGrad[2*Np+n]; - gqbar[0*Np+n] = 0.3333333333333333; + gqbar[0*Np+n] = 0.3333333333333333*p; gqbar[1*Np+n] = 0.055555555555555555*(p - 0.5*(chem*cg_x+Fx)); //double(100*n)+1.f; gqbar[2*Np+n] = 0.055555555555555555*(p - 0.5*(-chem*cg_x-Fx)); //double(100*n)+2.f; gqbar[3*Np+n] = 0.055555555555555555*(p - 0.5*(chem*cg_y+Fy)); //double(100*n)+3.f; @@ -40,6 +40,38 @@ extern "C" void ScaLBL_D3Q19_FreeLeeModel_Init(double *gqbar, double *mu_phi, do } } +extern "C" void ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(double *gqbar, double Fx, double Fy, double Fz, int Np) +{ + int n; + double p = 1.0;//NOTE: take initial pressure p=1.0 + + for (n=0; n 10Np => odd part of dist) + m1 = dist[nr1]; // reading the f1 data into register fq + + nr2 = neighborList[n+Np]; // neighbor 1 ( < 10Np => even part of dist) + m2 = dist[nr2]; // reading the f2 data into register fq + + // q=3 + nr3 = neighborList[n+2*Np]; // neighbor 4 + m3 = dist[nr3]; + + // q = 4 + nr4 = neighborList[n+3*Np]; // neighbor 3 + m4 = dist[nr4]; + + // q=5 + nr5 = neighborList[n+4*Np]; + m5 = dist[nr5]; + + // q = 6 + nr6 = neighborList[n+5*Np]; + m6 = dist[nr6]; + + // q=7 + nr7 = neighborList[n+6*Np]; + m7 = dist[nr7]; + + // q = 8 + nr8 = neighborList[n+7*Np]; + m8 = dist[nr8]; + + // q=9 + nr9 = neighborList[n+8*Np]; + m9 = dist[nr9]; + + // q = 10 + nr10 = neighborList[n+9*Np]; + m10 = dist[nr10]; + + // q=11 + nr11 = neighborList[n+10*Np]; + m11 = dist[nr11]; + + // q=12 + nr12 = neighborList[n+11*Np]; + m12 = dist[nr12]; + + // q=13 + nr13 = neighborList[n+12*Np]; + m13 = dist[nr13]; + + // q=14 + nr14 = neighborList[n+13*Np]; + m14 = dist[nr14]; + + // q=15 + nr15 = neighborList[n+14*Np]; + m15 = dist[nr15]; + + // q=16 + nr16 = neighborList[n+15*Np]; + m16 = dist[nr16]; + + // q=17 + nr17 = neighborList[n+16*Np]; + m17 = dist[nr17]; + + // q=18 + nr18 = neighborList[n+17*Np]; + m18 = dist[nr18]; + + //compute fluid velocity + ux = 3.0/rho0*(m1-m2+m7-m8+m9-m10+m11-m12+m13-m14+0.5*(Fx)/3.0); + uy = 3.0/rho0*(m3-m4+m7-m8-m9+m10+m15-m16+m17-m18+0.5*(Fy)/3.0); + uz = 3.0/rho0*(m5-m6+m11-m12-m13+m14+m15-m16-m17+m18+0.5*(Fz)/3.0); + //compute pressure + p = (m0+m2+m1+m4+m3+m6+m5+m8+m7+m10+m9+m12+m11+m14+m13+m16+m15+m18+m17); + + //------------------------------------------------- BCK collison ------------------------------------------------------------// + // q=0 + dist[n] = m0 + 0.5*(Fx*ux + Fy*uy + Fz*uz)*(-0.6666666666666666 + ux*ux + uy*uy + uz*uz) - + (m0 - 0.3333333333333333*p + 0.25*(Fx*ux + Fy*uy + Fz*uz)* + (-0.6666666666666666 + ux*ux + uy*uy + uz*uz) + 0.16666666666666666*rho0*(ux*ux + uy*uy + uz*uz))/ + tau; + + // q = 1 + dist[nr2] = m1 + 0.25*(Fx*(-1 + ux) + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) - + (m1 - 0.05555555555555555*p + 0.08333333333333333*rho0* + (-(ux*ux) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*uz)) + + 0.125*(Fx*(-1. + ux) + Fy*uy + Fz*uz)* + (-0.2222222222222222 - 1.*(ux*ux) + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*uz)))/tau; + + // q=2 + dist[nr1] = m2 + 0.25*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - ux*ux + + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) - + (m2 - 0.05555555555555555*p + 0.08333333333333333*rho0* + (-(ux*ux) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*uz)) + + 0.125*(Fx + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(ux*ux) + + 0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*uz)))/tau; + + // q = 3 + dist[nr4] = m3 + 0.25*(Fx*ux + Fy*(-1 + uy) + Fz*uz)*(-0.2222222222222222 - uy*uy + + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) - + (m3 - 0.05555555555555555*p + 0.08333333333333333*rho0* + (-(uy*uy) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*uz)) + + 0.125*(Fx*ux + Fy*(-1. + uy) + Fz*uz)* + (-0.2222222222222222 - 1.*(uy*uy) + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*uz)))/tau; + + // q = 4 + dist[nr3] = m4 + 0.25*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uy*uy + + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) - + (m4 - 0.05555555555555555*p + 0.08333333333333333*rho0* + (-(uy*uy) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*uz)) + + 0.125*(Fy + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uy*uy) + + 0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*uz)))/tau; + + // q = 5 + dist[nr6] = m5 + 0.25*(Fx*ux + Fy*uy + Fz*(-1 + uz))*(-0.2222222222222222 - uz*uz + + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) - + (m5 - 0.05555555555555555*p + 0.08333333333333333*rho0* + (-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2 + uz)*uz)) + + 0.125*(Fx*ux + Fy*uy + Fz*(-1. + uz))* + (-0.2222222222222222 - 1.*(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + (-2. + uz)*uz)))/tau; + + // q = 6 + dist[nr5] = m6 + 0.25*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - uz*uz + + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) - + (m6 - 0.05555555555555555*p + 0.08333333333333333*rho0* + (-(uz*uz) + 0.3333333333333333*(ux*ux + uy*uy + uz*(2 + uz))) + + 0.125*(Fz + Fx*ux + Fy*uy + Fz*uz)*(-0.2222222222222222 - 1.*(uz*uz) + + 0.3333333333333333*(ux*ux + uy*uy + uz*(2. + uz))))/tau; + + // q = 7 + dist[nr8] = m7 - 0.125*(Fx*(-1 + ux) + Fy*(-1 + uy) + Fz*uz)* + (0.2222222222222222 + (ux + uy)*(ux + uy) - 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\ + - (m7 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux + uy)*(ux + uy)) + 0.3333333333333333*(-2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + + 0.0625*(Fx*(-1. + ux) + Fy*(-1. + uy) + Fz*uz)* + (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + + 0.3333333333333333*(-2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau; + + // q = 8 + dist[nr7] = m8 + 0.125*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)* + (-0.2222222222222222 - (ux + uy)*(ux + uy) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz))\ + - (m8 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux + uy)*(ux + uy)) + 0.3333333333333333*(2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + + 0.0625*(Fx + Fy + Fx*ux + Fy*uy + Fz*uz)* + (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uy - 1.*(uy*uy) + + 0.3333333333333333*(2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau; + + // q = 9 + dist[nr10] = m9 + 0.125*(Fy + Fx*(-1 + ux) + Fy*uy + Fz*uz)* + (-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + - (m9 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux - uy)*(ux - uy)) + 0.3333333333333333*(-2*ux + ux*ux + 2*uy + uy*uy + uz*uz)) + + 0.0625*(Fy + Fx*(-1. + ux) + Fy*uy + Fz*uz)* + (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + + 0.3333333333333333*(-2.*ux + ux*ux + 2.*uy + uy*uy + uz*uz)))/tau; + + // q = 10 + dist[nr9] = m10 + 0.125*(Fx*(1 + ux) + Fy*(-1 + uy) + Fz*uz)* + (-0.2222222222222222 - (ux - uy)*(ux - uy) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz))\ + - (m10 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux - uy)*(ux - uy)) + 0.3333333333333333*(2*ux + ux*ux - 2*uy + uy*uy + uz*uz)) + + 0.0625*(Fx*(1 + ux) + Fy*(-1. + uy) + Fz*uz)* + (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uy - 1.*(uy*uy) + + 0.3333333333333333*(2.*ux + ux*ux - 2.*uy + uy*uy + uz*uz)))/tau; + + // q = 11 + dist[nr12] = m11 - 0.125*(Fx*(-1 + ux) + Fy*uy + Fz*(-1 + uz))* + (0.2222222222222222 + (ux + uz)*(ux + uz) - 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\ + - (m11 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux + uz)*(ux + uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + + 0.0625*(Fx*(-1. + ux) + Fy*uy + Fz*(-1. + uz))* + (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau; + + // q = 12 + dist[nr11] = m12 + 0.125*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)* + (-0.2222222222222222 - (ux + uz)*(ux + uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz)))\ + - (m12 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux + uz)*(ux + uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + uz*(2 + uz))) + + 0.0625*(Fx + Fz + Fx*ux + Fy*uy + Fz*uz)* + (-0.2222222222222222 - 1.*(ux*ux) - 2.*ux*uz - 1.*(uz*uz) + + 0.3333333333333333*(2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau; + + // q = 13 + dist[nr14] = m13 + 0.125*(Fz + Fx*(-1 + ux) + Fy*uy + Fz*uz)* + (-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz)))\ + - (m13 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux - uz)*(ux - uz)) + 0.3333333333333333*(-2*ux + ux*ux + uy*uy + uz*(2 + uz))) + + 0.0625*(Fz + Fx*(-1. + ux) + Fy*uy + Fz*uz)* + (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + + 0.3333333333333333*(-2.*ux + ux*ux + uy*uy + uz*(2. + uz))))/tau; + + // q= 14 + dist[nr13] = m14 + 0.125*(Fx*(1 + ux) + Fy*uy + Fz*(-1 + uz))* + (-0.2222222222222222 - (ux - uz)*(ux - uz) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz))\ + - (m14 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((ux - uz)*(ux - uz)) + 0.3333333333333333*(2*ux + ux*ux + uy*uy + (-2 + uz)*uz)) + + 0.0625*(Fx*(1 + ux) + Fy*uy + Fz*(-1. + uz))* + (-0.2222222222222222 - 1.*(ux*ux) + 2.*ux*uz - 1.*(uz*uz) + + 0.3333333333333333*(2.*ux + ux*ux + uy*uy + (-2. + uz)*uz)))/tau; + + // q = 15 + dist[nr16] = m15 - 0.125*(Fx*ux + Fy*(-1 + uy) + Fz*(-1 + uz))* + (0.2222222222222222 + (uy + uz)*(uy + uz) - 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz))\ + - (m15 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + (-2 + uz)*uz)) + + 0.0625*(Fx*ux + Fy*(-1. + uy) + Fz*(-1. + uz))* + (-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + (-2. + uz)*uz)))/tau; + + // q = 16 + dist[nr15] = m16 + 0.125*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)* + (-0.2222222222222222 - (uy + uz)*(uy + uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz)))\ + - (m16 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((uy + uz)*(uy + uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + uz*(2 + uz))) + + 0.0625*(Fy + Fz + Fx*ux + Fy*uy + Fz*uz)* + (-0.2222222222222222 - 1.*(uy*uy) - 2.*uy*uz - 1.*(uz*uz) + + 0.3333333333333333*(ux*ux + 2.*uy + uy*uy + uz*(2. + uz))))/tau; + + // q = 17 + dist[nr18] = m17 + 0.125*(Fz + Fx*ux + Fy*(-1 + uy) + Fz*uz)* + (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz)))\ + - (m17 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux - 2*uy + uy*uy + uz*(2 + uz))) + + 0.0625*(Fz + Fx*ux + Fy*(-1. + uy) + Fz*uz)* + (-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + + 0.3333333333333333*(ux*ux - 2.*uy + uy*uy + uz*(2. + uz))))/tau; + + // q = 18 + dist[nr17] = m18 + 0.125*(Fx*ux + Fy*(1 + uy) + Fz*(-1 + uz))* + (-0.2222222222222222 - (uy - uz)*(uy - uz) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz))\ + - (m18 - 0.027777777777777776*p + 0.041666666666666664*rho0* + (-((uy - uz)*(uy - uz)) + 0.3333333333333333*(ux*ux + 2*uy + uy*uy + (-2 + uz)*uz)) + + 0.0625*(Fx*ux + Fy*(1 + uy) + Fz*(-1. + uz))* + (-0.2222222222222222 - 1.*(uy*uy) + 2.*uy*uz - 1.*(uz*uz) + + 0.3333333333333333*(ux*ux + 2.*uy + uy*uy + (-2. + uz)*uz)))/tau; + //----------------------------------------------------------------------------------------------------------------------------------------// + + + //Update velocity on device + Vel[0*Np+n] = ux; + Vel[1*Np+n] = uy; + Vel[2*Np+n] = uz; + //Update pressure on device + Pressure[n] = p; + } +} + +extern "C" void ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(double *dist, double *Vel, double *Pressure, + double tau, double rho0, double Fx, double Fy, double Fz, int start, int finish, int Np){ + + int n; + double ux,uy,uz;//fluid velocity + double p;//pressure + // distribution functions + double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; + double m0,m3,m5,m7; + + for (int n=start; nkeyExists( "timestepMax" )){ timestepMax = freelee_db->getScalar( "timestepMax" ); } + if (freelee_db->keyExists( "tau" )){//only for single-fluid Lee model + tau = freelee_db->getScalar( "tau" ); + } if (freelee_db->keyExists( "tauA" )){ tauA = freelee_db->getScalar( "tauA" ); } @@ -54,6 +60,9 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){ if (freelee_db->keyExists( "tauM" )){ tauM = freelee_db->getScalar( "tauM" ); } + if (freelee_db->keyExists( "rho0" )){ + rho0 = freelee_db->getScalar( "rho0" ); + } if (freelee_db->keyExists( "rhoA" )){ rhoA = freelee_db->getScalar( "rhoA" ); } @@ -193,9 +202,9 @@ void ScaLBL_FreeLeeModel::ReadInput(){ } -void ScaLBL_FreeLeeModel::Create(){ +void ScaLBL_FreeLeeModel::Create_TwoFluid(){ /* - * This function creates the variables needed to run a LBM + * This function creates the variables needed to run two-fluid Lee model */ //......................................................... // Initialize communication structures in averaging domain @@ -207,7 +216,7 @@ void ScaLBL_FreeLeeModel::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)); + //ScaLBL_Comm_Regular = std::shared_ptr(new ScaLBL_Communicator(Mask)); ScaLBL_Comm_WideHalo = std::shared_ptr(new ScaLBLWideHalo_Communicator(Mask,2)); // create the layout for the LBM @@ -276,6 +285,51 @@ void ScaLBL_FreeLeeModel::Create(){ delete [] neighborList; } +void ScaLBL_FreeLeeModel::Create_SingleFluid(){ + /* + * This function creates the variables needed to run single-fluid Lee model + */ + //......................................................... + // Initialize communication structures in averaging domain + for (int i=0; iid[i] = Mask->id[i]; + Mask->CommInit(); + Np=Mask->PoreCount(); + //........................................................................... + if (rank==0) printf ("Create ScaLBL_Communicator \n"); + // 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)); + + // create the layout for the LBM + int Npad=(Np/16 + 2)*16; + if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); + Map.resize(Nx,Ny,Nz); Map.fill(-2); + auto neighborList= new int[18*Npad]; + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1); + comm.barrier(); + + //........................................................................... + // MAIN VARIABLES ALLOCATED HERE + //........................................................................... + // LBM variables + if (rank==0) printf ("Allocating distributions \n"); + //......................device distributions................................. + dist_mem_size = Np*sizeof(double); + neighborSize=18*(Np*sizeof(int)); + //........................................................................... + ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **) &gqbar, 19*dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); + //........................................................................... + // Update GPU data structures + if (rank==0) printf ("Setting up device map and neighbor list \n"); + // copy the neighbor list + ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); + comm.barrier(); + delete [] neighborList; +} + void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad() { double *phase; @@ -466,21 +520,31 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad() ScaLBL_CopyToDevice(mu_phi, mu_phi_host, Np*sizeof(double)); ScaLBL_Comm->Barrier(); comm.barrier(); + + //debug + //save the phase field and check it + //FILE *OUTFILE; + //sprintf(LocalRankFilename,"Phase_Init.%05i.raw",rank); + //OUTFILE = fopen(LocalRankFilename,"wb"); + //fwrite(phase,8,Nh,OUTFILE); + //fclose(OUTFILE); + + delete [] phase; delete [] ColorGrad_host; delete [] mu_phi_host; delete [] Dst; } -void ScaLBL_FreeLeeModel::Initialize(){ +void ScaLBL_FreeLeeModel::Initialize_TwoFluid(){ /* - * This function initializes model + * This function initializes two-fluid Lee model */ if (rank==0) printf ("Initializing phase field, chemical potential and color gradient\n"); AssignComponentLabels_ChemPotential_ColorGrad();//initialize phase field Phi if (rank==0) printf ("Initializing distributions for momentum transport\n"); - ScaLBL_D3Q19_FreeLeeModel_Init(gqbar, mu_phi, ColorGrad, Fx, Fy, Fz, Np); + ScaLBL_D3Q19_FreeLeeModel_TwoFluid_Init(gqbar, mu_phi, ColorGrad, Fx, Fy, Fz, Np); if (rank==0) printf ("Initializing density field and distributions for phase-field transport\n"); ScaLBL_FreeLeeModel_PhaseField_Init(dvcMap, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np); @@ -568,7 +632,84 @@ void ScaLBL_FreeLeeModel::Initialize(){ //ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double)); } -void ScaLBL_FreeLeeModel::Run(){ +void ScaLBL_FreeLeeModel::Initialize_SingleFluid(){ + /* + * This function initializes single-fluid Lee model + */ + if (rank==0) printf ("Initializing distributions for momentum transport\n"); + ScaLBL_D3Q19_FreeLeeModel_SingleFluid_Init(gqbar, Fx, Fy, Fz, Np); + + if (Restart == true){ + //TODO need to revise this function + //remove the phase-related part + + + +// if (rank==0){ +// printf("Reading restart file! \n"); +// } +// +// // Read in the restart file to CPU buffers +// int *TmpMap; +// TmpMap = new int[Np]; +// +// double *cPhi, *cDist, *cDen; +// cPhi = new double[N]; +// cDen = new double[2*Np]; +// cDist = new double[19*Np]; +// ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int)); +// //ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double)); +// +// ifstream File(LocalRestartFile,ios::binary); +// int idx; +// double value,va,vb; +// for (int n=0; nLastExterior(); n++){ +// va = cDen[n]; +// vb = cDen[Np + n]; +// value = (va-vb)/(va+vb); +// idx = TmpMap[n]; +// if (!(idx < 0) && idxFirstInterior(); nLastInterior(); n++){ +// va = cDen[n]; +// vb = cDen[Np + n]; +// value = (va-vb)/(va+vb); +// idx = TmpMap[n]; +// if (!(idx < 0) && idxBarrier(); +// comm.barrier(); +// +// if (rank==0) printf ("Initializing phase and density fields on device from Restart\n"); +// //TODO the following function is to be updated. +// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, 0, ScaLBL_Comm->LastExterior(), Np); +// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + } +} + +void ScaLBL_FreeLeeModel::Run_TwoFluid(){ int nprocs=nprocx*nprocy*nprocz; const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); @@ -694,8 +835,105 @@ void ScaLBL_FreeLeeModel::Run(){ // ************************************************************************ } +void ScaLBL_FreeLeeModel::Run_SingleFluid(){ + int nprocs=nprocx*nprocy*nprocz; + const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + + if (rank==0){ + printf("********************************************************\n"); + printf("No. of timesteps: %i \n", timestepMax); + fflush(stdout); + } -void ScaLBL_FreeLeeModel::WriteDebug(){ + //.......create and start timer............ + double starttime,stoptime,cputime; + ScaLBL_Comm->Barrier(); + comm.barrier(); + starttime = MPI_Wtime(); + //......................................... + + //************ MAIN ITERATION LOOP ***************************************/ + PROFILE_START("Loop"); + while (timestep < timestepMax ) { + //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); } + PROFILE_START("Update"); + // *************ODD TIMESTEP************* + timestep++; + //------------------------------------------------------------------------------------------------------------------- + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL + ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(NeighborList, gqbar, Velocity, Pressure, tau, rho0, Fx, Fy, Fz, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set boundary conditions + // TODO to be revised! + if (BoundaryCondition == 3){ + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, gqbar, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, gqbar, dout, timestep); + } + if (BoundaryCondition == 4){ + din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, gqbar, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, gqbar, dout, timestep); + } + else if (BoundaryCondition == 5){ + ScaLBL_Comm->D3Q19_Reflection_BC_z(gqbar); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar); + } + ScaLBL_D3Q19_AAodd_FreeLeeModel_SingleFluid_BGK(NeighborList, gqbar, Velocity, Pressure, tau, rho0, Fx, Fy, Fz, + 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + + // *************EVEN TIMESTEP************* + timestep++; + //------------------------------------------------------------------------------------------------------------------- + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL + ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(gqbar, Velocity, Pressure, tau, rho0, Fx, Fy, Fz, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set boundary conditions + // TODO to be revised! + if (BoundaryCondition == 3){ + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, gqbar, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, gqbar, dout, timestep); + } + else if (BoundaryCondition == 4){ + din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, gqbar, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, gqbar, dout, timestep); + } + else if (BoundaryCondition == 5){ + ScaLBL_Comm->D3Q19_Reflection_BC_z(gqbar); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(gqbar); + } + ScaLBL_D3Q19_AAeven_FreeLeeModel_SingleFluid_BGK(gqbar, Velocity, Pressure, tau, rho0, Fx, Fy, Fz, + 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + //************************************************************************ + PROFILE_STOP("Update"); + } + PROFILE_STOP("Loop"); + PROFILE_SAVE("lbpm_color_simulator",1); + //************************************************************************ + stoptime = MPI_Wtime(); + if (rank==0) printf("-------------------------------------------------------------------\n"); + // Compute the walltime per timestep + cputime = (stoptime - starttime)/timestep; + // Performance obtained from each node + double MLUPS = double(Np)/cputime/1000000; + + if (rank==0) printf("********************************************************\n"); + if (rank==0) printf("CPU time = %f \n", cputime); + if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; + if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + if (rank==0) printf("********************************************************\n"); + + // ************************************************************************ +} + +void ScaLBL_FreeLeeModel::WriteDebug_TwoFluid(){ // Copy back final phase indicator field and convert to regular layout DoubleArray PhaseData(Nxh,Nyh,Nzh); //ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); @@ -765,3 +1003,37 @@ void ScaLBL_FreeLeeModel::WriteDebug(){ fclose(CGZ_FILE); */ } + +void ScaLBL_FreeLeeModel::WriteDebug_SingleFluid(){ + + DoubleArray PhaseField(Nx,Ny,Nz); + + // Copy back final phase indicator field and convert to regular layout + ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField); + FILE *PFILE; + sprintf(LocalRankFilename,"Pressure.%05i.raw",rank); + PFILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,PFILE); + fclose(PFILE); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField); + FILE *VELX_FILE; + sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank); + VELX_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,VELX_FILE); + fclose(VELX_FILE); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField); + FILE *VELY_FILE; + sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank); + VELY_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,VELY_FILE); + fclose(VELY_FILE); + + ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField); + FILE *VELZ_FILE; + sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank); + VELZ_FILE = fopen(LocalRankFilename,"wb"); + fwrite(PhaseField.data(),8,N,VELZ_FILE); + fclose(VELZ_FILE); +} diff --git a/models/FreeLeeModel.h b/models/FreeLeeModel.h index 1b78792a..1e372f50 100644 --- a/models/FreeLeeModel.h +++ b/models/FreeLeeModel.h @@ -26,15 +26,20 @@ public: void ReadParams(std::shared_ptr db0); void SetDomain(); void ReadInput(); - void Create(); - void Initialize(); - void Run(); - void WriteDebug(); + void Create_TwoFluid(); + void Initialize_TwoFluid(); + void Run_TwoFluid(); + void WriteDebug_TwoFluid(); + void Create_SingleFluid(); + void Initialize_SingleFluid(); + void Run_SingleFluid(); + void WriteDebug_SingleFluid(); bool Restart,pBC; int timestep,timestepMax; int BoundaryCondition; double tauA,tauB,rhoA,rhoB; + double tau, rho0;//only for single-fluid Lee model double tauM;//relaxation time for phase field (or mass) double W,gamma,kappa,beta; double Fx,Fy,Fz,flux; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 63086219..8df4e6bd 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -7,6 +7,7 @@ ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator ) ADD_LBPM_EXECUTABLE( lbpm_greyscaleColor_simulator ) ADD_LBPM_EXECUTABLE( lbpm_electrokinetic_SingleFluid_simulator ) ADD_LBPM_EXECUTABLE( lbpm_freelee_simulator ) +ADD_LBPM_EXECUTABLE( lbpm_freelee_SingleFluidBGK_simulator ) #ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator ) #ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator ) ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator ) diff --git a/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp new file mode 100644 index 00000000..19d99b9c --- /dev/null +++ b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp @@ -0,0 +1,70 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "common/Utilities.h" +#include "models/FreeLeeModel.h" + +//******************************************************************* +// Implementation of Free-Energy Two-Phase LBM (Lee model) +//******************************************************************* + +int main( int argc, char **argv ) +{ + + // Initialize + Utilities::startup( argc, argv ); + + // Load the input database + auto db = std::make_shared( argv[1] ); + + { // Limit scope so variables that contain communicators will free before MPI_Finialize + + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + + if (rank == 0){ + printf("********************************************************\n"); + printf("Running Single-Fluid Solver based on Lee LBM \n"); + printf("********************************************************\n"); + } + // Initialize compute device + int device=ScaLBL_SetDevice(rank); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); + + PROFILE_ENABLE(1); + //PROFILE_ENABLE_TRACE(); + //PROFILE_ENABLE_MEMORY(); + PROFILE_SYNCHRONIZE(); + PROFILE_START("Main"); + Utilities::setErrorHandlers(); + + auto filename = argv[1]; + ScaLBL_FreeLeeModel LeeModel( rank,nprocs,comm ); + LeeModel.ReadParams( filename ); + LeeModel.SetDomain(); + LeeModel.ReadInput(); + LeeModel.Create_SingleFluid(); + LeeModel.Initialize_SingleFluid(); + LeeModel.Run_SingleFluid(); + LeeModel.WriteDebug_SingleFluid(); + + PROFILE_STOP("Main"); + auto file = db->getWithDefault( "TimerFile", "lbpm_freelee_SingleFluidBGK_simulator" ); + auto level = db->getWithDefault( "TimerLevel", 1 ); + PROFILE_SAVE( file,level ); + // **************************************************** + + + } // Limit scope so variables that contain communicators will free before MPI_Finialize + + Utilities::shutdown(); + return 0; +} diff --git a/tests/lbpm_freelee_simulator.cpp b/tests/lbpm_freelee_simulator.cpp index 3e9c372a..3663c4e9 100644 --- a/tests/lbpm_freelee_simulator.cpp +++ b/tests/lbpm_freelee_simulator.cpp @@ -51,10 +51,10 @@ int main( int argc, char **argv ) LeeModel.ReadParams( filename ); LeeModel.SetDomain(); LeeModel.ReadInput(); - LeeModel.Create(); - LeeModel.Initialize(); - LeeModel.Run(); - LeeModel.WriteDebug(); + LeeModel.Create_TwoFluid(); + LeeModel.Initialize_TwoFluid(); + LeeModel.Run_TwoFluid(); + LeeModel.WriteDebug_TwoFluid(); PROFILE_STOP("Main"); auto file = db->getWithDefault( "TimerFile", "lbpm_freelee_simulator" );