save the work;results needs to be validated

models/IonModel.cpp

@@ -787,23 +787,32 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
 }
 
 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; ic<number_ion_species; ic++){
 	    ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
         ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
         IonConcentration_LB_to_Phys(PhaseField);
 
-        //FILE *OUTFILE;
-        //sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
-        //OUTFILE = fopen(LocalRankFilename,"wb");
-        //fwrite(PhaseField.data(),8,N,OUTFILE);
-        //fclose(OUTFILE);
-
         sprintf(OutputFilename,"Ion%02i_Time_%i.raw",ic+1,timestep);
         Mask->AggregateLabels(OutputFilename,PhaseField);
     }
+}
 
+void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
+    //This function write out decomposed data
+    DoubleArray PhaseField(Nx,Ny,Nz);
+	for (int ic=0; ic<number_ion_species; ic++){
+	    ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
+        ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+        IonConcentration_LB_to_Phys(PhaseField);
+
+        FILE *OUTFILE;
+        sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank);
+        OUTFILE = fopen(LocalRankFilename,"wb");
+        fwrite(PhaseField.data(),8,N,OUTFILE);
+        fclose(OUTFILE);
+    }
 }
 
 void ScaLBL_IonModel::IonConcentration_LB_to_Phys(DoubleArray &Den_reg){

models/IonModel.h

@@ -31,6 +31,7 @@ public:
 	void Initialize();
 	void Run(double *Velocity, double *ElectricField);
     void getIonConcentration(int timestep);
+    void getIonConcentration_debug(int timestep);
     void DummyFluidVelocity();
     void DummyElectricField();
     double CalIonDenConvergence(vector<double> &ci_avg_previous);

models/PoissonSolver.cpp

@@ -544,21 +544,55 @@ void ScaLBL_Poisson::DummyChargeDensity(){
 	delete [] ChargeDensity_host;
 }
 
-void ScaLBL_Poisson::getElectricPotential(int timestep){
+void ScaLBL_Poisson::getElectricPotential_debug(int timestep){
+    //This function write out decomposed data
+    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);
+    FILE *OUTFILE;
+    sprintf(LocalRankFilename,"Electric_Potential_Time_%i.%05i.raw",timestep,rank);
+    OUTFILE = fopen(LocalRankFilename,"wb");
+    fwrite(PhaseField.data(),8,N,OUTFILE);
+    fclose(OUTFILE);
+}
 
-        DoubleArray PhaseField(Nx,Ny,Nz);
+void ScaLBL_Poisson::getElectricPotential(int timestep){
-	    //ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
+    //This function wirte out the data in a normal layout (by aggregating all decomposed domains)
-        ScaLBL_CopyToHost(PhaseField.data(),Psi,sizeof(double)*Nx*Ny*Nz);
+    DoubleArray PhaseField(Nx,Ny,Nz);
-        //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+	//ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
-        FILE *OUTFILE;
+    ScaLBL_CopyToHost(PhaseField.data(),Psi,sizeof(double)*Nx*Ny*Nz);
-        sprintf(LocalRankFilename,"Electric_Potential_Time_%i.%05i.raw",timestep,rank);
+    ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
-        OUTFILE = fopen(LocalRankFilename,"wb");
+    
-        fwrite(PhaseField.data(),8,N,OUTFILE);
+    sprintf(OutputFilename,"Electric_Potential_Time_%i.raw",timestep);
-        fclose(OUTFILE);
+    Mask->AggregateLabels(OutputFilename,PhaseField);
 }
 
 void ScaLBL_Poisson::getElectricField(int timestep){
 
+        DoubleArray PhaseField(Nx,Ny,Nz);
+
+	    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);
+        ElectricField_LB_to_Phys(PhaseField);
+        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_DeviceBarrier(); MPI_Barrier(comm);
+        sprintf(OutputFilename,"ElectricField_Z_Time_%i.raw",timestep);
+        Mask->AggregateLabels(OutputFilename,PhaseField);
+}
+
+void ScaLBL_Poisson::getElectricField_debug(int timestep){
+
         //ScaLBL_D3Q7_Poisson_getElectricField(fq,ElectricField,tau,Np);
         //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 

models/PoissonSolver.h

@@ -29,7 +29,9 @@ public:
 	void Initialize();
 	void Run(double *ChargeDensity);
     void getElectricPotential(int timestep);
+    void getElectricPotential_debug(int timestep);
     void getElectricField(int timestep);
+    void getElectricField_debug(int timestep);
     void DummyChargeDensity();//for debugging
 
 	//bool Restart,pBC;
@@ -76,6 +78,7 @@ private:
     char LocalRankString[8];
     char LocalRankFilename[40];
     char LocalRestartFile[40];
+    char OutputFilename[200];
    
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);    	

models/StokesModel.cpp

@@ -375,6 +375,32 @@ void ScaLBL_StokesModel::getVelocity(int timestep){
 	ScaLBL_D3Q19_Momentum(fq, Velocity, Np);
 	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 
+    DoubleArray PhaseField(Nx,Ny,Nz);
+
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
+    Velocity_LB_to_Phys(PhaseField);
+    ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+    sprintf(OutputFilename,"Velocity_X_Time_%i.raw",timestep);
+    Mask->AggregateLabels(OutputFilename,PhaseField);
+
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
+    Velocity_LB_to_Phys(PhaseField);
+    ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+    sprintf(OutputFilename,"Velocity_Y_Time_%i.raw",timestep);
+    Mask->AggregateLabels(OutputFilename,PhaseField);
+
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
+    Velocity_LB_to_Phys(PhaseField);
+    ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+    sprintf(OutputFilename,"Velocity_Z_Time_%i.raw",timestep);
+    Mask->AggregateLabels(OutputFilename,PhaseField);
+}
+
+void ScaLBL_StokesModel::getVelocity_debug(int timestep){
+    //get velocity in physical unit [m/sec]
+	ScaLBL_D3Q19_Momentum(fq, Velocity, Np);
+	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+
     DoubleArray PhaseField(Nx,Ny,Nz);
 	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
     Velocity_LB_to_Phys(PhaseField);

models/StokesModel.h

@@ -32,6 +32,7 @@ public:
 	void Run_Lite(double *ChargeDensity, double *ElectricField);
 	void VelocityField();
     void getVelocity(int timestep);
+    void getVelocity_debug(int timestep);
     double CalVelocityConvergence(double& flow_rate_previous,double *ChargeDensity, double *ElectricField);
 	
 	bool Restart,pBC;
@@ -79,6 +80,7 @@ private:
     char LocalRankString[8];
     char LocalRankFilename[40];
     char LocalRestartFile[40];
+    char OutputFilename[200];
    
     //int rank,nprocs;
     void LoadParams(std::shared_ptr<Database> db0);