diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 664e23e1..4f611204 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -862,6 +862,145 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, bool UseSlippingVelBC, int times } } +void ScaLBL_Poisson::Run(double *ChargeDensity, DoubleArray MembraneDistance, bool UseSlippingVelBC, int timestep_from_Study){ + + double error = 1.0; + double threshold = 10000000.0; + bool SET_THRESHOLD = false; + if (electric_db->keyExists( "rescale_at_distance" )){ + SET_THRESHOLD = true; + threshold = electric_db->getScalar( "rescale_at_distance" ); + } + if (BoundaryConditionInlet > 0) SET_THRESHOLD = false; + if (BoundaryConditionOutlet > 0) SET_THRESHOLD = false; + + double *host_Error; + host_Error = new double [Np]; + + timestep=0; + auto t1 = std::chrono::system_clock::now(); + while (timestep < timestepMax && error > tolerance) { + //************************************************************************/ + // *************ODD TIMESTEP*************// + timestep++; + //SolveElectricPotentialAAodd(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential + SolvePoissonAAodd(ChargeDensity, UseSlippingVelBC);//perform collision + ScaLBL_Comm->Barrier(); comm.barrier(); + + // *************EVEN TIMESTEP*************// + timestep++; + //SolveElectricPotentialAAeven(timestep_from_Study,ChargeDensity, UseSlippingVelBC);//update electric potential + SolvePoissonAAeven(ChargeDensity, UseSlippingVelBC);//perform collision + ScaLBL_Comm->Barrier(); comm.barrier(); + //************************************************************************/ + + // Check convergence of steady-state solution + if (timestep==2){ + //save electric potential for convergence check + } + if (timestep%analysis_interval==0){ + /* get the elecric potential */ + ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz); + if (rank==0) printf(" ... getting Poisson solver error \n"); + double err = 0.0; + double max_error = 0.0; + ScaLBL_CopyToHost(host_Error,ResidualError,sizeof(double)*Np); + for (int idx=0; idx max_error ){ + max_error = err; + } + } + error=Dm->Comm.maxReduce(max_error); + + if (error > tolerance & SET_THRESHOLD){ + /* don't use this with an external BC */ + // cpompute the far-field electric potential + double inside_local = 0.0; + double outside_local = 0.0; + double inside_count_local = 0.0; + double outside_count_local = 0.0; + /* global values */ + double inside_global = 0.0; + double outside_global = 0.0; + double inside_count_global = 0.0; + double outside_count_global = 0.0; + for (int k=1; k threshold && distance < (threshold + 1.0)){ + outside_count_local += 1.0; + outside_local += Psi_host(n); + } + else if (distance < (-1.0)*threshold && distance > (-1.0)*(threshold + 1.0)){ + inside_count_local += 1.0; + inside_local += Psi_host(n); + } + } + } + } + inside_count_global = Dm->Comm.sumReduce(inside_count_local); + outside_count_global = Dm->Comm.sumReduce(outside_count_local); + outside_global = Dm->Comm.sumReduce(outside_local); + inside_global = Dm->Comm.sumReduce(inside_local); + outside_global /= outside_count_global; + inside_global /= inside_count_global; + + if (rank==0) printf(" Rescaling far-field electric potential to value (outside): %f \n",outside_global); + if (rank==0) printf(" Rescaling far-field electric potential to value (inside): %f \n",inside_global); + + // rescale the far-field electric potential + for (int k=1; k (threshold + 1.0)){ + Psi_host(n) = outside_global; + } + else if ( distance < (-1.0)*(threshold + 1.0)){ + Psi_host(n) = inside_global; + } + } + } + } + ScaLBL_CopyToDevice(Psi,Psi_host.data(),sizeof(double)*Nx*Ny*Nz); + } + /* compute the eletric field */ + //ScaLBL_D3Q19_Poisson_getElectricField(fq, ElectricField, tau, Np); + } + } + if (rank == 0) + printf("---------------------------------------------------------------" + "----\n"); + // Compute the walltime per timestep + auto t2 = std::chrono::system_clock::now(); + double cputime = std::chrono::duration(t2 - t1).count() / timestep; + // Performance obtained from each node + double MLUPS = double(Np) / cputime / 1000000; + + if (rank == 0) + printf("********************************************************\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"); + delete [] host_Error; + + //************************************************************************/ + + if(WriteLog==true){ + getConvergenceLog(timestep,error); + } +} + void ScaLBL_Poisson::getConvergenceLog(int timestep,double error){ if ( rank == 0 ) { fprintf(TIMELOG,"%i %.5g\n",timestep,error); diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 8ab60140..faee8c01 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -35,6 +35,8 @@ public: void Create(); void Initialize(double time_conv_from_Study); void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study); + void Run(double *ChargeDensity, DoubleArray MembraneDistance, + 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/tests/lbpm_nernst_planck_cell_simulator.cpp b/tests/lbpm_nernst_planck_cell_simulator.cpp index 3ac2065c..f98ed1b7 100644 --- a/tests/lbpm_nernst_planck_cell_simulator.cpp +++ b/tests/lbpm_nernst_planck_cell_simulator.cpp @@ -125,7 +125,7 @@ int main(int argc, char **argv) while (timestep < Study.timestepMax){ timestep++; - PoissonSolver.Run(IonModel.ChargeDensity,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental + PoissonSolver.Run(IonModel.ChargeDensity,IonModel.MembraneDistance,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental //comm.barrier(); //if (rank == 0) printf(" Poisson step %i \n",timestep); //StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity