Merge branch 'electrokinetic' of github.com:JamesEMcClure/LBPM-WIA into electrokinetic

common/Domain.cpp

@@ -1448,3 +1448,94 @@ void Domain::ReadFromFile(const std::string& Filename,const std::string& Datatyp
 	//Comm.barrier();
 	MPI_Barrier(Comm);
 }
+
+void Domain::AggregateLabels( const std::string& filename, DoubleArray &UserData ){
+	
+	int nx = Nx;
+	int ny = Ny;
+	int nz = Nz;
+	
+	int npx = nprocx();
+	int npy = nprocy();
+	int npz = nprocz();
+	
+	int ipx = iproc();
+	int ipy = jproc();
+	int ipz = kproc();		
+	
+	int nprocs = nprocx()*nprocy()*nprocz();
+		
+	int full_nx = npx*(nx-2);
+	int full_ny = npy*(ny-2);
+	int full_nz = npz*(nz-2);
+	int local_size = (nx-2)*(ny-2)*(nz-2);
+	unsigned long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
+	
+	double *LocalID;
+	LocalID = new double [local_size];
+		
+	//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
+	// assign the ID for the local sub-region
+	for (int k=1; k<nz-1; k++){
+		for (int j=1; j<ny-1; j++){
+			for (int i=1; i<nx-1; i++){
+				int n = k*nx*ny+j*nx+i;
+				double local_id_val = UserData(i,j,k); 
+				LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
+			}
+		}
+	}
+	MPI_Barrier(Comm);
+
+	// populate the FullID 
+	if (rank() == 0){
+		double *FullID;
+		FullID = new double [full_size];
+		// first handle local ID for rank 0
+		for (int k=1; k<nz-1; k++){
+			for (int j=1; j<ny-1; j++){
+				for (int i=1; i<nx-1; i++){
+					int x = i-1;
+					int y = j-1;
+					int z = k-1;
+					int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
+					unsigned long int n_full = z*long(full_nx)*long(full_ny) + y*long(full_nx) + x;
+					FullID[n_full] = LocalID[n_local];
+				}
+			}
+		}
+		// next get the local ID from the other ranks
+		for (int rnk = 1; rnk<nprocs; rnk++){
+			ipz = rnk / (npx*npy);
+			ipy = (rnk - ipz*npx*npy) / npx;
+			ipx = (rnk - ipz*npx*npy - ipy*npx); 
+			//printf("ipx=%i ipy=%i ipz=%i\n", ipx, ipy, ipz);
+			int tag = 15+rnk;
+			MPI_Recv(LocalID,local_size,MPI_DOUBLE,rnk,tag,Comm,MPI_STATUS_IGNORE);
+			for (int k=1; k<nz-1; k++){
+				for (int j=1; j<ny-1; j++){
+					for (int i=1; i<nx-1; i++){
+						int x = i-1 + ipx*(nx-2);
+						int y = j-1 + ipy*(ny-2);
+						int z = k-1 + ipz*(nz-2);
+						int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
+						unsigned long int n_full = z*long(full_nx)*long(full_ny) + y*long(full_nx) + x;
+						FullID[n_full] = LocalID[n_local];
+					}
+				}
+			}
+		}
+		// write the output
+		FILE *OUTFILE = fopen(filename.c_str(),"wb");
+		fwrite(FullID,8,full_size,OUTFILE);
+		fclose(OUTFILE);
+	}
+	else{
+		// send LocalID to rank=0
+		int tag = 15+ rank();
+		int dstrank = 0;
+		MPI_Send(LocalID,local_size,MPI_DOUBLE,dstrank,tag,Comm);
+	}
+	MPI_Barrier(Comm);
+
+}

common/Domain.h

@@ -198,6 +198,7 @@ public: // Public variables (need to create accessors instead)
     void CommInit(); 
     int PoreCount();
     void AggregateLabels( const std::string& filename );
+    void AggregateLabels( const std::string& filename, DoubleArray &UserData );
 
 private:
 

models/IonModel.cpp

@@ -787,7 +787,20 @@ 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);
 
+        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);
@@ -800,7 +813,6 @@ void ScaLBL_IonModel::getIonConcentration(int timestep){
         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);
@@ -85,6 +86,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/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);    	

tests/TestIonModel.cpp

@@ -9,6 +9,7 @@
 
 #include "models/IonModel.h"
 #include "models/MultiPhysController.h"
+#include "common/Utilities.h"
 
 using namespace std;
 
