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..f51c8323 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,9 @@ 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; idx = Map[n]; psi = Psi[idx]; @@ -175,8 +176,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 +189,9 @@ 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; 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/cuda/Poisson.cu b/cuda/Poisson.cu index 84a78330..dd0cb4bc 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,9 @@ __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 +307,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){ diff --git a/hip/Poisson.cu b/hip/Poisson.cu index 34975f58..de6cc8eb 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,9 @@ __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 +307,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){ 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..ca97054d 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, @@ -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/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/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