From 9826ef5624f65f05ae853399d129d654786506d7 Mon Sep 17 00:00:00 2001 From: James McClure Date: Tue, 29 Dec 2020 14:04:43 -0500 Subject: [PATCH] adding silo vis capabilities to electrochem --- analysis/ElectroChemistry.cpp | 119 ++++++++++++++++-- analysis/ElectroChemistry.h | 13 +- models/IonModel.cpp | 16 +-- models/IonModel.h | 7 +- models/PoissonSolver.cpp | 32 ++--- models/PoissonSolver.h | 8 +- models/StokesModel.cpp | 17 ++- models/StokesModel.h | 6 +- tests/TestIonModel.cpp | 2 +- tests/TestNernstPlanck.cpp | 6 +- tests/TestPNP_Stokes.cpp | 8 +- tests/TestPoissonSolver.cpp | 4 +- ...m_electrokinetic_SingleFluid_simulator.cpp | 17 +-- 13 files changed, 176 insertions(+), 79 deletions(-) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index 84d618b1..b052d459 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -1,7 +1,11 @@ #include "analysis/ElectroChemistry.h" ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr dm): - Dm(dm){ + Dm(dm), + fillData(dm->Comm,dm->rank_info,{dm->Nx-2,dm->Ny-2,dm->Nz-2},{1,1,1},0,1) +{ + + MPI_Comm_dup(dm->Comm,&comm); Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz; Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0; @@ -14,15 +18,6 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr dm): Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0); SDs.resize(Nx,Ny,Nz); SDs.fill(0); - DoubleArray Rho; // density field - DoubleArray ChemicalPotential; // density field - DoubleArray ElectricalPotential; // density field - DoubleArray Pressure; // pressure field - DoubleArray Vel_x; // velocity field - DoubleArray Vel_y; - DoubleArray Vel_z; - DoubleArray SDs; - if (Dm->rank()==0){ bool WriteHeader=false; TIMELOG = fopen("electrokinetic.csv","r"); @@ -36,7 +31,7 @@ ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr dm): { // If timelog is empty, write a short header to list the averages //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"sw krw krn vw vn pw pn\n"); + fprintf(TIMELOG,"TBD TBD\n"); } } @@ -50,10 +45,108 @@ void ElectroChemistryAnalyzer::SetParams(){ } -void ElectroChemistryAnalyzer::Basic(){ +void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes){ + + Poisson.getElectricPotential(ElectricalPotential); + for (int ion=0; ion input_db, int timestep){ + auto vis_db = input_db->getDatabase( "Visualization" ); + char VisName[40]; + + IO::initialize("","silo","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 ); + auto ElectricPotential = std::make_shared(); + auto IonConcentration = std::make_shared(); + auto VxVar = std::make_shared(); + auto VyVar = std::make_shared(); + auto VzVar = std::make_shared(); + + if (vis_db->getWithDefault( "save_electric_potential", true )){ + ElectricPotential->name = "ElectricPotential"; + ElectricPotential->type = IO::VariableType::VolumeVariable; + ElectricPotential->dim = 1; + ElectricPotential->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(ElectricPotential); + } + + if (vis_db->getWithDefault( "save_concentration", true )){ + for (int ion=0; ionname = VisName; + IonConcentration->type = IO::VariableType::VolumeVariable; + IonConcentration->dim = 1; + IonConcentration->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(IonConcentration); + } + + } + if (vis_db->getWithDefault( "save_velocity", false )){ + VxVar->name = "Velocity_x"; + VxVar->type = IO::VariableType::VolumeVariable; + VxVar->dim = 1; + VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VxVar); + VyVar->name = "Velocity_y"; + VyVar->type = IO::VariableType::VolumeVariable; + VyVar->dim = 1; + VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VyVar); + VzVar->name = "Velocity_z"; + VzVar->type = IO::VariableType::VolumeVariable; + VzVar->dim = 1; + VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + visData[0].vars.push_back(VzVar); + } + + 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 (int ion=0; ionname = VisName; + ASSERT(visData[0].vars[1]->name==VisName); + Array& IonConcentrationData = visData[0].vars[1]->data; + Ion.getIonConcentration(Rho,ion); + fillData.copy(Rho,IonConcentrationData); + } + } + + if (vis_db->getWithDefault( "save_velocity", false )){ + ASSERT(visData[0].vars[2]->name=="Velocity_x"); + ASSERT(visData[0].vars[3]->name=="Velocity_y"); + ASSERT(visData[0].vars[4]->name=="Velocity_z"); + Stokes.getVelocity(Vel_x,Vel_y,Vel_z); + Array& VelxData = visData[0].vars[2]->data; + Array& VelyData = visData[0].vars[3]->data; + Array& VelzData = visData[0].vars[4]->data; + fillData.copy(Vel_x,VelxData); + fillData.copy(Vel_y,VelyData); + fillData.copy(Vel_z,VelzData); + } + + if (vis_db->getWithDefault( "write_silo", true )) + IO::writeData( timestep, visData, comm ); + +/* if (vis_db->getWithDefault( "save_8bit_raw", true )){ + char CurrentIDFilename[40]; + sprintf(CurrentIDFilename,"id_t%d.raw",timestep); + Averages.AggregateLabels(CurrentIDFilename); + } +*/ } diff --git a/analysis/ElectroChemistry.h b/analysis/ElectroChemistry.h index ad3a3578..fa404fb0 100644 --- a/analysis/ElectroChemistry.h +++ b/analysis/ElectroChemistry.h @@ -1,5 +1,5 @@ /* - * Sub-phase averaging tools + * averaging tools for electrochemistry */ #ifndef ElectroChem_INC @@ -16,9 +16,14 @@ #include "IO/MeshDatabase.h" #include "IO/Reader.h" #include "IO/Writer.h" +#include "models/IonModel.h" +#include "models/PoissonSolver.h" +#include "models/StokesModel.h" class ElectroChemistryAnalyzer{ public: + MPI_Comm comm; + int tag; std::shared_ptr Dm; double Volume; // input variables @@ -42,10 +47,12 @@ public: ~ElectroChemistryAnalyzer(); void SetParams(); - void Basic(); - void Write(int time); + void Basic( ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes); + void WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes, std::shared_ptr input_db, int timestep); private: + std::vector visData; + fillHalo fillData; FILE *TIMELOG; }; #endif diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 852bc194..6fbb627e 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -900,17 +900,13 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){ //if (rank==0) printf("********************************************************\n"); } -void ScaLBL_IonModel::getIonConcentration(int timestep){ - //This function wirte out the data in a normal layout (by aggregating all decomposed domains) - DoubleArray PhaseField(Nx,Ny,Nz); - for (int ic=0; icRegularLayout(Map,&Ci[ic*Np],PhaseField); - ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - IonConcentration_LB_to_Phys(PhaseField); +void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const int ic){ + //This function wirte out the data in a normal layout (by aggregating all decomposed domains) + + ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration); + ScaLBL_DeviceBarrier(); MPI_Barrier(comm); + IonConcentration_LB_to_Phys(IonConcentration); - sprintf(OutputFilename,"Ion%02i_Time_%i.raw",ic+1,timestep); - Mask->AggregateLabels(OutputFilename,PhaseField); - } } void ScaLBL_IonModel::getIonConcentration_debug(int timestep){ diff --git a/models/IonModel.h b/models/IonModel.h index 59382002..4b370978 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -1,6 +1,10 @@ /* * Ion transporte LB Model */ + +#ifndef ScaLBL_IonModel_INC +#define ScaLBL_IonModel_INC + #include #include #include @@ -30,7 +34,7 @@ public: void Create(); void Initialize(); void Run(double *Velocity, double *ElectricField); - void getIonConcentration(int timestep); + void getIonConcentration(DoubleArray &IonConcentration, const int ic); void getIonConcentration_debug(int timestep); void DummyFluidVelocity(); void DummyElectricField(); @@ -95,3 +99,4 @@ private: void AssignIonConcentration_FromFile(double *Ci,const vector &File_ion); void IonConcentration_LB_to_Phys(DoubleArray &Den_reg); }; +#endif diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index c50d8816..b0dde2c7 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -587,38 +587,26 @@ void ScaLBL_Poisson::getElectricPotential_debug(int timestep){ fclose(OUTFILE); } -void ScaLBL_Poisson::getElectricPotential(int timestep){ +void ScaLBL_Poisson::getElectricPotential(DoubleArray &ReturnValues){ //This function wirte out the data in a normal layout (by aggregating all decomposed domains) - DoubleArray PhaseField(Nx,Ny,Nz); //ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField); - ScaLBL_CopyToHost(PhaseField.data(),Psi,sizeof(double)*Nx*Ny*Nz); - ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - - sprintf(OutputFilename,"Electric_Potential_Time_%i.raw",timestep); - Mask->AggregateLabels(OutputFilename,PhaseField); + ScaLBL_CopyToHost(ReturnValues.data(),Psi,sizeof(double)*Nx*Ny*Nz); } -void ScaLBL_Poisson::getElectricField(int timestep){ +void ScaLBL_Poisson::getElectricField(DoubleArray &Values_x, DoubleArray &Values_y, DoubleArray &Values_z){ - DoubleArray PhaseField(Nx,Ny,Nz); - - ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],PhaseField); - ElectricField_LB_to_Phys(PhaseField); + ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],Values_x); + ElectricField_LB_to_Phys(Values_x); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - sprintf(OutputFilename,"ElectricField_X_Time_%i.raw",timestep); - Mask->AggregateLabels(OutputFilename,PhaseField); - ScaLBL_Comm->RegularLayout(Map,&ElectricField[1*Np],PhaseField); - ElectricField_LB_to_Phys(PhaseField); + ScaLBL_Comm->RegularLayout(Map,&ElectricField[1*Np],Values_y); + ElectricField_LB_to_Phys(Values_y); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - sprintf(OutputFilename,"ElectricField_Y_Time_%i.raw",timestep); - Mask->AggregateLabels(OutputFilename,PhaseField); - ScaLBL_Comm->RegularLayout(Map,&ElectricField[2*Np],PhaseField); - ElectricField_LB_to_Phys(PhaseField); + ScaLBL_Comm->RegularLayout(Map,&ElectricField[2*Np],Values_z); + ElectricField_LB_to_Phys(Values_z); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - sprintf(OutputFilename,"ElectricField_Z_Time_%i.raw",timestep); - Mask->AggregateLabels(OutputFilename,PhaseField); + } void ScaLBL_Poisson::getElectricField_debug(int timestep){ diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 241e871a..74abd775 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -16,6 +16,9 @@ #include "analysis/Minkowski.h" #include "ProfilerApp.h" +#ifndef ScaLBL_POISSON_INC +#define ScaLBL_POISSON_INC + class ScaLBL_Poisson{ public: ScaLBL_Poisson(int RANK, int NP, MPI_Comm COMM); @@ -29,9 +32,9 @@ public: void Create(); void Initialize(); void Run(double *ChargeDensity); - void getElectricPotential(int timestep); + void getElectricPotential(DoubleArray &ReturnValues); void getElectricPotential_debug(int timestep); - void getElectricField(int timestep); + void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y, DoubleArray &Values_z); void getElectricField_debug(int timestep); void DummyChargeDensity();//for debugging @@ -96,3 +99,4 @@ private: void getConvergenceLog(int timestep,double error); }; +#endif diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp index a365c46e..086e3633 100644 --- a/models/StokesModel.cpp +++ b/models/StokesModel.cpp @@ -370,25 +370,22 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){ } } -void ScaLBL_StokesModel::getVelocity(int timestep){ +void ScaLBL_StokesModel::getVelocity(DoubleArray &Vel_x, DoubleArray &Vel_y, DoubleArray &Vel_z){ //get velocity in physical unit [m/sec] ScaLBL_D3Q19_Momentum(fq, Velocity, Np); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x); - Velocity_LB_to_Phys(Velocity_x); + ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Vel_x); + Velocity_LB_to_Phys(Vel_x); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - sprintf(OutputFilename,"Velocity_X_Time_%i.raw",timestep); - ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y); - Velocity_LB_to_Phys(Velocity_y); + ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Vel_y); + Velocity_LB_to_Phys(Vel_y); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - sprintf(OutputFilename,"Velocity_Y_Time_%i.raw",timestep); - ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z); - Velocity_LB_to_Phys(Velocity_z); + ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Vel_z); + Velocity_LB_to_Phys(Vel_z); ScaLBL_DeviceBarrier(); MPI_Barrier(comm); - sprintf(OutputFilename,"Velocity_Z_Time_%i.raw",timestep); } void ScaLBL_StokesModel::getVelocity_debug(int timestep){ diff --git a/models/StokesModel.h b/models/StokesModel.h index 6375d4ff..8da373bd 100644 --- a/models/StokesModel.h +++ b/models/StokesModel.h @@ -1,6 +1,9 @@ /* * Multi-relaxation time LBM Model */ +#ifndef ScaLBL_StokesModel_INC +#define ScaLBL_StokesModel_INC + #include #include #include @@ -31,7 +34,7 @@ public: void Run(); void Run_Lite(double *ChargeDensity, double *ElectricField); void VelocityField(); - void getVelocity(int timestep); + void getVelocity(DoubleArray &Velx, DoubleArray &Vel_y, DoubleArray &Vel_z); void getVelocity_debug(int timestep); double CalVelocityConvergence(double& flow_rate_previous,double *ChargeDensity, double *ElectricField); @@ -86,3 +89,4 @@ private: void Velocity_LB_to_Phys(DoubleArray &Vel_reg); vector computeElectricForceAvg(double *ChargeDensity, double *ElectricField); }; +#endif \ No newline at end of file diff --git a/tests/TestIonModel.cpp b/tests/TestIonModel.cpp index 0b57ff1c..2a0a02a9 100644 --- a/tests/TestIonModel.cpp +++ b/tests/TestIonModel.cpp @@ -76,7 +76,7 @@ int main(int argc, char **argv) error = IonModel.CalIonDenConvergence(ci_avg_previous); } } - IonModel.getIonConcentration(timestep); + IonModel.getIonConcentration_debug(timestep); if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); if (rank==0) printf("*************************************************************\n"); diff --git a/tests/TestNernstPlanck.cpp b/tests/TestNernstPlanck.cpp index 96d9c388..def67d5b 100644 --- a/tests/TestNernstPlanck.cpp +++ b/tests/TestNernstPlanck.cpp @@ -87,9 +87,9 @@ int main(int argc, char **argv) } } - PoissonSolver.getElectricPotential(timestep); - PoissonSolver.getElectricField(timestep); - IonModel.getIonConcentration(timestep); + PoissonSolver.getElectricPotential_debug(timestep); + PoissonSolver.getElectricField_debug(timestep); + IonModel.getIonConcentration_debug(timestep); if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); if (rank==0) printf("*************************************************************\n"); diff --git a/tests/TestPNP_Stokes.cpp b/tests/TestPNP_Stokes.cpp index 0e0c3e81..bf05f73c 100644 --- a/tests/TestPNP_Stokes.cpp +++ b/tests/TestPNP_Stokes.cpp @@ -107,10 +107,10 @@ int main(int argc, char **argv) } } - PoissonSolver.getElectricPotential(timestep); - PoissonSolver.getElectricField(timestep); - IonModel.getIonConcentration(timestep); - StokesModel.getVelocity(timestep); + PoissonSolver.getElectricPotential_debug(timestep); + PoissonSolver.getElectricField_debug(timestep); + IonModel.getIonConcentration_debug(timestep); + StokesModel.getVelocity_debug(timestep); if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); if (rank==0) printf("*************************************************************\n"); diff --git a/tests/TestPoissonSolver.cpp b/tests/TestPoissonSolver.cpp index 8753a7b0..32353f65 100644 --- a/tests/TestPoissonSolver.cpp +++ b/tests/TestPoissonSolver.cpp @@ -59,8 +59,8 @@ int main(int argc, char **argv) PoissonSolver.DummyChargeDensity(); PoissonSolver.Run(PoissonSolver.ChargeDensityDummy); - PoissonSolver.getElectricPotential(1); - PoissonSolver.getElectricField(1); + PoissonSolver.getElectricPotential_debug(1); + PoissonSolver.getElectricField_debug(1); if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); if (rank==0) printf("*************************************************************\n"); diff --git a/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp b/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp index 0c1054a2..2b3726a4 100644 --- a/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp +++ b/tests/lbpm_electrokinetic_SingleFluid_simulator.cpp @@ -12,6 +12,7 @@ #include "models/PoissonSolver.h" #include "models/MultiPhysController.h" #include "common/Utilities.h" +#include "analysis/ElectroChemistry.h" using namespace std; @@ -53,7 +54,7 @@ int main(int argc, char **argv) ScaLBL_IonModel IonModel(rank,nprocs,comm); ScaLBL_Poisson PoissonSolver(rank,nprocs,comm); ScaLBL_Multiphys_Controller Study(rank,nprocs,comm);//multiphysics controller coordinating multi-model coupling - + // Load controller information Study.ReadParams(filename); @@ -68,7 +69,10 @@ int main(int argc, char **argv) IonModel.SetDomain(); IonModel.ReadInput(); - IonModel.Create(); + IonModel.Create(); + + // Create analysis object + ElectroChemistryAnalyzer Analysis(IonModel.Dm); // Get internal iteration number StokesModel.timestepMax = Study.getStokesNumIter_PNP_coupling(StokesModel.time_conv,IonModel.time_conv); @@ -96,18 +100,17 @@ int main(int argc, char **argv) timestep++;//AA operations if (timestep%Study.visualization_interval==0){ - PoissonSolver.getElectricPotential(timestep); + Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep); + /* PoissonSolver.getElectricPotential(timestep); PoissonSolver.getElectricField(timestep); IonModel.getIonConcentration(timestep); StokesModel.getVelocity(timestep); + */ } } if (rank==0) printf("Save simulation raw data at maximum timestep\n"); - PoissonSolver.getElectricPotential(timestep); - PoissonSolver.getElectricField(timestep); - IonModel.getIonConcentration(timestep); - StokesModel.getVelocity(timestep); + Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep); if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n"); if (rank==0) printf("*************************************************************\n");