From ee527a42cca9f216207a111f4a27d0a733feb8c9 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 18 Jan 2022 16:46:09 +1100 Subject: [PATCH 1/5] update Poisson and Stoke solvers to include slipping velocity BC; to be built --- common/ScaLBL.h | 13 +++++--- cpu/Poisson.cpp | 20 +++++++----- cpu/Stokes.cpp | 31 ++++++++++++------- models/PoissonSolver.cpp | 18 +++++------ models/PoissonSolver.h | 2 +- models/StokesModel.cpp | 14 ++++----- ...m_electrokinetic_SingleFluid_simulator.cpp | 2 +- 7 files changed, 58 insertions(+), 42 deletions(-) diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 2be0121e..cb1da09d 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -296,7 +296,7 @@ extern "C" void ScaLBL_D3Q7_Ion_ChargeDensity(double *Den, double *ChargeDensity * @param finish - lattice node to finish loop * @param Np - size of local sub-domain (derived from Domain structure) */ -extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB, +extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np); /** @@ -307,19 +307,22 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList,int *Map, double *di * @param Psi - * @param ElectricField - electric field * @param tau - relaxation time -* @param epsilon_LB - +* @param epsilon_LB - dielectric constant of medium +* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve * @param start - lattice node to start loop * @param finish - lattice node to finish loop * @param Np - size of local sub-domain (derived from Domain structure) */ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, - double epsilon_LB, int start, int finish, int Np); + double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np); /** * \brief Poisson-Boltzmann solver / solve electric potential based on AA odd access pattern for D3Q7 * @param neighborList - neighbors based on D3Q19 lattice structure * @param Map - mapping between sparse and dense representations * @param dist - D3Q7 distributions * @param Psi - +* @param epsilon_LB - dielectric constant of medium +* @param UseSlippingVelBC - flag indicating the use of Helmholtz-Smoluchowski slipping velocity equation when EDL is too small to resolve * @param start - lattice node to start loop * @param finish - lattice node to finish loop * @param Np - size of local sub-domain (derived from Domain structure) @@ -350,10 +353,10 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in // LBM Stokes Model (adapted from MRT model) extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, - double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np); + double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np); extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, - double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np); + double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np); extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np); diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp index 92b89eae..ae3a9eb6 100644 --- a/cpu/Poisson.cpp +++ b/cpu/Poisson.cpp @@ -95,7 +95,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential( extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, - double tau, double epsilon_LB, + double tau, double epsilon_LB, bool UseSlippingVelBC, int start, int finish, int Np) { int n; double psi; //electric potential @@ -109,8 +109,11 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, for (n = start; n < finish; n++) { //Load data - rho_e = Den_charge[n]; - rho_e = rho_e / epsilon_LB; + //rho_e = Den_charge[n]; + //rho_e = rho_e / epsilon_LB; + //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral + //and thus the net space charge density is zero. + rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; idx = Map[n]; psi = Psi[idx]; @@ -175,8 +178,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, - double epsilon_LB, int start, - int finish, int Np) { + double epsilon_LB, bool UseSlippingVelBC, + int start, int finish, int Np) { int n; double psi; //electric potential double Ex, Ey, Ez; //electric field @@ -188,8 +191,11 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, for (n = start; n < finish; n++) { //Load data - rho_e = Den_charge[n]; - rho_e = rho_e / epsilon_LB; + //rho_e = Den_charge[n]; + //rho_e = rho_e / epsilon_LB; + //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral + //and thus the net space charge density is zero. + rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; idx = Map[n]; psi = Psi[idx]; diff --git a/cpu/Stokes.cpp b/cpu/Stokes.cpp index 5955bb8a..ce1ff7f8 100644 --- a/cpu/Stokes.cpp +++ b/cpu/Stokes.cpp @@ -4,7 +4,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT( double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, - double time_conv, int start, int finish, int Np) { + double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) { double fq; // conserved momemnts double rho, jx, jy, jz; @@ -38,13 +38,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT( Ey = ElectricField[n + 1 * Np]; Ez = ElectricField[n + 2 * Np]; //compute total body force, including input body force (Gx,Gy,Gz) - Fx = - Gx + - rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / - den_scale; //the extra factors at the end necessarily convert unit from phys to LB - Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / + Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; //the extra factors at the end necessarily convert unit from phys to LB + Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / den_scale; - Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / + Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / den_scale; // q=0 @@ -479,7 +477,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT( int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, - double time_conv, int start, int finish, int Np) { + double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) { double fq; // conserved momemnts double rho, jx, jy, jz; @@ -513,12 +511,21 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT( Ex = ElectricField[n + 0 * Np]; Ey = ElectricField[n + 1 * Np]; Ez = ElectricField[n + 2 * Np]; + //compute total body force, including input body force (Gx,Gy,Gz) - Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + //Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + // den_scale; //the extra factors at the end necessarily convert unit from phys to LB + //Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / + // den_scale; + //Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / + // den_scale; + //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral + //and body force induced by external efectric field is reduced to slipping velocity BC. + Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; //the extra factors at the end necessarily convert unit from phys to LB + Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / den_scale; - Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / - den_scale; - Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / + Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / den_scale; // q=0 diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 8658ac07..e3dccfd2 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -523,7 +523,7 @@ void ScaLBL_Poisson::Initialize(double time_conv_from_Study){ //} } -void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){ +void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study){ //.......create and start timer............ //double starttime,stoptime,cputime; @@ -537,13 +537,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){ // *************ODD TIMESTEP*************// timestep++; SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential - SolvePoissonAAodd(ChargeDensity);//perform collision + SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); // *************EVEN TIMESTEP*************// timestep++; SolveElectricPotentialAAeven(timestep_from_Study);//update electric potential - SolvePoissonAAeven(ChargeDensity);//perform collision + SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision ScaLBL_Comm->Barrier(); comm.barrier(); //************************************************************************/ @@ -724,9 +724,9 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np); } -void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){ - ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np); +void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC){ + ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); //TODO: perhaps add another ScaLBL_Comm routine to update Psi values on solid boundary nodes. //something like: //ScaLBL_Comm->SolidDirichletBoundaryUpdates(Psi, Psi_BCLabel, timestep); @@ -739,9 +739,9 @@ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){ //} } -void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity){ - ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np); +void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC){ + ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, UseSlippingVelBC, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel); //if (BoundaryConditionSolid==1){ // ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi); diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index bb021218..e815951c 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -34,7 +34,7 @@ public: void ReadInput(); void Create(); void Initialize(double time_conv_from_Study); - void Run(double *ChargeDensity, int timestep_from_Study); + void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study); void getElectricPotential(DoubleArray &ReturnValues); void getElectricPotential_debug(int timestep); void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y, diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp index 1e395ba9..2098be5f 100644 --- a/models/StokesModel.cpp +++ b/models/StokesModel.cpp @@ -593,7 +593,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL ScaLBL_D3Q19_AAodd_StokesMRT( NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, - rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv, + rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE // Set boundary conditions @@ -610,8 +610,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, } ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, - Fz, rho0, den_scale, h, time_conv, 0, - ScaLBL_Comm->LastExterior(), Np); + Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC, + 0, ScaLBL_Comm->LastExterior(), Np); if (UseSlippingVelBC == true) { ScaLBL_Comm->SolidSlippingVelocityBCD3Q19( @@ -626,8 +626,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL ScaLBL_D3Q19_AAeven_StokesMRT( fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, - Fy, Fz, rho0, den_scale, h, time_conv, ScaLBL_Comm->FirstInterior(), - ScaLBL_Comm->LastInterior(), Np); + Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE // Set boundary conditions if (BoundaryCondition == 3) { @@ -643,8 +643,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, } ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, - Fz, rho0, den_scale, h, time_conv, 0, - ScaLBL_Comm->LastExterior(), Np); + Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC, + 0, ScaLBL_Comm->LastExterior(), Np); if (UseSlippingVelBC == true) { ScaLBL_Comm->SolidSlippingVelocityBCD3Q19( fq, ZetaPotentialSolid, ElectricField, SolidGrad, epsilon_LB, diff --git a/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp b/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp index caeef89a..e814f915 100644 --- a/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp +++ b/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp @@ -94,7 +94,7 @@ int main(int argc, char **argv) while (timestep < Study.timestepMax){ timestep++; - PoissonSolver.Run(IonModel.ChargeDensity,timestep);//solve Poisson equtaion to get steady-state electrical potental + PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,timestep);//solve Poisson equtaion to get steady-state electrical potental StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential From 169a102f6c48378114fda0cd3093d978bae338c9 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 18 Jan 2022 01:15:35 -0500 Subject: [PATCH 2/5] fix a few typos and bug and build passed --- models/PoissonSolver.h | 4 ++-- tests/TestNernstPlanck.cpp | 2 +- tests/TestPNP_Stokes.cpp | 2 +- tests/TestPoissonSolver.cpp | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index e815951c..ca97054d 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -111,8 +111,8 @@ private: void SolveElectricPotentialAAodd(int timestep_from_Study); void SolveElectricPotentialAAeven(int timestep_from_Study); //void SolveElectricField(); - void SolvePoissonAAodd(double *ChargeDensity); - void SolvePoissonAAeven(double *ChargeDensity); + void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC); + void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC); void getConvergenceLog(int timestep,double error); double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step); diff --git a/tests/TestNernstPlanck.cpp b/tests/TestNernstPlanck.cpp index f0f82e52..b9492d9d 100644 --- a/tests/TestNernstPlanck.cpp +++ b/tests/TestNernstPlanck.cpp @@ -74,7 +74,7 @@ int main(int argc, char **argv) while (timestep < Study.timestepMax && error > Study.tolerance){ timestep++; - PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental + PoissonSolver.Run(IonModel.ChargeDensity,false,0);//solve Poisson equtaion to get steady-state electrical potental IonModel.Run(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField); //solve for ion transport and electric potential timestep++;//AA operations diff --git a/tests/TestPNP_Stokes.cpp b/tests/TestPNP_Stokes.cpp index 65b796f7..19324f8f 100644 --- a/tests/TestPNP_Stokes.cpp +++ b/tests/TestPNP_Stokes.cpp @@ -94,7 +94,7 @@ int main(int argc, char **argv) while (timestep < Study.timestepMax && error > Study.tolerance){ timestep++; - PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental + PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,0);//solve Poisson equtaion to get steady-state electrical potental StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential diff --git a/tests/TestPoissonSolver.cpp b/tests/TestPoissonSolver.cpp index 38f242c6..7d39889a 100644 --- a/tests/TestPoissonSolver.cpp +++ b/tests/TestPoissonSolver.cpp @@ -69,7 +69,7 @@ int main(int argc, char **argv) int timeSave = int(PoissonSolver.TestPeriodicSaveInterval/PoissonSolver.TestPeriodicTimeConv); while (timestep Date: Thu, 27 Jan 2022 14:44:41 +1100 Subject: [PATCH 3/5] udpate cuda with corrected slipping vel BC --- cuda/Poisson.cu | 26 ++++++++++++++++---------- cuda/Stokes.cu | 30 ++++++++++++++++++------------ 2 files changed, 34 insertions(+), 22 deletions(-) diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index 84a78330..3d637e7d 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -104,7 +104,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub } } -__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){ +__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){ int n; double psi;//electric potential @@ -122,8 +122,11 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub if (n>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np); + dvc_ScaLBL_D3Q7_AAodd_Poisson<<>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -305,10 +311,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q7_AAeven_Poisson<<>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np); + dvc_ScaLBL_D3Q7_AAeven_Poisson<<>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ diff --git a/cuda/Stokes.cu b/cuda/Stokes.cu index d091b0b4..e6a2055a 100644 --- a/cuda/Stokes.cu +++ b/cuda/Stokes.cu @@ -5,7 +5,7 @@ #define NBLOCKS 1024 #define NTHREADS 256 -__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){ +__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){ int n; double fq; @@ -46,9 +46,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis Ey = ElectricField[n+1*Np]; Ez = ElectricField[n+2*Np]; //compute total body force, including input body force (Gx,Gy,Gz) - Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; - Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; - Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; + Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; //the extra factors at the end necessarily convert unit from phys to LB + Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; + Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; // q=0 fq = dist[n]; @@ -510,7 +513,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis } } -__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){ +__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){ int n; double fq; @@ -550,9 +553,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit Ey = ElectricField[n+1*Np]; Ez = ElectricField[n+2*Np]; //compute total body force, including input body force (Gx,Gy,Gz) - Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB - Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; - Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; + Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; //the extra factors at the end necessarily convert unit from phys to LB + Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; + Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; // q=0 fq = dist[n]; @@ -969,10 +975,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit } } -extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np); + dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ @@ -981,10 +987,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np); + dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np); cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ From 02932e26c56d0d025b740518b86d78fdd56c4162 Mon Sep 17 00:00:00 2001 From: Li Rex Date: Thu, 27 Jan 2022 16:49:05 +1100 Subject: [PATCH 4/5] update slipping vel bc also in HIP version --- hip/Poisson.cu | 26 ++++++++++++++++---------- hip/Stokes.cu | 30 ++++++++++++++++++------------ 2 files changed, 34 insertions(+), 22 deletions(-) diff --git a/hip/Poisson.cu b/hip/Poisson.cu index 34975f58..e869a17b 100644 --- a/hip/Poisson.cu +++ b/hip/Poisson.cu @@ -104,7 +104,7 @@ __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub } } -__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){ +__global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){ int n; double psi;//electric potential @@ -122,8 +122,11 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub if (n>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np); + dvc_ScaLBL_D3Q7_AAodd_Poisson<<>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ @@ -305,10 +311,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q7_AAeven_Poisson<<>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np); + dvc_ScaLBL_D3Q7_AAeven_Poisson<<>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ diff --git a/hip/Stokes.cu b/hip/Stokes.cu index a6a05fba..79ea9c4a 100644 --- a/hip/Stokes.cu +++ b/hip/Stokes.cu @@ -6,7 +6,7 @@ #define NTHREADS 256 -__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){ +__global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){ int n; double fq; @@ -47,9 +47,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis Ey = ElectricField[n+1*Np]; Ez = ElectricField[n+2*Np]; //compute total body force, including input body force (Gx,Gy,Gz) - Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; - Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; - Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; + Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; //the extra factors at the end necessarily convert unit from phys to LB + Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; + Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; // q=0 fq = dist[n]; @@ -511,7 +514,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis } } -__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){ +__global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC,int start, int finish, int Np){ int n; double fq; @@ -551,9 +554,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit Ey = ElectricField[n+1*Np]; Ez = ElectricField[n+2*Np]; //compute total body force, including input body force (Gx,Gy,Gz) - Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB - Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; - Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale; + Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; //the extra factors at the end necessarily convert unit from phys to LB + Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; + Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) / + den_scale; // q=0 fq = dist[n]; @@ -970,10 +976,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit } } -extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np); + dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ @@ -982,10 +988,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do //cudaProfilerStop(); } -extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){ +extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){ //cudaProfilerStart(); - dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np); + dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ From 71fdedaa2fb0bed93f3d2bcf229429726eeb5a42 Mon Sep 17 00:00:00 2001 From: Li Rex Date: Fri, 28 Jan 2022 16:26:34 +1100 Subject: [PATCH 5/5] clean up code for the updated slipping vel BC --- cpu/Poisson.cpp | 4 ---- cuda/Poisson.cu | 4 ---- hip/Poisson.cu | 4 ---- 3 files changed, 12 deletions(-) diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp index ae3a9eb6..f51c8323 100644 --- a/cpu/Poisson.cpp +++ b/cpu/Poisson.cpp @@ -109,8 +109,6 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, for (n = start; n < finish; n++) { //Load data - //rho_e = Den_charge[n]; - //rho_e = rho_e / epsilon_LB; //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral //and thus the net space charge density is zero. rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; @@ -191,8 +189,6 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, for (n = start; n < finish; n++) { //Load data - //rho_e = Den_charge[n]; - //rho_e = rho_e / epsilon_LB; //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral //and thus the net space charge density is zero. rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB; diff --git a/cuda/Poisson.cu b/cuda/Poisson.cu index 3d637e7d..dd0cb4bc 100644 --- a/cuda/Poisson.cu +++ b/cuda/Poisson.cu @@ -122,8 +122,6 @@ __global__ void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub if (n