@@ -18,24 +19,27 @@ using namespace std;
 
 int main(int argc, char **argv)
 {
-    // Initialize MPI
+    // Initialize MPI and error handlers
-    int provided_thread_support = -1;
+    Utilities::startup( argc, argv );
-    MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+
-    MPI_Comm comm;
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+
-    int rank = comm_rank(comm);
+        MPI_Comm comm;
-    int nprocs = comm_size(comm);
+        MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-    if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
+        int rank = comm_rank(comm);
-      std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+        int nprocs = comm_size(comm);
-    }
 
-    // Limit scope so variables that contain communicators will free before MPI_Finialize
-    { 
         if (rank == 0){
             printf("**************************************\n");
             printf("Running Test for Ion Transport \n");
             printf("**************************************\n");
         }
+        // Initialize compute device
+        ScaLBL_SetDevice(rank);
+        ScaLBL_DeviceBarrier();
+        MPI_Barrier(comm);
+
+        PROFILE_ENABLE(1);
         //PROFILE_ENABLE_TRACE();
         //PROFILE_ENABLE_MEMORY();
         PROFILE_SYNCHRONIZE();
@@ -80,10 +84,13 @@ int main(int argc, char **argv)
         PROFILE_STOP("Main");
         PROFILE_SAVE("TestIonModel",1);
         // ****************************************************
+        
         MPI_Barrier(comm);
+        MPI_Comm_free(&comm);
+
     } // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_free(&comm);
+
-    MPI_Finalize();
+    Utilities::shutdown();
 }
 
 

tests/TestNernstPlanck.cpp

@@ -10,6 +10,7 @@
 #include "models/IonModel.h"
 #include "models/PoissonSolver.h"
 #include "models/MultiPhysController.h"
+#include "common/Utilities.h"
 
 using namespace std;
 
@@ -20,23 +21,27 @@ using namespace std;
 int main(int argc, char **argv)
 {
     // Initialize MPI
-    int provided_thread_support = -1;
+    Utilities::startup( argc, argv );
-    MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+
-    MPI_Comm comm;
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+
-    int rank = comm_rank(comm);
+        MPI_Comm comm;
-    int nprocs = comm_size(comm);
+        MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-    if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
+        int rank = comm_rank(comm);
-      std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+        int nprocs = comm_size(comm);
-    }
 
-    // Limit scope so variables that contain communicators will free before MPI_Finialize
-    { 
         if (rank == 0){
             printf("********************************************************\n");
             printf("Running Test for LB-Poisson-Ion Coupling \n");
             printf("********************************************************\n");
         }
+
+        // Initialize compute device
+        ScaLBL_SetDevice(rank);
+        ScaLBL_DeviceBarrier();
+        MPI_Barrier(comm);
+
+        PROFILE_ENABLE(1);
         //PROFILE_ENABLE_TRACE();
         //PROFILE_ENABLE_MEMORY();
         PROFILE_SYNCHRONIZE();
@@ -90,12 +95,14 @@ int main(int argc, char **argv)
         if (rank==0) printf("*************************************************************\n");
 
         PROFILE_STOP("Main");
-        PROFILE_SAVE("lbpm_electrokinetic_simulator",1);
+        PROFILE_SAVE("TestNernstPlanck",1);
         // ****************************************************
+        //
         MPI_Barrier(comm);
+        MPI_Comm_free(&comm);
     } // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_free(&comm);
+
-    MPI_Finalize();
+    Utilities::shutdown();
 }
 
 

tests/TestPNP_Stokes.cpp

@@ -11,6 +11,7 @@
 #include "models/StokesModel.h"
 #include "models/PoissonSolver.h"
 #include "models/MultiPhysController.h"
+#include "common/Utilities.h"
 
 using namespace std;
 
@@ -20,24 +21,27 @@ using namespace std;
 
 int main(int argc, char **argv)
 {
-    // Initialize MPI
+    // Initialize MPI and error handlers
-    int provided_thread_support = -1;
+    Utilities::startup( argc, argv );
-    MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+
-    MPI_Comm comm;
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+
-    int rank = comm_rank(comm);
+        MPI_Comm comm;
-    int nprocs = comm_size(comm);
+        MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-    if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
+        int rank = comm_rank(comm);
-      std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+        int nprocs = comm_size(comm);
-    }
 
