adding silo vis capabilities to electrochem

analysis/ElectroChemistry.cpp

@@ -1,7 +1,11 @@
 #include "analysis/ElectroChemistry.h"
 
 ElectroChemistryAnalyzer::ElectroChemistryAnalyzer(std::shared_ptr <Domain> 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 <Domain> 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 <Domain> 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<Ion.number_ion_species; ion++){
+		Ion.getIonConcentration(Rho,ion);
+	}
+	
 	
 }
 
-void ElectroChemistryAnalyzer::Write(int time){
+void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes, std::shared_ptr<Database> 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<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
+    auto ElectricPotential = std::make_shared<IO::Variable>();
+    auto IonConcentration = std::make_shared<IO::Variable>();
+    auto VxVar = std::make_shared<IO::Variable>();
+    auto VyVar = std::make_shared<IO::Variable>();
+    auto VzVar = std::make_shared<IO::Variable>();
+    
+    if (vis_db->getWithDefault<bool>( "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<bool>( "save_concentration", true )){
+    	for (int ion=0; ion<Ion.number_ion_species; ion++){
+    		sprintf(VisName,"IonConcentration_%i",ion);
+    		IonConcentration->name = 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<bool>( "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<bool>( "save_electric_potential", true )){
+    	ASSERT(visData[0].vars[0]->name=="ElectricPotential");
+    	Poisson.getElectricPotential(ElectricalPotential);
+    	Array<double>& ElectricPotentialData = visData[0].vars[0]->data;
+    	fillData.copy(ElectricalPotential,ElectricPotentialData);
+    }
+
+    if (vis_db->getWithDefault<bool>( "save_concentration", true )){
+    	for (int ion=0; ion<Ion.number_ion_species; ion++){
+    		sprintf(VisName,"IonConcentration_%i",ion);
+    		IonConcentration->name = VisName;
+    		ASSERT(visData[0].vars[1]->name==VisName);
+    		Array<double>& IonConcentrationData = visData[0].vars[1]->data;
+    		Ion.getIonConcentration(Rho,ion);
+    		fillData.copy(Rho,IonConcentrationData);
+    	}
+    }
+
+    if (vis_db->getWithDefault<bool>( "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<double>& VelxData = visData[0].vars[2]->data;
+    	Array<double>& VelyData = visData[0].vars[3]->data;
+    	Array<double>& 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<bool>( "write_silo", true ))
+    	IO::writeData( timestep, visData, comm );
+
+/*    if (vis_db->getWithDefault<bool>( "save_8bit_raw", true )){
+    	char CurrentIDFilename[40];
+    	sprintf(CurrentIDFilename,"id_t%d.raw",timestep);
+    	Averages.AggregateLabels(CurrentIDFilename);
+    }
+*/	
 }

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 <Domain> Dm;
 	double Volume;
 	// input variables
@@ -42,10 +47,12 @@ public:
 	~ElectroChemistryAnalyzer();
 	
 	void SetParams();
-	void Basic();
+	void Basic( ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes);
-	void Write(int time);
+	void WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poisson, ScaLBL_StokesModel &Stokes, std::shared_ptr<Database> input_db, int timestep);
 
 private:
+    std::vector<IO::MeshDataStruct> visData;
+    fillHalo<double> fillData;
 	FILE *TIMELOG;
 };
 #endif

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){
+void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const int ic){
-    //This function wirte out the data in a normal layout (by aggregating all decomposed domains)
+	//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
-    DoubleArray PhaseField(Nx,Ny,Nz);
+
-	for (int ic=0; ic<number_ion_species; ic++){
+	ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration);
-	    ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
+	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
-        ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+	IonConcentration_LB_to_Phys(IonConcentration);
-        IonConcentration_LB_to_Phys(PhaseField);
 
-        sprintf(OutputFilename,"Ion%02i_Time_%i.raw",ic+1,timestep);
-        Mask->AggregateLabels(OutputFilename,PhaseField);
-    }
 }
 
 void ScaLBL_IonModel::getIonConcentration_debug(int timestep){

models/IonModel.h

@@ -1,6 +1,10 @@
 /*
  * Ion transporte LB Model
  */
+
+#ifndef ScaLBL_IonModel_INC
+#define ScaLBL_IonModel_INC
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <sys/stat.h>
@@ -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<std::string> &File_ion);
     void IonConcentration_LB_to_Phys(DoubleArray &Den_reg);
 };
