From 9f76e7b1e8251a9ab87ca36b4e1ad8b93afd1b4d Mon Sep 17 00:00:00 2001 From: James McClure Date: Fri, 6 May 2022 19:03:19 -0400 Subject: [PATCH] update to cell vis routine --- analysis/ElectroChemistry.cpp | 407 +++++++++++++++++++- analysis/ElectroChemistry.h | 5 + tests/lbpm_cell_simulator.cpp | 1 - tests/lbpm_nernst_planck_cell_simulator.cpp | 22 +- 4 files changed, 421 insertions(+), 14 deletions(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index fcd4d2a2..99880455 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -49,7 +49,7 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr dm) IonFluxElectrical_y.fill(0); IonFluxElectrical_z.resize(Nx, Ny, Nz); IonFluxElectrical_z.fill(0); - + if (Dm->rank() == 0) { bool WriteHeader = false; TIMELOG = fopen("electrokinetic.csv", "r"); @@ -595,3 +595,408 @@ void ElectroChemistryAnalyzer::WriteVis(ScaLBL_IonModel &Ion, } */ } + +void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, + ScaLBL_Poisson &Poisson, + int timestep) { + + int i, j, k; + double Vin = 0.0; + double Vout = 0.0; + Poisson.getElectricPotential(ElectricalPotential); + + /* local sub-domain averages */ + double *rho_avg_local; + double *rho_mu_avg_local; + double *rho_mu_fluctuation_local; + double *rho_psi_avg_local; + double *rho_psi_fluctuation_local; + /* global averages */ + double *rho_avg_global; + double *rho_mu_avg_global; + double *rho_mu_fluctuation_global; + double *rho_psi_avg_global; + double *rho_psi_fluctuation_global; + + /* Get the distance to the membrane */ + if (Ion.USE_MEMBRANE){ + //Ion.MembraneDistance; + } + + /* local sub-domain averages */ + rho_avg_local = new double[Ion.number_ion_species]; + rho_mu_avg_local = new double[Ion.number_ion_species]; + rho_mu_fluctuation_local = new double[Ion.number_ion_species]; + rho_psi_avg_local = new double[Ion.number_ion_species]; + rho_psi_fluctuation_local = new double[Ion.number_ion_species]; + /* global averages */ + rho_avg_global = new double[Ion.number_ion_species]; + rho_mu_avg_global = new double[Ion.number_ion_species]; + rho_mu_fluctuation_global = new double[Ion.number_ion_species]; + rho_psi_avg_global = new double[Ion.number_ion_species]; + rho_psi_fluctuation_global = new double[Ion.number_ion_species]; + + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + rho_avg_local[ion] = 0.0; + rho_mu_avg_local[ion] = 0.0; + rho_psi_avg_local[ion] = 0.0; + Ion.getIonConcentration(Rho, ion); + /* Compute averages for each ion */ + for (k = 1; k < Nz; k++) { + for (j = 1; j < Ny; j++) { + for (i = 1; i < Nx; i++) { + rho_avg_local[ion] += Rho(i, j, k); + rho_mu_avg_local[ion] += Rho(i, j, k) * Rho(i, j, k); + rho_psi_avg_local[ion] += + Rho(i, j, k) * ElectricalPotential(i, j, k); + } + } + } + rho_avg_global[ion] = Dm->Comm.sumReduce(rho_avg_local[ion]) / Volume; + rho_mu_avg_global[ion] = + Dm->Comm.sumReduce(rho_mu_avg_local[ion]) / Volume; + rho_psi_avg_global[ion] = + Dm->Comm.sumReduce(rho_psi_avg_local[ion]) / Volume; + + if (rho_avg_global[ion] > 0.0) { + rho_mu_avg_global[ion] /= rho_avg_global[ion]; + rho_psi_avg_global[ion] /= rho_avg_global[ion]; + } + } + + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + rho_mu_fluctuation_local[ion] = 0.0; + rho_psi_fluctuation_local[ion] = 0.0; + /* Compute averages for each ion */ + for (k = 1; k < Nz; k++) { + for (j = 1; j < Ny; j++) { + for (i = 1; i < Nx; i++) { + rho_mu_fluctuation_local[ion] += + (Rho(i, j, k) * Rho(i, j, k) - rho_mu_avg_global[ion]); + rho_psi_fluctuation_local[ion] += + (Rho(i, j, k) * ElectricalPotential(i, j, k) - + rho_psi_avg_global[ion]); + } + } + } + rho_mu_fluctuation_global[ion] = + Dm->Comm.sumReduce(rho_mu_fluctuation_local[ion]); + rho_psi_fluctuation_global[ion] = + Dm->Comm.sumReduce(rho_psi_fluctuation_local[ion]); + } + + if (Dm->rank() == 0) { + fprintf(TIMELOG, "%i ", timestep); + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + fprintf(TIMELOG, "%.8g ", rho_avg_global[ion]); + fprintf(TIMELOG, "%.8g ", rho_mu_avg_global[ion]); + fprintf(TIMELOG, "%.8g ", rho_psi_avg_global[ion]); + fprintf(TIMELOG, "%.8g ", rho_mu_fluctuation_global[ion]); + fprintf(TIMELOG, "%.8g ", rho_psi_fluctuation_global[ion]); + } + fprintf(TIMELOG, "%.8g %.8g\n", Vin, Vout); + fflush(TIMELOG); + } + /* else{ + fprintf(TIMELOG,"%i ",timestep); + for (int ion=0; ion input_db, + int timestep) { + + auto vis_db = input_db->getDatabase("Visualization"); + char VisName[40]; + auto format = vis_db->getWithDefault( "format", "hdf5" ); + + std::vector visData; + fillHalo fillData(Dm->Comm, Dm->rank_info, + {Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2}, {1, 1, 1}, + 0, 1); + + IO::initialize("",format,"false"); + // Create the MeshDataStruct + visData.resize(1); + + visData[0].meshName = "domain"; + visData[0].mesh = + std::make_shared(Dm->rank_info, Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2, Dm->Lx, Dm->Ly, Dm->Lz); + //electric potential + auto ElectricPotentialVar = std::make_shared(); + //electric field + auto ElectricFieldVar_x = std::make_shared(); + auto ElectricFieldVar_y = std::make_shared(); + auto ElectricFieldVar_z = std::make_shared(); + + //ion concentration + std::vector> IonConcentration; + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + IonConcentration.push_back(std::make_shared()); + } + + // diffusive ion flux + std::vector> IonFluxDiffusive; + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + //push in x-,y-, and z-component for each ion species + IonFluxDiffusive.push_back(std::make_shared()); + IonFluxDiffusive.push_back(std::make_shared()); + IonFluxDiffusive.push_back(std::make_shared()); + } + + // electro-migrational ion flux + std::vector> IonFluxElectrical; + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + //push in x-,y-, and z-component for each ion species + IonFluxElectrical.push_back(std::make_shared()); + IonFluxElectrical.push_back(std::make_shared()); + IonFluxElectrical.push_back(std::make_shared()); + } + //-------------------------------------------------------------------------------------------------------------------- + + //-------------------------------------Create Names for Variables------------------------------------------------------ + if (vis_db->getWithDefault("save_electric_potential", true)) { + ElectricPotentialVar->name = "ElectricPotential"; + ElectricPotentialVar->type = IO::VariableType::VolumeVariable; + ElectricPotentialVar->dim = 1; + ElectricPotentialVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2); + visData[0].vars.push_back(ElectricPotentialVar); + } + + if (vis_db->getWithDefault("save_concentration", true)) { + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + sprintf(VisName, "IonConcentration_%zu", ion + 1); + IonConcentration[ion]->name = VisName; + IonConcentration[ion]->type = IO::VariableType::VolumeVariable; + IonConcentration[ion]->dim = 1; + IonConcentration[ion]->data.resize(Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2); + visData[0].vars.push_back(IonConcentration[ion]); + } + } + + if (vis_db->getWithDefault("save_ion_flux_diffusive", false)) { + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + // x-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxDiffusive_x", ion + 1); + IonFluxDiffusive[3 * ion + 0]->name = VisName; + IonFluxDiffusive[3 * ion + 0]->type = + IO::VariableType::VolumeVariable; + IonFluxDiffusive[3 * ion + 0]->dim = 1; + IonFluxDiffusive[3 * ion + 0]->data.resize(Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2); + visData[0].vars.push_back(IonFluxDiffusive[3 * ion + 0]); + // y-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxDiffusive_y", ion + 1); + IonFluxDiffusive[3 * ion + 1]->name = VisName; + IonFluxDiffusive[3 * ion + 1]->type = + IO::VariableType::VolumeVariable; + IonFluxDiffusive[3 * ion + 1]->dim = 1; + IonFluxDiffusive[3 * ion + 1]->data.resize(Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2); + visData[0].vars.push_back(IonFluxDiffusive[3 * ion + 1]); + // z-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxDiffusive_z", ion + 1); + IonFluxDiffusive[3 * ion + 2]->name = VisName; + IonFluxDiffusive[3 * ion + 2]->type = + IO::VariableType::VolumeVariable; + IonFluxDiffusive[3 * ion + 2]->dim = 1; + IonFluxDiffusive[3 * ion + 2]->data.resize(Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2); + visData[0].vars.push_back(IonFluxDiffusive[3 * ion + 2]); + } + } + + + if (vis_db->getWithDefault("save_ion_flux_electrical", false)) { + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + // x-component of electro-migrational flux + sprintf(VisName, "Ion%zu_FluxElectrical_x", ion + 1); + IonFluxElectrical[3 * ion + 0]->name = VisName; + IonFluxElectrical[3 * ion + 0]->type = + IO::VariableType::VolumeVariable; + IonFluxElectrical[3 * ion + 0]->dim = 1; + IonFluxElectrical[3 * ion + 0]->data.resize(Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2); + visData[0].vars.push_back(IonFluxElectrical[3 * ion + 0]); + // y-component of electro-migrational flux + sprintf(VisName, "Ion%zu_FluxElectrical_y", ion + 1); + IonFluxElectrical[3 * ion + 1]->name = VisName; + IonFluxElectrical[3 * ion + 1]->type = + IO::VariableType::VolumeVariable; + IonFluxElectrical[3 * ion + 1]->dim = 1; + IonFluxElectrical[3 * ion + 1]->data.resize(Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2); + visData[0].vars.push_back(IonFluxElectrical[3 * ion + 1]); + // z-component of electro-migrational flux + sprintf(VisName, "Ion%zu_FluxElectrical_z", ion + 1); + IonFluxElectrical[3 * ion + 2]->name = VisName; + IonFluxElectrical[3 * ion + 2]->type = + IO::VariableType::VolumeVariable; + IonFluxElectrical[3 * ion + 2]->dim = 1; + IonFluxElectrical[3 * ion + 2]->data.resize(Dm->Nx - 2, Dm->Ny - 2, + Dm->Nz - 2); + visData[0].vars.push_back(IonFluxElectrical[3 * ion + 2]); + } + } + + if (vis_db->getWithDefault("save_electric_field", false)) { + ElectricFieldVar_x->name = "ElectricField_x"; + ElectricFieldVar_x->type = IO::VariableType::VolumeVariable; + ElectricFieldVar_x->dim = 1; + ElectricFieldVar_x->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2); + visData[0].vars.push_back(ElectricFieldVar_x); + ElectricFieldVar_y->name = "ElectricField_y"; + ElectricFieldVar_y->type = IO::VariableType::VolumeVariable; + ElectricFieldVar_y->dim = 1; + ElectricFieldVar_y->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2); + visData[0].vars.push_back(ElectricFieldVar_y); + ElectricFieldVar_z->name = "ElectricField_z"; + ElectricFieldVar_z->type = IO::VariableType::VolumeVariable; + ElectricFieldVar_z->dim = 1; + ElectricFieldVar_z->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2); + visData[0].vars.push_back(ElectricFieldVar_z); + } + //-------------------------------------------------------------------------------------------------------------------- + + //------------------------------------Save All Variables-------------------------------------------------------------- + if (vis_db->getWithDefault("save_electric_potential", true)) { + ASSERT(visData[0].vars[0]->name == "ElectricPotential"); + Poisson.getElectricPotential(ElectricalPotential); + Array &ElectricPotentialData = visData[0].vars[0]->data; + fillData.copy(ElectricalPotential, ElectricPotentialData); + } + + if (vis_db->getWithDefault("save_concentration", true)) { + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + sprintf(VisName, "IonConcentration_%zu", ion + 1); + //IonConcentration[ion]->name = VisName; + ASSERT(visData[0].vars[1 + ion]->name == VisName); + Array &IonConcentrationData = + visData[0].vars[1 + ion]->data; + Ion.getIonConcentration(Rho, ion); + fillData.copy(Rho, IonConcentrationData); + } + } + + if (vis_db->getWithDefault("save_ion_flux_diffusive", false)) { + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + + // x-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxDiffusive_x", ion + 1); + //IonFluxDiffusive[3*ion+0]->name = VisName; + ASSERT(visData[0] + .vars[4 + Ion.number_ion_species + 3 * ion + 0] + ->name == VisName); + // y-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxDiffusive_y", ion + 1); + //IonFluxDiffusive[3*ion+1]->name = VisName; + ASSERT(visData[0] + .vars[4 + Ion.number_ion_species + 3 * ion + 1] + ->name == VisName); + // z-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxDiffusive_z", ion + 1); + //IonFluxDiffusive[3*ion+2]->name = VisName; + ASSERT(visData[0] + .vars[4 + Ion.number_ion_species + 3 * ion + 2] + ->name == VisName); + + Array &IonFluxData_x = + visData[0].vars[4 + Ion.number_ion_species + 3 * ion + 0]->data; + Array &IonFluxData_y = + visData[0].vars[4 + Ion.number_ion_species + 3 * ion + 1]->data; + Array &IonFluxData_z = + visData[0].vars[4 + Ion.number_ion_species + 3 * ion + 2]->data; + Ion.getIonFluxDiffusive(IonFluxDiffusive_x, IonFluxDiffusive_y, + IonFluxDiffusive_z, ion); + fillData.copy(IonFluxDiffusive_x, IonFluxData_x); + fillData.copy(IonFluxDiffusive_y, IonFluxData_y); + fillData.copy(IonFluxDiffusive_z, IonFluxData_z); + } + } + + if (vis_db->getWithDefault("save_ion_flux_electrical", false)) { + for (size_t ion = 0; ion < Ion.number_ion_species; ion++) { + + // x-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxElectrical_x", ion + 1); + //IonFluxDiffusive[3*ion+0]->name = VisName; + ASSERT(visData[0] + .vars[4 + Ion.number_ion_species * (1 + 6) + 3 * ion + 0] + ->name == VisName); + // y-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxElectrical_y", ion + 1); + //IonFluxDiffusive[3*ion+1]->name = VisName; + ASSERT(visData[0] + .vars[4 + Ion.number_ion_species * (1 + 6) + 3 * ion + 1] + ->name == VisName); + // z-component of diffusive flux + sprintf(VisName, "Ion%zu_FluxElectrical_z", ion + 1); + //IonFluxDiffusive[3*ion+2]->name = VisName; + ASSERT(visData[0] + .vars[4 + Ion.number_ion_species * (1 + 6) + 3 * ion + 2] + ->name == VisName); + + Array &IonFluxData_x = + visData[0] + .vars[4 + Ion.number_ion_species * (1 + 6) + 3 * ion + 0] + ->data; + Array &IonFluxData_y = + visData[0] + .vars[4 + Ion.number_ion_species * (1 + 6) + 3 * ion + 1] + ->data; + Array &IonFluxData_z = + visData[0] + .vars[4 + Ion.number_ion_species * (1 + 6) + 3 * ion + 2] + ->data; + Ion.getIonFluxElectrical(IonFluxElectrical_x, IonFluxElectrical_y, + IonFluxElectrical_z, ion); + fillData.copy(IonFluxElectrical_x, IonFluxData_x); + fillData.copy(IonFluxElectrical_y, IonFluxData_y); + fillData.copy(IonFluxElectrical_z, IonFluxData_z); + } + } + + if (vis_db->getWithDefault("save_electric_field", false)) { + ASSERT( + visData[0].vars[4 + Ion.number_ion_species * (1 + 9) + 0]->name == + "ElectricField_x"); + ASSERT( + visData[0].vars[4 + Ion.number_ion_species * (1 + 9) + 1]->name == + "ElectricField_y"); + ASSERT( + visData[0].vars[4 + Ion.number_ion_species * (1 + 9) + 2]->name == + "ElectricField_z"); + Poisson.getElectricField(ElectricalField_x, ElectricalField_y, + ElectricalField_z); + Array &ElectricalFieldxData = + visData[0].vars[4 + Ion.number_ion_species * (1 + 9) + 0]->data; + Array &ElectricalFieldyData = + visData[0].vars[4 + Ion.number_ion_species * (1 + 9) + 1]->data; + Array &ElectricalFieldzData = + visData[0].vars[4 + Ion.number_ion_species * (1 + 9) + 2]->data; + fillData.copy(ElectricalField_x, ElectricalFieldxData); + fillData.copy(ElectricalField_y, ElectricalFieldyData); + fillData.copy(ElectricalField_z, ElectricalFieldzData); + } + + if (vis_db->getWithDefault("write_silo", true)) + IO::writeData(timestep, visData, Dm->Comm); + //-------------------------------------------------------------------------------------------------------------------- + /* if (vis_db->getWithDefault( "save_8bit_raw", true )){ + char CurrentIDFilename[40]; + sprintf(CurrentIDFilename,"id_t%d.raw",timestep); + Averages.AggregateLabels(CurrentIDFilename); + } +*/ +} \ No newline at end of file diff --git a/analysis/ElectroChemistry.h b/analysis/ElectroChemistry.h index 1b739611..34a12501 100644 --- a/analysis/ElectroChemistry.h +++ b/analysis/ElectroChemistry.h @@ -29,6 +29,8 @@ public: double nu_n, nu_w; double gamma_wn, beta; double Fx, Fy, Fz; + + bool USE_MEMBRANE; //........................................................................... int Nx, Ny, Nz; @@ -62,6 +64,9 @@ public: void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes, std::shared_ptr input_db, int timestep); + void Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, int timestep); + void WriteVis(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, + std::shared_ptr input_db, int timestep); private: FILE *TIMELOG; diff --git a/tests/lbpm_cell_simulator.cpp b/tests/lbpm_cell_simulator.cpp index 7e3d255b..dd0b0d52 100644 --- a/tests/lbpm_cell_simulator.cpp +++ b/tests/lbpm_cell_simulator.cpp @@ -138,7 +138,6 @@ int main(int argc, char **argv) PoissonSolver.getElectricPotential_debug(timestep); PoissonSolver.getElectricField_debug(timestep); IonModel.getIonConcentration_debug(timestep); - } } diff --git a/tests/lbpm_nernst_planck_cell_simulator.cpp b/tests/lbpm_nernst_planck_cell_simulator.cpp index 488da489..7e80e357 100644 --- a/tests/lbpm_nernst_planck_cell_simulator.cpp +++ b/tests/lbpm_nernst_planck_cell_simulator.cpp @@ -81,7 +81,7 @@ int main(int argc, char **argv) fflush(stdout); // Create analysis object - //ElectroChemistryAnalyzer Analysis(IonModel.Dm); + ElectroChemistryAnalyzer Analysis(IonModel.Dm); // Get internal iteration number //StokesModel.timestepMax = Study.getStokesNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv); @@ -137,24 +137,22 @@ int main(int argc, char **argv) fflush(stdout); - //if (timestep%Study.analysis_interval==0){ - // Analysis.Basic(IonModel,PoissonSolver,StokesModel,timestep); - //} + if (timestep%Study.analysis_interval==0){ + Analysis.Basic(IonModel,PoissonSolver,timestep); + } if (timestep%Study.visualization_interval==0){ - //Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep); - // PoissonSolver.getElectricPotential(timestep); - //PoissonSolver.getElectricField(timestep); - //IonModel.getIonConcentration(timestep); + + Analysis.WriteVis(IonModel,PoissonSolver,Study.db,timestep); //StokesModel.getVelocity(timestep); - PoissonSolver.getElectricPotential_debug(timestep); - PoissonSolver.getElectricField_debug(timestep); - IonModel.getIonConcentration_debug(timestep); + //PoissonSolver.getElectricPotential_debug(timestep); + // PoissonSolver.getElectricField_debug(timestep); + //IonModel.getIonConcentration_debug(timestep); } } if (rank==0) printf("Save simulation raw data at maximum timestep\n"); - //Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep); + Analysis.WriteVis(IonModel,PoissonSolver,Study.db,timestep); if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); if (rank==0) printf("*************************************************************\n");