From cd2bcfba85d37b27bfc31fd6d7ea3338a078dc8a Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Sat, 15 May 2021 02:27:02 -0400 Subject: [PATCH] build pass; TODO: verify the slippingBC model --- common/ScaLBL.cpp | 6 +++--- cuda/D3Q7BC.cu | 7 ++++--- models/StokesModel.cpp | 30 ++++++++++++++++-------------- models/StokesModel.h | 3 +++ 4 files changed, 26 insertions(+), 20 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index ae94dd34..ec13efc3 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -973,7 +973,7 @@ int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborLis } -void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC=false) +void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC) { int idx,i,j,k; @@ -1362,10 +1362,10 @@ void ScaLBL_Communicator::SolidNeumannD3Q7(double *fq, double *BoundaryValue){ } void ScaLBL_Communicator::SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_potential, double *ElectricField, double *SolidGrad, - double epslion_LB, double tau, double rho0, double den_scale,double h, double time_conv){ + double epsilon_LB, double tau, double rho0, double den_scale,double h, double time_conv){ // fq is a D3Q19 distribution // BoundaryValues is a list of values to assign at bounce-back solid sites - ScaLBL_Solid_SlippingVelocityBC_D3Q19(fq,zeta_potential,Electricfield,SolidGrad,epsilon_LB,tau,rho0,den_scale,h,time_conv, + ScaLBL_Solid_SlippingVelocityBC_D3Q19(fq,zeta_potential,ElectricField,SolidGrad,epsilon_LB,tau,rho0,den_scale,h,time_conv, bb_dist,bb_interactions,fluid_boundary,lattice_weight,lattice_cx,lattice_cy,lattice_cz,n_bb_d3q19,N); } diff --git a/cuda/D3Q7BC.cu b/cuda/D3Q7BC.cu index 53b9a810..d5f8091b 100644 --- a/cuda/D3Q7BC.cu +++ b/cuda/D3Q7BC.cu @@ -49,6 +49,7 @@ __global__ void dvc_ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double * double value_b,value_q; double Ex,Ey,Ez; double Etx,Ety,Etz;//tangential part of electric field + double E_mag_normal; double nsx,nsy,nsz;//unit normal solid gradient double ubx,uby,ubz;//slipping velocity at fluid boundary nodes float cx,cy,cz;//lattice velocity (D3Q19) @@ -75,9 +76,9 @@ __global__ void dvc_ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double * Etx = Ex - E_mag_normal*nsx; Ety = Ey - E_mag_normal*nsy; Etz = Ez - E_mag_normal*nsz; - ubx = -eplison_LB*value_b*Etx/(nu_LB*rho0)*time_conv*time_conv/h/h/den_scale; - uby = -eplison_LB*value_b*Ety/(nu_LB*rho0)*time_conv*time_conv/h/h/den_scale; - ubz = -eplison_LB*value_b*Etz/(nu_LB*rho0)*time_conv*time_conv/h/h/den_scale; + ubx = -epsilon_LB*value_b*Etx/(nu_LB*rho0)*time_conv*time_conv/h/h/den_scale; + uby = -epsilon_LB*value_b*Ety/(nu_LB*rho0)*time_conv*time_conv/h/h/den_scale; + ubz = -epsilon_LB*value_b*Etz/(nu_LB*rho0)*time_conv*time_conv/h/h/den_scale; //compute bounce-back distribution LB_weight = lattice_weight[idx]; diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp index a9c5e040..7d70928d 100644 --- a/models/StokesModel.cpp +++ b/models/StokesModel.cpp @@ -95,8 +95,8 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){ if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){ UseSlippingVelBC = stokes_db->getScalar( "UseElectroosmoticVelocityBC" ); } - if (electric_db->keyExists( "epsilonR" )){ - epsilonR = electric_db->getScalar( "epsilonR" ); + if (stokes_db->keyExists( "epsilonR" )){ + epsilonR = stokes_db->getScalar( "epsilonR" ); } // Re-calculate model parameters due to parameter read @@ -182,8 +182,8 @@ void ScaLBL_StokesModel::ReadParams(string filename){ if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){ UseSlippingVelBC = stokes_db->getScalar( "UseElectroosmoticVelocityBC" ); } - if (electric_db->keyExists( "epsilonR" )){ - epsilonR = electric_db->getScalar( "epsilonR" ); + if (stokes_db->keyExists( "epsilonR" )){ + epsilonR = stokes_db->getScalar( "epsilonR" ); } // Re-calculate model parameters due to parameter read @@ -282,14 +282,14 @@ void ScaLBL_StokesModel::ReadInput(){ if (rank == 0) cout << " Domain set." << endl; } -void ScaLBL_IonModel::AssignZetaPotentialSolid(double *zeta_potential_solid) +void ScaLBL_StokesModel::AssignZetaPotentialSolid(double *zeta_potential_solid) { size_t NLABELS=0; signed char VALUE=0; double AFFINITY=0.f; - auto LabelList = ion_db->getVector( "SolidLabels" ); - auto AffinityList = ion_db->getVector( "ZetaPotentialSolidList" ); + auto LabelList = stokes_db->getVector( "SolidLabels" ); + auto AffinityList = stokes_db->getVector( "ZetaPotentialSolidList" ); NLABELS=LabelList.size(); if (NLABELS != AffinityList.size()){ @@ -335,9 +335,8 @@ void ScaLBL_IonModel::AssignZetaPotentialSolid(double *zeta_potential_solid) } } -void ScaLBL_IonModel::AssignSolidGrad(double *solid_grad) +void ScaLBL_StokesModel::AssignSolidGrad(double *solid_grad) { - //TODO need to normalize the computed solid grad!!! double *Dst; Dst = new double [3*3*3]; for (int kk=0; kk<3; kk++){ @@ -424,9 +423,12 @@ void ScaLBL_IonModel::AssignSolidGrad(double *solid_grad) } } } - solid_grad[idx+0*Np] = phi_x; - solid_grad[idx+1*Np] = phi_y; - solid_grad[idx+2*Np] = phi_z; + //solid_grad normalization + double phi_mag=sqrt(phi_x*phi_x+phi_y*phi_y+phi_z*phi_z); + if (phi_mag==0.0) phi_mag=1.0; + solid_grad[idx+0*Np] = phi_x/phi_mag; + solid_grad[idx+1*Np] = phi_y/phi_mag; + solid_grad[idx+2*Np] = phi_z/phi_mag; } } } @@ -543,7 +545,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){ if (UseSlippingVelBC==true){ ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad, - epslion_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv); + epsilon_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv); } ScaLBL_Comm->Barrier(); comm.barrier(); @@ -570,7 +572,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){ 0, ScaLBL_Comm->LastExterior(), Np); if (UseSlippingVelBC==true){ ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad, - epslion_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv); + epsilon_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv); } ScaLBL_Comm->Barrier(); comm.barrier(); //************************************************************************/ diff --git a/models/StokesModel.h b/models/StokesModel.h index 86bdc8b6..1e80f269 100644 --- a/models/StokesModel.h +++ b/models/StokesModel.h @@ -73,6 +73,7 @@ public: double *Velocity; double *Pressure; double *ZetaPotentialSolid; + double *SolidGrad; //Minkowski Morphology; DoubleArray Velocity_x; @@ -91,5 +92,7 @@ private: void LoadParams(std::shared_ptr db0); void Velocity_LB_to_Phys(DoubleArray &Vel_reg); vector computeElectricForceAvg(double *ChargeDensity, double *ElectricField); + void AssignSolidGrad(double *solid_grad); + void AssignZetaPotentialSolid(double *zeta_potential_solid); }; #endif