-    // Limit scope so variables that contain communicators will free before MPI_Finialize
-    { 
         if (rank == 0){
             printf("********************************************************\n");
             printf("Running Test for LB-Poisson-Ion Coupling \n");
             printf("********************************************************\n");
         }
+        // Initialize compute device
+        ScaLBL_SetDevice(rank);
+        ScaLBL_DeviceBarrier();
+        MPI_Barrier(comm);
+
+        PROFILE_ENABLE(1);
         //PROFILE_ENABLE_TRACE();
         //PROFILE_ENABLE_MEMORY();
         PROFILE_SYNCHRONIZE();
@@ -114,10 +118,13 @@ int main(int argc, char **argv)
         PROFILE_STOP("Main");
         PROFILE_SAVE("TestPNP_Stokes",1);
         // ****************************************************
+        
         MPI_Barrier(comm);
+        MPI_Comm_free(&comm);
+
     } // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_free(&comm);
+
-    MPI_Finalize();
+    Utilities::shutdown();
 }
 
 

tests/TestPoissonSolver.cpp

@@ -8,6 +8,7 @@
 #include <math.h>
 
 #include "models/PoissonSolver.h"
+#include "common/Utilities.h"
 
 using namespace std;
 
@@ -18,23 +19,26 @@ using namespace std;
 int main(int argc, char **argv)
 {
     // Initialize MPI
-    int provided_thread_support = -1;
+    Utilities::startup( argc, argv );
-    MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+    
-    MPI_Comm comm;
+    {// Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+
-    int rank = comm_rank(comm);
+        MPI_Comm comm;
-    int nprocs = comm_size(comm);
+        MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-    if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
+        int rank = comm_rank(comm);
-      std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+        int nprocs = comm_size(comm);
-    }
 
-    // Limit scope so variables that contain communicators will free before MPI_Finialize
-    { 
         if (rank == 0){
             printf("********************************************************\n");
             printf("Running Test for LB-Poisson Solver \n");
             printf("********************************************************\n");
         }
+        // Initialize compute device
+        ScaLBL_SetDevice(rank);
+        ScaLBL_DeviceBarrier();
+        MPI_Barrier(comm);
+
+        PROFILE_ENABLE(1);
         //PROFILE_ENABLE_TRACE();
         //PROFILE_ENABLE_MEMORY();
         PROFILE_SYNCHRONIZE();
@@ -64,10 +68,13 @@ int main(int argc, char **argv)
         PROFILE_STOP("Main");
         PROFILE_SAVE("TestPoissonSolver",1);
         // ****************************************************
+        
         MPI_Barrier(comm);
+        MPI_Comm_free(&comm);
+
     } // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_free(&comm);
+
-    MPI_Finalize();
+    Utilities::shutdown();
 }
 
 

tests/lbpm_electrokinetic_SingleFluid_simulator.cpp

@@ -11,6 +11,7 @@
 #include "models/StokesModel.h"
 #include "models/PoissonSolver.h"
 #include "models/MultiPhysController.h"
+#include "common/Utilities.h"
 
 using namespace std;
 
@@ -20,24 +21,27 @@ using namespace std;
 
 int main(int argc, char **argv)
 {
-    // Initialize MPI
+    // Initialize MPI and error handlers
-    int provided_thread_support = -1;
+    Utilities::startup( argc, argv );
-    MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+
-    MPI_Comm comm;
+    { // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_dup(MPI_COMM_WORLD,&comm);
+
-    int rank = comm_rank(comm);
+        MPI_Comm comm;
-    int nprocs = comm_size(comm);
+        MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-    if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE ){
+        int rank = comm_rank(comm);
-      std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+        int nprocs = comm_size(comm);
-    }
 
-    // Limit scope so variables that contain communicators will free before MPI_Finialize
-    { 
         if (rank == 0){
             printf("********************************************************\n");
             printf("Running LBPM electrokinetic single-fluid solver \n");
             printf("********************************************************\n");
         }
+        // Initialize compute device
+        ScaLBL_SetDevice(rank);
+        ScaLBL_DeviceBarrier();
+        MPI_Barrier(comm);
+
+        PROFILE_ENABLE(1);
         //PROFILE_ENABLE_TRACE();
         //PROFILE_ENABLE_MEMORY();
         PROFILE_SYNCHRONIZE();
@@ -111,10 +115,13 @@ int main(int argc, char **argv)
         PROFILE_STOP("Main");
         PROFILE_SAVE("lbpm_electrokinetic_SingleFluid_simulator",1);
         // ****************************************************
+        
         MPI_Barrier(comm);
+        MPI_Comm_free(&comm);
+
     } // Limit scope so variables that contain communicators will free before MPI_Finialize
-    MPI_Comm_free(&comm);
+
-    MPI_Finalize();
+  Utilities::shutdown();
 }