+#endif

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_CopyToHost(ReturnValues.data(),Psi,sizeof(double)*Nx*Ny*Nz);
-    ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
-    
-    sprintf(OutputFilename,"Electric_Potential_Time_%i.raw",timestep);
-    Mask->AggregateLabels(OutputFilename,PhaseField);
 }
 
-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],Values_x);
-
+        ElectricField_LB_to_Phys(Values_x);
-	    ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],PhaseField);
-        ElectricField_LB_to_Phys(PhaseField);
         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);
+	    ScaLBL_Comm->RegularLayout(Map,&ElectricField[1*Np],Values_y);
-        ElectricField_LB_to_Phys(PhaseField);
+        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);
+	    ScaLBL_Comm->RegularLayout(Map,&ElectricField[2*Np],Values_z);
-        ElectricField_LB_to_Phys(PhaseField);
+        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){

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

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);
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Vel_x);
-    Velocity_LB_to_Phys(Velocity_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);
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Vel_y);
-    Velocity_LB_to_Phys(Velocity_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);
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Vel_z);
-    Velocity_LB_to_Phys(Velocity_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){

models/StokesModel.h

@@ -1,6 +1,9 @@
 /*
  * Multi-relaxation time LBM Model
  */
+#ifndef ScaLBL_StokesModel_INC
+#define ScaLBL_StokesModel_INC
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <sys/stat.h>
@@ -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<double> computeElectricForceAvg(double *ChargeDensity, double *ElectricField);
 };
+#endif

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");

tests/TestNernstPlanck.cpp

@@ -87,9 +87,9 @@ int main(int argc, char **argv)
             }
         }
 
-        PoissonSolver.getElectricPotential(timestep);
+        PoissonSolver.getElectricPotential_debug(timestep);
-        PoissonSolver.getElectricField(timestep);
+        PoissonSolver.getElectricField_debug(timestep);
-        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");

tests/TestPNP_Stokes.cpp

@@ -107,10 +107,10 @@ int main(int argc, char **argv)
             }
         }
 
-        PoissonSolver.getElectricPotential(timestep);
+        PoissonSolver.getElectricPotential_debug(timestep);
-        PoissonSolver.getElectricField(timestep);
+        PoissonSolver.getElectricField_debug(timestep);
-        IonModel.getIonConcentration(timestep);
+        IonModel.getIonConcentration_debug(timestep);
-        StokesModel.getVelocity(timestep);
+        StokesModel.getVelocity_debug(timestep);
 
         if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
         if (rank==0) printf("*************************************************************\n");

tests/TestPoissonSolver.cpp

@@ -59,8 +59,8 @@ int main(int argc, char **argv)
         PoissonSolver.DummyChargeDensity();   
 
         PoissonSolver.Run(PoissonSolver.ChargeDensityDummy);
-        PoissonSolver.getElectricPotential(1);
+        PoissonSolver.getElectricPotential_debug(1);
-        PoissonSolver.getElectricField(1);
+        PoissonSolver.getElectricField_debug(1);
 
         if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
         if (rank==0) printf("*************************************************************\n");

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);
+    	Analysis.WriteVis(IonModel,PoissonSolver,StokesModel,Study.db,timestep);
-        PoissonSolver.getElectricField(timestep);
-        IonModel.getIonConcentration(timestep);
-        StokesModel.getVelocity(timestep);
 
         if (rank==0) printf("Maximum timestep is reached and the simulation is completed\n");
         if (rank==0) printf("*************************************************************\n");