diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 6ae49cca..d9b75c43 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1182,7 +1182,7 @@ void ScaLBL_Communicator::SendD3Q7AA(double *Aq, int Component){ // NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2 if (Lock==true){ - ERROR("ScaLBL Error (SendD3Q19): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?"); + ERROR("ScaLBL Error (SendD3Q7): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?"); } else{ Lock=true; @@ -1288,7 +1288,7 @@ void ScaLBL_Communicator::BiSendD3Q7AA(double *Aq, double *Bq){ // NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2 if (Lock==true){ - ERROR("ScaLBL Error (SendD3Q19): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?"); + ERROR("ScaLBL Error (BiSendD3Q7): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?"); } else{ Lock=true; @@ -1415,7 +1415,7 @@ void ScaLBL_Communicator::TriSendD3Q7AA(double *Aq, double *Bq, double *Cq){ // NOTE: the center distribution f0 must NOT be at the start of feven, provide offset to start of f2 if (Lock==true){ - ERROR("ScaLBL Error (SendD3Q19): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?"); + ERROR("ScaLBL Error (TriSendD3Q7): ScaLBL_Communicator is locked -- did you forget to match Send/Recv calls?"); } else{ Lock=true; diff --git a/cpu/Poisson.cpp b/cpu/Poisson.cpp index 355e4223..d24bc12f 100644 --- a/cpu/Poisson.cpp +++ b/cpu/Poisson.cpp @@ -1,65 +1,18 @@ -extern "C" void ScaLBL_D3Q7_AAeven_Poisson(double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz){ + +extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, double *dist, double *Den_charge, double *Psi, double *ElectricField, double rlx, double epsilon_LB,double deltaT, + int start, int finish, int Np){ int n; - // conserved momemnts - double rho,ux,uy,uz,uu; - // non-conserved moments - double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18; + double psi;//electric potential + double Ex,Ey,Ez;//electrical field + double rho_e;//local charge density + double f0,f1,f2,f3,f4,f5,f6; + int nr1,nr2,nr3,nr4,nr5,nr6; for (int n=start; n( filename ); domain_db = db->getDatabase( "Domain" ); - electric_db = db->getDatabase( "Electric" ); + electric_db = db->getDatabase( "Poisson" ); - tau = 1.0; + k2_inv = 4.5;//the inverse of 2nd-rank moment of D3Q7 lattice + deltaT = 0.3;//time step of LB-Poisson equation + tau = 0.5+k2_inv*deltaT; timestepMax = 100000; - tolerance = 1.0e-8; - Fx = Fy = 0.0; - Fz = 1.0e-5; + tolerance = 1.0e-6;//stopping criterion for obtaining steady-state electricla potential + h = 1.0;//resolution; unit: um/lu + epsilon0 = 8.85e-12;//electrical permittivity of vaccum; unit:[C/(V*m)] + epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)] + epsilonR = 78.4;//default dielectric constant for water + epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity + analysis_interval = 1000; - // Color Model parameters + // LB-Poisson Model parameters if (electric_db->keyExists( "timestepMax" )){ timestepMax = electric_db->getScalar( "timestepMax" ); } - + if (electric_db->keyExists( "analysis_interval" )){ + analysis_interval = electric_db->getScalar( "analysis_interval" ); + } + if (electric_db->keyExists( "tolerance" )){ + tolerance = electric_db->getScalar( "tolerance" ); + } + if (electric_db->keyExists( "deltaT" )){ + deltaT = electric_db->getScalar( "deltaT" ); + } + if (electric_db->keyExists( "epsilonR" )){ + epsilonR = electric_db->getScalar( "epsilonR" ); + } // Read domain parameters + if (domain_db->keyExists( "voxel_length" )){//default unit: um/lu + h = domain_db->getScalar( "voxel_length" ); + } if (domain_db->keyExists( "BC" )){ BoundaryCondition = domain_db->getScalar( "BC" ); } + //Re-calcualte model parameters if user updates input + epsilon0_LB = epsilon0*(h*1.0e-6);//unit:[C/(V*lu)] + epsilon_LB = epsilon0_LB*epsilonR;//electrical permittivity + tau = 0.5+k2_inv*deltaT; - mu=(tau-0.5)/3.0; } void ScaLBL_Poisson::SetDomain(){ Dm = std::shared_ptr(new Domain(domain_db,comm)); // full domain for analysis @@ -54,6 +77,7 @@ void ScaLBL_Poisson::SetDomain(){ N = Nx*Ny*Nz; Distance.resize(Nx,Ny,Nz); + Psi_host.resize(Nx,Ny,Nz); for (int i=0; iid[i] = 1; // initialize this way //Averages = std::shared_ptr ( new TwoPhase(Dm) ); // TwoPhase analysis object @@ -158,6 +182,7 @@ void ScaLBL_Poisson::Create(){ ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); ScaLBL_AllocateDeviceMemory((void **) &fq, 7*dist_mem_size); ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np); //........................................................................... // Update GPU data structures if (rank==0) printf ("Setting up device map and neighbor list \n"); @@ -171,54 +196,95 @@ void ScaLBL_Poisson::Initialize(){ /* * This function initializes model */ - if (rank==0) printf ("Initializing distributions \n"); - ScaLBL_D3Q19_Init(fq, Np); + if (rank==0) printf ("Initializing D3Q7 distributions for LB-Poisson solver\n"); + ScaLBL_D3Q7_Poisson_Init(fq, Np); } void ScaLBL_Poisson::Run(double *ChargeDensity){ - double rlx = 1.0/tau; + + //LB-related parameter + double rlx = 1.0/tau; + //.......create and start timer............ double starttime,stoptime,cputime; ScaLBL_DeviceBarrier(); MPI_Barrier(comm); starttime = MPI_Wtime(); - if (rank==0) printf("Beginning AA timesteps, timestepMax = %i \n", timestepMax); - if (rank==0) printf("********************************************************\n"); + if (rank==0) printf("***************************************************************************\n"); + if (rank==0) printf("LB-Poisson Solver: timestepMax = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance); + if (rank==0) printf("***************************************************************************\n"); timestep=0; double error = 1.0; - double flow_rate_previous = 0.0; + double psi_avg_previous = 0.0; while (timestep < timestepMax && error > tolerance) { //************************************************************************/ - timestep++; + // *************ODD TIMESTEP*************// + timestep++; ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FROM NORMAL - ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz); + ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE // Set boundary conditions /* ... */ - ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz); + ScaLBL_D3Q7_AAodd_Poisson(NeighborList, fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + + // *************EVEN TIMESTEP*************// timestep++; ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FORM NORMAL - ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, Fx, Fy, Fz); + ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE // Set boundary conditions /* ... */ - ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, 0, ScaLBL_Comm->LastExterior(), Np, rlx, Fx, Fy, Fz); + ScaLBL_D3Q7_AAeven_Poisson(fq, ChargeDensity, Psi, ElectricField, rlx, epsilon_LB, deltaT, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); //************************************************************************/ + + // Check convergence of steady-state solution + if (timestep%analysis_interval==0){ + + ScaLBL_Comm->RegularLayout(Map,&Psi,Psi_host); + double count_loc=0; + double count; + double psi_avg; + double psi_loc=0.f; + + for (int k=1; k 0){ + psi_loc += Psi_host(i,j,k); + count_loc+=1.0; + } + } + } + } + MPI_Allreduce(&psi_loc,&psi_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm); + MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm); + + psi_avg /= count; + double psi_avg_mag=psi_avg; + if (psi_avg==0.0) psi_avg_mag=1.0; + error = fabs(psi_avg-psi_avg_previous)/fabs(psi_avg_mag); + psi_avg_previous = psi_avg; + } } //************************************************************************/ stoptime = MPI_Wtime(); - if (rank==0) printf("-------------------------------------------------------------------\n"); + if (rank==0) printf("LB-Poission Solver: a steady-state solution is obtained\n"); + 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("******************* LB-Poisson Solver Statistics ********************\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"); + if (rank==0) printf("*********************************************************************\n"); } + +//void ScaLBL_Poisson::get_ElectricField(){ +//// ??? +//} diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 63e10df0..7eb7ac83 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -29,17 +29,19 @@ public: void Initialize(); void Run(double *ChargeDensity); - bool Restart,pBC; + //bool Restart,pBC; int timestep,timestepMax; + int analysis_interval; int BoundaryCondition; - double tau,mu; - double Fx,Fy,Fz,flux; - double din,dout; + double tau; double tolerance; + double k2_inv,deltaT; + double epsilon0,epsilon0_LB,epsilonR,epsilon_LB; int Nx,Ny,Nz,N,Np; int rank,nprocx,nprocy,nprocz,nprocs; double Lx,Ly,Lz; + double h;//image resolution std::shared_ptr Dm; // this domain is for analysis std::shared_ptr Mask; // this domain is for lbm @@ -51,9 +53,11 @@ public: IntArray Map; DoubleArray Distance; + DoubleArray Psi_host; int *NeighborList; double *fq; double *Psi; + double *ElectricField; private: MPI_Comm comm; diff --git a/tests/lbpm_electrokinetic_dfh_simulator.cpp b/tests/lbpm_electrokinetic_dfh_simulator.cpp index 5209565c..9725dda4 100644 --- a/tests/lbpm_electrokinetic_dfh_simulator.cpp +++ b/tests/lbpm_electrokinetic_dfh_simulator.cpp @@ -38,7 +38,7 @@ int main(int argc, char **argv) if (rank == 0){ printf("********************************************************\n"); - printf("Running Color LBM \n"); + printf("Running Electrokinetic LBM Simulator \n"); printf("********************************************************\n"); } PROFILE_ENABLE(1);