diff --git a/models/IonModel.cpp b/models/IonModel.cpp index d1bcfe5a..f976b1f0 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -61,6 +61,9 @@ void ScaLBL_IonModel::ReadParams(string filename, vector &num_iter) { Ex_dummy = 0.0; //for debugging, unit [V/m] Ey_dummy = 0.0; //for debugging, unit [V/m] Ez_dummy = 0.0; //for debugging, unit [V/m] + + sprintf(LocalRankString, "%05d", rank); + sprintf(LocalRestartFile, "%s%s", "Restart.", LocalRankString); //--------------------------------------------------------------------------// BoundaryConditionSolid = 0; @@ -329,6 +332,9 @@ void ScaLBL_IonModel::ReadParams(string filename) { Ex_dummy = 0.0; //for debugging, unit [V/m] Ey_dummy = 0.0; //for debugging, unit [V/m] Ez_dummy = 0.0; //for debugging, unit [V/m] + + sprintf(LocalRankString, "%05d", rank); + sprintf(LocalRestartFile, "%s%s", "Restart.", LocalRankString); //--------------------------------------------------------------------------// // Read domain parameters @@ -345,6 +351,9 @@ void ScaLBL_IonModel::ReadParams(string filename) { if (ion_db->keyExists("tolerance")) { tolerance = ion_db->getScalar("tolerance"); } + if (ion_db->keyExists("Restart")) { + Restart = ion_db->getScalar("Restart"); + } if (ion_db->keyExists("temperature")) { T = ion_db->getScalar("temperature"); //re-calculate thermal voltage @@ -1114,6 +1123,43 @@ void ScaLBL_IonModel::Initialize() { IonConcentration[ic], Np); } } + /** RESTART **/ + if (Restart == true) { + if (rank == 0) { + printf(" ION MODEL: Reading restart file! \n"); + } + double*cDist; + double *Ci_host; + cDist = new double[7 * number_ion_species * Np]; + Ci_host = new double[number_ion_species * Np]; + + ifstream File(LocalRestartFile, ios::binary); + int idx; + double value,sum; + // Read the distributions + for (size_t ic = 0; ic < number_ion_species; ic++){ + for (int n = 0; n < Np; n++) { + sum = 0.0; + for (int q = 0; q < 7; q++) { + File.read((char *)&value, sizeof(value)); + cDist[ic * q * Np + n] = value; + sum += value; + } + Ci_host[ic * Np + n] = sum; + } + } + File.close(); + + // Copy the restart data to the GPU + ScaLBL_CopyToDevice(Ci, Ci_host, Np * number_ion_species* sizeof(double)); + ScaLBL_CopyToDevice(fq, cDist, 7 * Np * number_ion_species *sizeof(double)); + ScaLBL_Comm->Barrier(); + comm.barrier(); + delete[] Ci_host; + delete[] cDist; + } + /** END RESTART **/ + if (rank == 0) printf("LB Ion Solver: initializing charge density\n"); for (size_t ic = 0; ic < number_ion_species; ic++) { @@ -1594,6 +1640,32 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl //if (rank==0) printf("********************************************************\n"); } +void ScaLBL_IonModel::Checkpoint(){ + + if (rank == 0) { + printf(" ION MODEL: Writing restart file! \n"); + } + int idx; + double value,sum; + double*cDist; + cDist = new double[7 * number_ion_species * Np]; + ScaLBL_CopyToHost(cDist, fq, 7 * Np * number_ion_species *sizeof(double)); + + ofstream File(LocalRestartFile, ios::binary); + for (size_t ic = 0; ic < number_ion_species; ic++){ + for (int n = 0; n < Np; n++) { + // Write the distributions + for (int q = 0; q < 7; q++) { + value = cDist[ic * q * Np + n]; + File.write((char *)&value, sizeof(value)); + } + } + } + File.close(); + + delete[] cDist; + +} void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const size_t ic) { diff --git a/models/IonModel.h b/models/IonModel.h index f3f844ec..9a6756e6 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -49,9 +49,10 @@ public: void getIonFluxElectrical_debug(int timestep); void DummyFluidVelocity(); void DummyElectricField(); + void Checkpoint(); double CalIonDenConvergence(vector &ci_avg_previous); - //bool Restart,pBC; + bool Restart; int timestep; vector timestepMax; int BoundaryConditionSolid; diff --git a/models/MultiPhysController.cpp b/models/MultiPhysController.cpp index dc361717..b0d3c817 100644 --- a/models/MultiPhysController.cpp +++ b/models/MultiPhysController.cpp @@ -16,6 +16,7 @@ void ScaLBL_Multiphys_Controller::ReadParams(string filename) { // Default parameters timestepMax = 10000; Restart = false; + restart_interval = 100000; num_iter_Stokes = 1; num_iter_Ion.push_back(1); analysis_interval = 500; @@ -35,6 +36,9 @@ void ScaLBL_Multiphys_Controller::ReadParams(string filename) { visualization_interval = study_db->getScalar("visualization_interval"); } + if (study_db->keyExists("restart_interval")) { + restart_interval = study_db->getScalar("restart_interval"); + } if (study_db->keyExists("tolerance")) { tolerance = study_db->getScalar("tolerance"); } diff --git a/models/MultiPhysController.h b/models/MultiPhysController.h index d41a6db8..baee9966 100644 --- a/models/MultiPhysController.h +++ b/models/MultiPhysController.h @@ -38,6 +38,7 @@ public: vector num_iter_Ion; int analysis_interval; int visualization_interval; + int restart_interval; double tolerance; double time_conv_max; double time_conv_MainLoop; diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index 3c41a151..e7fb631d 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -66,8 +66,12 @@ void ScaLBL_Poisson::ReadParams(string filename){ TestPeriodicTime = 1.0;//unit: [sec] TestPeriodicTimeConv = 0.01; //unit [sec/lt] TestPeriodicSaveInterval = 0.1; //unit [sec] + Restart = "false"; // LB-Poisson Model parameters + if (electric_db->keyExists( "Restart" )){ + Restart = electric_db->getScalar("Restart"); + } if (electric_db->keyExists( "timestepMax" )){ timestepMax = electric_db->getScalar( "timestepMax" ); } @@ -134,6 +138,10 @@ void ScaLBL_Poisson::ReadParams(string filename){ //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 + + /* restart string */ + sprintf(LocalRankString, "%05d", rank); + sprintf(LocalRestartFile, "%s%s", "Psi.", LocalRankString); if (rank==0) printf("***********************************************************************************\n"); if (rank==0) printf("LB-Poisson Solver: steady-state MaxTimeStep = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance); @@ -203,7 +211,7 @@ void ScaLBL_Poisson::ReadInput(){ sprintf(LocalRankString,"%05d",Dm->rank()); sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); - sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString); + sprintf(LocalRestartFile,"%s%s","Psi.",LocalRankString); if (domain_db->keyExists( "Filename" )){ @@ -566,6 +574,21 @@ void ScaLBL_Poisson::Potential_Init(double *psi_init){ else{//mixed periodic and non-periodic BCs are not supported! ERROR("Error: check the type of inlet and outlet boundary condition! Mixed periodic and non-periodic BCs are found!\n"); } + /** RESTART **/ + if (Restart == true) { + if (rank == 0) { + printf(" POISSON MODEL: Reading restart file! \n"); + } + ifstream File(LocalRestartFile, ios::binary); + double value; + // Read the distributions + for (int n = 0; n < Nx*Ny*Nz; n++) { + File.read((char *)&value, sizeof(value)); + psi_init[ n] = value; + } + File.close(); + } + /** END RESTART **/ } double ScaLBL_Poisson::getBoundaryVoltagefromPeriodicBC(double V0, double freq, double phase_shift, int time_step){ @@ -1040,6 +1063,27 @@ void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingV } } +void ScaLBL_Poisson::Checkpoint(){ + + if (rank == 0) { + printf(" POISSON MODEL: Writing restart file! \n"); + } + double value; + double *cPsi; + cPsi = new double[Nx*Ny*Nz]; + ScaLBL_CopyToHost(cPsi, Psi, Nx*Ny*Nz *sizeof(double)); + + ofstream File(LocalRestartFile, ios::binary); + for (int n = 0; n < Nx*Ny*Nz; n++) { + value = cPsi[n]; + File.write((char *)&value, sizeof(value)); + } + + File.close(); + + delete[] cPsi; +} + void ScaLBL_Poisson::DummyChargeDensity(){ double *ChargeDensity_host; ChargeDensity_host = new double[Np]; diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 1e162c9c..8ab60140 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -40,9 +40,11 @@ public: void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y, DoubleArray &Values_z); void getElectricField_debug(int timestep); + void Checkpoint(); + void DummyChargeDensity(); //for debugging - //bool Restart,pBC; + bool Restart; int timestep, timestepMax; int analysis_interval; int BoundaryConditionInlet; diff --git a/tests/lbpm_nernst_planck_cell_simulator.cpp b/tests/lbpm_nernst_planck_cell_simulator.cpp index 2cb206c5..3ac2065c 100644 --- a/tests/lbpm_nernst_planck_cell_simulator.cpp +++ b/tests/lbpm_nernst_planck_cell_simulator.cpp @@ -149,6 +149,11 @@ int main(int argc, char **argv) //IonModel.getIonConcentration_debug(timestep); } + + if (timestep%Study.restart_interval==0){ + IonModel.Checkpoint(); + PoissonSolver.Checkpoint(); + } } if (rank==0) printf("Save simulation raw data at maximum timestep\n");