From fee3d9eadc65e61cc07ba03d597785d8abf68e97 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Tue, 11 May 2021 22:25:05 -0400 Subject: [PATCH] save unfinished work; to be built and tested --- common/ScaLBL.cpp | 167 +++++++++++++++++++++++++++++- common/ScaLBL.h | 13 ++- cuda/D3Q7BC.cu | 67 +++++++++++++ models/PoissonSolver.cpp | 2 - models/StokesModel.cpp | 212 ++++++++++++++++++++++++++++++++++++++- models/StokesModel.h | 3 + 6 files changed, 457 insertions(+), 7 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 4726bae6..ae94dd34 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) +void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC=false) { int idx,i,j,k; @@ -1051,6 +1051,23 @@ void ScaLBL_Communicator::SetupBounceBackList(IntArray &Map, signed char *id, in int *bb_interactions_tmp = new int [local_count]; ScaLBL_AllocateDeviceMemory((void **) &bb_dist, sizeof(int)*local_count); ScaLBL_AllocateDeviceMemory((void **) &bb_interactions, sizeof(int)*local_count); + int *fluid_boundary_tmp; + double *lattice_weight_tmp; + float *lattice_cx_tmp; + float *lattice_cy_tmp; + float *lattice_cz_tmp; + if(SlippingVelBC==true){ + fluid_boundary_tmp = new int [local_count]; + lattice_weight_tmp = new double [local_count]; + lattice_cx_tmp = new float [local_count]; + lattice_cy_tmp = new float [local_count]; + lattice_cz_tmp = new float [local_count]; + ScaLBL_AllocateDeviceMemory((void **) &fluid_boundary, sizeof(int)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_weight, sizeof(double)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_cx, sizeof(float)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_cy, sizeof(float)*local_count); + ScaLBL_AllocateDeviceMemory((void **) &lattice_cz, sizeof(float)*local_count); + } local_count=0; for (k=1;k>>(dist, zeta_potential, ElectricField, SolidGrad, + epsilon_LB, tau, rho0, den_scale, h, time_conv, + BounceBackDist_list, BounceBackSolid_list, FluidBoundary_list, + lattice_weight, lattice_cx, lattice_cy, lattice_cz, + count, Np); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_Solid_SlippingVelocityBC_D3Q19 (kernel): %s \n",cudaGetErrorString(err)); + } +} + extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){ int GRID = count / 512 + 1; dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vin, count, Np); diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 25a31600..7f4d7cdc 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -98,9 +98,7 @@ void ScaLBL_Poisson::ReadParams(string filename){ h = domain_db->getScalar( "voxel_length" ); } - //Re-calcualte model parameters if user updates input - epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)] epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity if (rank==0) printf("***********************************************************************************\n"); diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp index fe6b0c92..a9c5e040 100644 --- a/models/StokesModel.cpp +++ b/models/StokesModel.cpp @@ -8,6 +8,7 @@ ScaLBL_StokesModel::ScaLBL_StokesModel(int RANK, int NP, const Utilities::MPI& COMM): rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0), Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),mu(0),h(0),nu_phys(0),rho_phys(0),rho0(0),den_scale(0),time_conv(0),tolerance(0), +epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),UseSlippingVelBC(0), Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM) { @@ -38,6 +39,12 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){ tolerance = 1.0e-8; Fx = Fy = 0.0; Fz = 1.0e-5; + //Stokes solver also needs the following parameters for slipping velocity BC + epsilon0 = 8.85e-12;//electric permittivity of vaccum; unit:[C/(V*m)] + epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)] + epsilonR = 78.4;//default dielectric constant of water + epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity + UseSlippingVelBC = false; //--------------------------------------------------------------------------// // Read domain parameters @@ -85,12 +92,18 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){ if (stokes_db->keyExists( "flux" )){ flux = stokes_db->getScalar( "flux" ); } + if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){ + UseSlippingVelBC = stokes_db->getScalar( "UseElectroosmoticVelocityBC" ); + } + if (electric_db->keyExists( "epsilonR" )){ + epsilonR = electric_db->getScalar( "epsilonR" ); + } // Re-calculate model parameters due to parameter read mu=(tau-0.5)/3.0; time_conv = (h*h*1.0e-12)*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt] den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density - + epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity } void ScaLBL_StokesModel::ReadParams(string filename){ @@ -100,7 +113,6 @@ void ScaLBL_StokesModel::ReadParams(string filename){ db = std::make_shared( filename ); domain_db = db->getDatabase( "Domain" ); stokes_db = db->getDatabase( "Stokes" ); - //---------------------- Default model parameters --------------------------// rho_phys = 1000.0; //by default use water density; unit [kg/m^3] @@ -114,6 +126,12 @@ void ScaLBL_StokesModel::ReadParams(string filename){ tolerance = 1.0e-8; Fx = Fy = 0.0; Fz = 1.0e-5; + //Stokes solver also needs the following parameters for slipping velocity BC + epsilon0 = 8.85e-12;//electric permittivity of vaccum; unit:[C/(V*m)] + epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)] + epsilonR = 78.4;//default dielectric constant of water + epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity + UseSlippingVelBC = false; //--------------------------------------------------------------------------// // Read domain parameters @@ -161,12 +179,18 @@ void ScaLBL_StokesModel::ReadParams(string filename){ if (stokes_db->keyExists( "flux" )){ flux = stokes_db->getScalar( "flux" ); } + if (stokes_db->keyExists( "UseElectroosmoticVelocityBC" )){ + UseSlippingVelBC = stokes_db->getScalar( "UseElectroosmoticVelocityBC" ); + } + if (electric_db->keyExists( "epsilonR" )){ + epsilonR = electric_db->getScalar( "epsilonR" ); + } // Re-calculate model parameters due to parameter read mu=(tau-0.5)/3.0; time_conv = (h*h*1.0e-12)*mu/nu_phys;//time conversion factor from physical to LB unit; [sec/lt] den_scale = rho_phys/rho0*(h*h*h*1.0e-18);//scale factor for density - + epsilon_LB = epsilon0_LB*epsilonR;//electric permittivity } void ScaLBL_StokesModel::SetDomain(){ @@ -258,6 +282,157 @@ void ScaLBL_StokesModel::ReadInput(){ if (rank == 0) cout << " Domain set." << endl; } +void ScaLBL_IonModel::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" ); + + NLABELS=LabelList.size(); + if (NLABELS != AffinityList.size()){ + ERROR("Error: LB Stokes Solver: SolidLabels and ZetaPotentialSolidList must be the same length! \n"); + } + + double label_count[NLABELS]; + double label_count_global[NLABELS]; + + for (size_t idx=0; idxid[n]; + AFFINITY=0.f; + // Assign the affinity from the paired list + for (unsigned int idx=0; idx < NLABELS; idx++){ + if (VALUE == LabelList[idx]){ + AFFINITY=AffinityList[idx];//no need to convert unit for zeta potential (i.e. volt) + label_count[idx] += 1.0; + idx = NLABELS; + } + } + zeta_potential_solid[n] = AFFINITY; + } + } + } + + for (size_t idx=0; idxComm.sumReduce( label_count[idx]); + + if (rank==0){ + printf("LB Stokes Solver: number of solid labels: %lu \n",NLABELS); + for (unsigned int idx=0; idxid[nn]); + //Since the solid unit normal vector is wanted, treat + //wet node as 0.0 and solid node as 1.0 + GWNS = (GWNS>0.0) ? 0.0:1.0; + phi_x += GWNS*weight*vec_x; + phi_y += GWNS*weight*vec_y; + phi_z += GWNS*weight*vec_z; + } + } + } + solid_grad[idx+0*Np] = phi_x; + solid_grad[idx+1*Np] = phi_y; + solid_grad[idx+2*Np] = phi_z; + } + } + } + } +} + void ScaLBL_StokesModel::Create(){ /* * This function creates the variables needed to run a LBM @@ -301,6 +476,26 @@ void ScaLBL_StokesModel::Create(){ ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); comm.barrier(); + if (UseSlippingVelBC==true){ + ScaLBL_Comm->SetupBounceBackList(Map, Mask->id.data(), Np,1); + comm.barrier(); + + //For slipping velocity BC, need zeta potential and solid unit normal vector + ScaLBL_AllocateDeviceMemory((void **) &ZetaPotentialSolid, sizeof(double)*Nx*Ny*Nz); + ScaLBL_AllocateDeviceMemory((void **) &SolidGrad, sizeof(double)*Np); //unit normal vector of solid nodes + + double *ZetaPotentialSolid_host; + ZetaPotentialSolid_host = new double[Nx*Ny*Nz]; + AssignZetaPotentialSolid(ZetaPotentialSolid_host); + double *SolidGrad_host; + SolidGrad_host = new double[3*Np]; + AssignSolidGrad(SolidGrad_host); + ScaLBL_CopyToDevice(ZetaPotentialSolid, ZetaPotentialSolid_host, Nx*Ny*Nz*sizeof(double)); + ScaLBL_CopyToDevice(SolidGrad, SolidGrad_host, 3*Np*sizeof(double)); + ScaLBL_Comm->Barrier(); + delete [] ZetaPotentialSolid_host; + delete [] SolidGrad_host; + } } void ScaLBL_StokesModel::Initialize(){ @@ -324,6 +519,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){ timestep = 0; while (timestep < timestepMax) { //************************************************************************/ + //**************ODD TIMESTEP*************// timestep++; 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, @@ -344,8 +540,14 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){ } 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); + + if (UseSlippingVelBC==true){ + ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad, + epslion_LB, 1.0/rlx_setA, rho0, den_scale, h, time_conv); + } ScaLBL_Comm->Barrier(); comm.barrier(); + //**************EVEN TIMESTEP*************// timestep++; 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, @@ -366,6 +568,10 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){ } 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); + if (UseSlippingVelBC==true){ + ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(fq, ZetaPotentialSolid, ElectricField, SolidGrad, + epslion_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 31784a4f..86bdc8b6 100644 --- a/models/StokesModel.h +++ b/models/StokesModel.h @@ -51,6 +51,8 @@ public: double time_conv; double h;//image resolution double den_scale;//scale factor for density + double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;//Stokes solver also needs this for slipping velocity BC + bool UseSlippingVelBC; int Nx,Ny,Nz,N,Np; int rank,nprocx,nprocy,nprocz,nprocs; @@ -70,6 +72,7 @@ public: double *fq; double *Velocity; double *Pressure; + double *ZetaPotentialSolid; //Minkowski Morphology; DoubleArray Velocity_x;