Updating StackTrace and improving performance converting uCT data

StackTrace/StackTrace.cpp

@@ -7,6 +7,7 @@
 
 #include <algorithm>
 #include <atomic>
+#include <cerrno>
 #include <csignal>
 #include <cstring>
 #include <iostream>
@@ -348,8 +349,11 @@ static inline int exec3( const char *cmd, FUNCTION &fun )
         if ( buffer[0] != 0 )
             fun( buffer );
     }
-    auto status = pclose( pipe );
+    int code = pclose( pipe );
-    int code    = WEXITSTATUS( status );
+    if ( errno == ECHILD ) {
+        errno = 0;
+        code  = 0;
+    }
     std::this_thread::yield(); // Allow any signals to process
     resetSignal( SIGCHLD );    // Clear child exited
     return code;
@@ -1741,7 +1745,7 @@ std::vector<int> StackTrace::defaultSignalsToCatch()
  *  Set the signal handlers                                                  *
  ****************************************************************************/
 static std::function<void( const StackTrace::abort_error &err )> abort_fun;
-static StackTrace::abort_error rethrow()
+StackTrace::abort_error rethrow()
 {
     StackTrace::abort_error error;
 #ifdef USE_LINUX
@@ -1775,14 +1779,14 @@ static StackTrace::abort_error rethrow()
     }
     return error;
 }
-static void term_func_abort( int sig )
+void StackTrace::terminateFunctionSignal( int sig )
 {
     StackTrace::abort_error err;
     err.type      = StackTrace::terminateType::signal;
     err.signal    = sig;
     err.bytes     = StackTrace::Utilities::getMemoryUsage();
     err.stack     = StackTrace::backtrace();
-    err.stackType = StackTrace::printStackType::global;
+    err.stackType = StackTrace::getDefaultStackType();
     abort_fun( err );
 }
 static bool signals_set[256] = { false };
@@ -1829,7 +1833,7 @@ void StackTrace::setErrorHandler( std::function<void( const StackTrace::abort_er
 {
     abort_fun = abort;
     std::set_terminate( term_func );
-    setSignals( defaultSignalsToCatch(), &term_func_abort );
+    setSignals( defaultSignalsToCatch(), &terminateFunctionSignal );
     std::set_unexpected( term_func );
 }
 void StackTrace::clearErrorHandler()
@@ -2215,7 +2219,7 @@ void StackTrace::cleanupStackTrace( multi_stack_info &stack )
             // Remove callstack (and all children) for threads that are just contributing
             bool test = function.find( "_callstack_signal_handler" ) != npos ||
                         function.find( "getGlobalCallStacks" ) != npos ||
-                        function.find( "(" ) == npos;
+                        function.find( "backtrace" ) != npos || function.find( "(" ) == npos;
             if ( test ) {
                 it = stack.children.erase( it );
                 continue;
@@ -2515,3 +2519,11 @@ const char *StackTrace::abort_error::what() const noexcept
             d_msg.erase( i, 1 );
     return d_msg.c_str();
 }
+
+
+/****************************************************************************
+ * Get/Set default stack type                                                *
+ ****************************************************************************/
+static StackTrace::printStackType abort_stackType = StackTrace::printStackType::global;
+void StackTrace::setDefaultStackType( StackTrace::printStackType type ) { abort_stackType = type; }
+StackTrace::printStackType StackTrace::getDefaultStackType() { return abort_stackType; }

StackTrace/StackTrace.h

@@ -246,6 +246,10 @@ void clearSignals();
 void raiseSignal( int signal );
 
 
+//! Default function to abort after catching a signal
+void terminateFunctionSignal( int signal );
+
+
 //! Return a list of all signals that can be caught
 std::vector<int> allSignalsToCatch();
 
@@ -289,6 +293,13 @@ multi_stack_info generateFromString( const std::vector<std::string> &str );
 multi_stack_info generateFromString( const std::string &str );
 
 
+//! Set default stack type
+void setDefaultStackType( StackTrace::printStackType );
+
+//! Get default stack type
+StackTrace::printStackType getDefaultStackType();
+
+
 } // namespace StackTrace
 
 

StackTrace/Utilities.cpp

@@ -8,8 +8,10 @@
 #include <cstring>
 #include <fstream>
 #include <iostream>
+#include <mutex>
 #include <sstream>
 #include <stdexcept>
+#include <typeinfo>
 
 #ifdef USE_MPI
 #include "mpi.h"
@@ -19,6 +21,10 @@
 #include "MemoryApp.h"
 #endif
 
+#ifdef USE_GCOV
+extern "C" void __gcov_flush( void );
+#endif
+
 
 #define perr std::cerr
 
@@ -65,6 +71,12 @@
 // clang-format on
 
 
+#ifdef __GNUC__
+#define USE_ABI
+#include <cxxabi.h>
+#endif
+
+
 namespace StackTrace {
 
 
@@ -96,13 +108,12 @@ inline size_t findfirst( const std::vector<TYPE> &X, TYPE Y )
 /****************************************************************************
  *  Function to terminate the program                                        *
  ****************************************************************************/
-static bool abort_throwException      = false;
+static bool abort_throwException = false;
-static printStackType abort_stackType = printStackType::global;
+static int force_exit            = 0;
-static int force_exit                 = 0;
 void Utilities::setAbortBehavior( bool throwException, int stackType )
 {
     abort_throwException = throwException;
-    abort_stackType      = static_cast<printStackType>( stackType );
+    StackTrace::setDefaultStackType( static_cast<printStackType>( stackType ) );
 }
 void Utilities::abort( const std::string &message, const std::string &filename, const int line )
 {
@@ -112,16 +123,28 @@ void Utilities::abort( const std::string &message, const std::string &filename,
     err.type      = terminateType::abort;
     err.line      = line;
     err.bytes     = Utilities::getMemoryUsage();
-    err.stackType = abort_stackType;
+    err.stackType = StackTrace::getDefaultStackType();
     err.stack     = StackTrace::backtrace();
     throw err;
 }
-static void terminate( const StackTrace::abort_error &err )
+static std::mutex terminate_mutex;
+static inline void callAbort()
 {
+#ifdef USE_GCOV
+    __gcov_flush();
+#endif
+    terminate_mutex.unlock();
+    std::abort();
+}
+void Utilities::terminate( const StackTrace::abort_error &err )
+{
+    // Lock mutex to ensure multiple threads do not try to abort simultaneously
+    terminate_mutex.lock();
+    // Clear the error handlers
     clearErrorHandler();
     // Print the message and abort
     if ( force_exit > 1 ) {
-        std::abort();
+        callAbort();
     } else if ( !abort_throwException ) {
         // Use MPI_abort (will terminate all processes)
         force_exit = 2;
@@ -135,10 +158,11 @@ static void terminate( const StackTrace::abort_error &err )
             MPI_Abort( MPI_COMM_WORLD, -1 );
         }
 #endif
-        std::abort();
+        callAbort();
     } else {
         perr << err.what();
-        std::abort();
+        perr.flush();
+        callAbort();
     }
 }
 
@@ -149,7 +173,7 @@ static void terminate( const StackTrace::abort_error &err )
 static void setTerminateErrorHandler()
 {
     // Set the terminate routine for runtime errors
-    StackTrace::setErrorHandler( terminate );
+    StackTrace::setErrorHandler( Utilities::terminate );
 }
 void Utilities::setErrorHandlers()
 {
@@ -293,4 +317,18 @@ std::string Utilities::exec( const string_view &cmd, int &exit_code )
 }
 
 
+/****************************************************************************
+ *  Get the type name                                                        *
+ ****************************************************************************/
+std::string Utilities::getTypeName( const std::type_info &id )
+{
+    std::string name = id.name();
+#if defined( USE_ABI )
+    int status;
+    name = abi::__cxa_demangle( name.c_str(), 0, 0, &status );
+#endif
+    return name;
+}
+
+
 } // namespace StackTrace

StackTrace/Utilities.h

@@ -4,6 +4,7 @@
 #include <stdexcept>
 #include <string>
 #include <thread>
+#include <typeinfo>
 
 #include "StackTrace/StackTrace.h"
 #include "StackTrace/string_view.h"
@@ -28,9 +29,14 @@ void abort( const std::string &message, const std::string &filename, const int l
 void setAbortBehavior( bool throwException, int stackType = 2 );
 
 
+//! Function to terminate the application
+void terminate( const StackTrace::abort_error &err );
+
+
 //! Function to set the error handlers
 void setErrorHandlers();
 
+
 //! Function to clear the error handlers
 void clearErrorHandlers();
 
@@ -92,6 +98,18 @@ void cause_segfault();
 std::string exec( const StackTrace::string_view &cmd, int &exit_code );
 
 
+//! Return the hopefully demangled name of the given type
+std::string getTypeName( const std::type_info &id );
+
+
+//! Return the hopefully demangled name of the given type
+template<class TYPE>
+inline std::string getTypeName()
+{
+    return getTypeName( typeid( TYPE ) );
+}
+
+
 } // namespace Utilities
 } // namespace StackTrace
 

StackTrace/string_view.h

@@ -119,7 +119,7 @@ public:
         int result = 0;
         for ( int i = 0; i < N && result == 0; i++ )
             if ( d_data[i] != other[i] )
-                result = d_data[i] < other[i] ? -i : i;
+                result = d_data[i] < other[i] ? -( i + 1 ) : ( i + 1 );
         if ( result == 0 )
             result = size() == other.size() ? 0 : size() < other.size() ? -1 : 1;
         return result;

analysis/runAnalysis.cpp

@@ -767,6 +767,8 @@ void runAnalysis::run(int timestep, std::shared_ptr<Database> input_db, TwoPhase
         double *Pressure, double *Velocity, double *fq, double *Den)
 {
     int N = d_N[0]*d_N[1]*d_N[2];
+    NULL_USE( N );
+    NULL_USE( Phi );
     
 	auto db = input_db->getDatabase( "Analysis" );
     //int timestep = db->getWithDefault<int>( "timestep", 0 );
@@ -937,8 +939,6 @@ void runAnalysis::run(int timestep, std::shared_ptr<Database> input_db, TwoPhase
  ******************************************************************/
 void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
 {
-    int N = d_N[0]*d_N[1]*d_N[2];
-
     // Check which analysis steps we need to perform
 	auto color_db =  input_db->getDatabase( "Color" );
 	auto vis_db =  input_db->getDatabase( "Visualization" );
@@ -954,7 +954,7 @@ void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPha
         finish();
     }
 
-    PROFILE_START("run");
+    PROFILE_START("basic");
 
     // Copy the appropriate variables to the host (so we can spawn new threads)
     ScaLBL_DeviceBarrier();
@@ -983,7 +983,6 @@ void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPha
     }
     PROFILE_STOP("Copy data to host");
 
-    PROFILE_START("run",1);
     // Spawn threads to do the analysis work
     //if (timestep%d_restart_interval==0){
     // if ( matches(type,AnalysisType::ComputeAverages) ) {
@@ -1036,12 +1035,11 @@ void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPha
         d_wait_vis = d_tpool.add_work(work);
     }
 
-    PROFILE_STOP("run");
+    PROFILE_STOP("basic");
 }
 
 void runAnalysis::WriteVisData(int timestep, std::shared_ptr<Database> input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
 {
-    int N = d_N[0]*d_N[1]*d_N[2];
 	auto color_db =  input_db->getDatabase( "Color" );
 	auto vis_db =  input_db->getDatabase( "Visualization" );
     //int timestep = color_db->getWithDefault<int>( "timestep", 0 );
@@ -1068,7 +1066,6 @@ void runAnalysis::WriteVisData(int timestep, std::shared_ptr<Database> input_db,
     d_wait_vis = d_tpool.add_work(work2);
 
     //Averages.WriteVis = false;
-   // }
     
     PROFILE_STOP("write vis");
 }

common/Communication.hpp

@@ -44,9 +44,9 @@ Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src
     if ( !src_data.empty() ) {
         int i1[3] = { src_size[0] * src_rank.ix, src_size[1] * src_rank.jy, src_size[2] * src_rank.kz };
         int i2[3] = { i1[0] + src_size[0] - 1, i1[1] + src_size[1] - 1, i1[2] + src_size[2] - 1 };
-        for ( size_t i=0; i<dst_rank.nx; i++ ) {
+        for ( int i=0; i<dst_rank.nx; i++ ) {
-            for ( size_t j=0; j<dst_rank.ny; j++ ) {
+            for ( int j=0; j<dst_rank.ny; j++ ) {
-                for ( size_t k=0; k<dst_rank.nz; k++ ) {
+                for ( int k=0; k<dst_rank.nz; k++ ) {
                     int j1[3] = { i * dst_size[0], j * dst_size[1], k * dst_size[2] };
                     int j2[3] = { j1[0] + dst_size[0] - 1, j1[1] + dst_size[1] - 1, j1[2] + dst_size[2] - 1 };
                     auto index = calcOverlap( i1, i2, j1, j2 );
@@ -65,9 +65,9 @@ Array<TYPE> redistribute( const RankInfoStruct& src_rank, const Array<TYPE>& src
     Array<TYPE> dst_data( dst_size[0], dst_size[1], dst_size[2] );
     int i1[3] = { dst_size[0] * dst_rank.ix, dst_size[1] * dst_rank.jy, dst_size[2] * dst_rank.kz };
     int i2[3] = { i1[0] + dst_size[0] - 1, i1[1] + dst_size[1] - 1, i1[2] + dst_size[2] - 1 };
-    for ( size_t i=0; i<src_rank.nx; i++ ) {
+    for ( int i=0; i<src_rank.nx; i++ ) {
-        for ( size_t j=0; j<src_rank.ny; j++ ) {
+        for ( int j=0; j<src_rank.ny; j++ ) {
-            for ( size_t k=0; k<src_rank.nz; k++ ) {
+            for ( int k=0; k<src_rank.nz; k++ ) {
                 int j1[3] = { i * src_size[0], j * src_size[1], k * src_size[2] };
                 int j2[3] = { j1[0] + src_size[0] - 1, j1[1] + src_size[1] - 1, j1[2] + src_size[2] - 1 };
                 auto index = calcOverlap( i1, i2, j1, j2 );

common/ReadMicroCT.cpp

@@ -70,8 +70,6 @@ Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
     auto n = domain.getVector<int>( "n" );
     int rank = comm_rank(MPI_COMM_WORLD);
     auto nproc = domain.getVector<int>( "nproc" );
-	auto ReadValues = domain.getVector<int>( "ReadValues" );
-	auto WriteValues = domain.getVector<int>( "WriteValues" );
     RankInfoStruct rankInfo( rank, nproc[0], nproc[1], nproc[2] );
     
     // Determine the largest file number to get
@@ -95,29 +93,26 @@ Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
     		ERROR( "Invalid name for first file" );
     	}
     	data = readMicroCT( filename );
-
-    	// Relabel the data
-    	for (int k = 0; k<1024; k++){
-    		for (int j = 0; j<1024; j++){
-    			for (int i = 0; i<1024; i++){
-    				//n = k*Nfx*Nfy + j*Nfx + i;
-    				//char locval = loc_id[n];
-    				char locval = data(i,j,k);
-    				for (int idx=0; idx<ReadValues.size(); idx++){
-    					signed char oldvalue=ReadValues[idx];
-    					signed char newvalue=WriteValues[idx];
-    					if (locval == oldvalue){
-    						data(i,j,k) = newvalue;
-    						idx = ReadValues.size();
-    					}
-    				}
-    			}
-    		}
-    	}
     }
 
     // Redistribute the data
     data = redistribute( srcRankInfo, data, rankInfo, { n[0], n[1], n[2] }, comm );
 
+	// Relabel the data
+    auto ReadValues = domain.getVector<int>( "ReadValues" );
+    auto WriteValues = domain.getVector<int>( "WriteValues" );
+    ASSERT( ReadValues.size() == WriteValues.size() );
+    int readMaxValue = 0;
+    for ( auto v : ReadValues )
+        readMaxValue = std::max( data.max()+1, v );
+    std::vector<int> map( readMaxValue + 1, -1 );
+    for ( size_t i=0; i<ReadValues.size(); i++ )
+        map[ReadValues[i]] = WriteValues[i];
+    for ( size_t i=0; i<data.length(); i++ ) {
+        int t = data(i);
+        ASSERT( t >= 0 && t <= readMaxValue );
+        data(i) = map[t];
+    }
+
     return data;
 }

common/Utilities.cpp

@@ -1,10 +1,116 @@
 #include "common/Utilities.h"
+#include "StackTrace/StackTrace.h"
+#include "StackTrace/ErrorHandlers.h"
+
+#ifdef USE_TIMER
+#include "MemoryApp.h"
+#include "ProfilerApp.h"
+#endif
+
+#ifdef USE_MPI
+#include "mpi.h"
+#endif
 
-#include <math.h>
 #include <algorithm>
+#include <math.h>
+#include <mutex>
 
 
-// Factor a number into it's prime factors
+// Mutex for Utility functions
+static std::mutex Utilities_mutex;
+
+
+/****************************************************************************
+ *  Function to perform the default startup/shutdown sequences               *
+ ****************************************************************************/
+void Utilities::startup( int argc, char **argv )
+{
+    NULL_USE( argc );
+    NULL_USE( argv );
+    // Disable OpenMP
+    Utilities::setenv( "OMP_NUM_THREADS", "1" );
+    Utilities::setenv( "MKL_NUM_THREADS", "1" );
+    // Start MPI
+#ifdef USE_MPI
+    int provided;
+    MPI_Init_thread( &argc, &argv, MPI_THREAD_MULTIPLE, &provided );
+    if ( provided < MPI_THREAD_MULTIPLE ) {
+        int rank;
+        MPI_Comm_rank( MPI_COMM_WORLD, &rank );
+        if ( rank == 0 )
+            std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+    }
+    StackTrace::globalCallStackInitialize( MPI_COMM_WORLD );
+#endif
+    // Set the error handlers
+    Utilities::setAbortBehavior( true, 3 );
+    Utilities::setErrorHandlers();
+}
+void Utilities::shutdown()
+{
+    // Clear the error handlers
+    Utilities::clearErrorHandlers();
+    StackTrace::clearSignals();
+    StackTrace::clearSymbols();
+    int rank = 0;
+#ifdef USE_MPI
+    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
+    StackTrace::globalCallStackFinalize();
+    MPI_Barrier( MPI_COMM_WORLD );
+    MPI_Finalize();
+#endif
+#ifdef USE_TIMER
+    PROFILE_DISABLE();
+    auto memory = MemoryApp::getMemoryStats();
+    if ( rank == 0 && memory.N_new > memory.N_delete )
+        MemoryApp::print( std::cout );
+#endif
+}
+
+
+/****************************************************************************
+ *  Function to set an environemental variable                               *
+ ****************************************************************************/
+void Utilities::setenv( const std::string &name, const std::string &value )
+{
+    Utilities_mutex.lock();
+#if defined( USE_LINUX ) || defined( USE_MAC )
+    bool pass = false;
+    if ( !value.empty() )
+        pass = ::setenv( name.data(), value.data(), 1 ) == 0;
+    else
+        pass = ::unsetenv( name.data() ) == 0;
+#elif defined( USE_WINDOWS )
+    bool pass = SetEnvironmentVariable( name.data(), value.data() ) != 0;
+#else
+#error Unknown OS
+#endif
+    Utilities_mutex.unlock();
+    if ( !pass ) {
+        char msg[1024];
+        if ( !value.empty() )
+            sprintf(
+                msg, "Error setting enviornmental variable: %s=%s\n", name.data(), value.data() );
+        else
+            sprintf( msg, "Error clearing enviornmental variable: %s\n", name.data() );
+        ERROR( msg );
+    }
+}
+std::string Utilities::getenv( const std::string &name )
+{
+    std::string var;
+    Utilities_mutex.lock();
+    auto tmp = std::getenv( name.data() );
+    if ( tmp )
+        var = std::string( tmp );
+    Utilities_mutex.unlock();
+    return var;
+}
+
+
+/****************************************************************************
+ *  Factor a number into it's prime factors                                  *
+ ****************************************************************************/
 std::vector<int> Utilities::factor(size_t number)
 {
     if ( number<=3 ) 
@@ -54,9 +160,13 @@ std::vector<int> Utilities::factor(size_t number)
 }
 
 
-// Dummy function to prevent compiler from optimizing away variable
+/****************************************************************************
+ *  Dummy function to prevent compiler from optimizing away variable         *
+ ****************************************************************************/
 void Utilities::nullUse( void* data )
 {
     NULL_USE(data);
 }
 
+
+

common/Utilities.h

@@ -25,6 +25,37 @@ using StackTrace::Utilities::sleep_ms;
 using StackTrace::Utilities::sleep_s;
 
 
+/*!
+ * \brief Start MPI, error handlers
+ * \details This routine will peform the default startup sequence
+ * \param argc              argc from main
+ * \param argv              argv from main
+ */
+void startup( int argc, char **argv );
+
+/*!
+ * \brief Stop MPI, error handlers
+ * \details This routine will peform the default shutdown sequence to match startup
+ */
+void shutdown();
+
+
+/*!
+ * Get an environmental variable
+ * @param name              The name of the environmental variable
+ * @return                  The value of the enviornmental variable
+ */
+std::string getenv( const std::string &name );
+
+
+/*!
+ * Set an environmental variable
+ * @param name              The name of the environmental variable
+ * @param value             The value to set
+ */
+void setenv( const std::string &name, const std::string &value );
+
+
 //! std::string version of sprintf
 inline std::string stringf( const char *format, ... );
 

tests/lbpm_color_simulator.cpp

@@ -7,6 +7,7 @@
 #include <fstream>
 
 #include "models/ColorModel.h"
+#include "common/Utilities.h"
 
 //#define WRE_SURFACES
 
@@ -15,7 +16,6 @@
  * James E. McClure 2013-2014
  */
 
-using namespace std;
 
 //*************************************************************************
 // Implementation of Two-Phase Immiscible LBM using CUDA
@@ -23,27 +23,26 @@ 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;
-  MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-  int rank = comm_rank(comm);
-  int nprocs = comm_size(comm);
-  if ( rank==0 && provided_thread_support<MPI_THREAD_MULTIPLE )
-    std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
   { // Limit scope so variables that contain communicators will free before MPI_Finialize
 
-	if (rank == 0){
+    MPI_Comm comm;
-		printf("********************************************************\n");
+    MPI_Comm_dup(MPI_COMM_WORLD,&comm);
-		printf("Running Color LBM	\n");
+    int rank = comm_rank(comm);
-		printf("********************************************************\n");
+    int nprocs = comm_size(comm);
-	}
+
-	// Initialize compute device
+    if (rank == 0){
-	int device=ScaLBL_SetDevice(rank);
+	    printf("********************************************************\n");
-	ScaLBL_DeviceBarrier();
+	    printf("Running Color LBM	\n");
-	MPI_Barrier(comm);
+	    printf("********************************************************\n");
-	
+    }
+    // Initialize compute device
+    ScaLBL_SetDevice(rank);
+    ScaLBL_DeviceBarrier();
+    MPI_Barrier(comm);
+
     PROFILE_ENABLE(1);
     //PROFILE_ENABLE_TRACE();
     //PROFILE_ENABLE_MEMORY();
@@ -51,23 +50,26 @@ int main(int argc, char **argv)
     PROFILE_START("Main");
     Utilities::setErrorHandlers();
 
-	auto filename = argv[1];
+    auto filename = argv[1];
-	ScaLBL_ColorModel ColorModel(rank,nprocs,comm);
+    ScaLBL_ColorModel ColorModel(rank,nprocs,comm);
-	ColorModel.ReadParams(filename);
+    ColorModel.ReadParams(filename);
-	ColorModel.SetDomain();    
+    ColorModel.SetDomain();    
-	ColorModel.ReadInput();    
+    ColorModel.ReadInput();    
-	ColorModel.Create();       // creating the model will create data structure to match the pore structure and allocate variables
+    ColorModel.Create();       // creating the model will create data structure to match the pore structure and allocate variables
-	ColorModel.Initialize();   // initializing the model will set initial conditions for variables
+    ColorModel.Initialize();   // initializing the model will set initial conditions for variables
-	ColorModel.Run();	       
+    ColorModel.Run();	       
-	//ColorModel.WriteDebug();
+    //ColorModel.WriteDebug();
-	
+
     PROFILE_STOP("Main");
     PROFILE_SAVE("lbpm_color_simulator",1);
-	// ****************************************************
+    // ****************************************************
-	MPI_Barrier(comm);
+
+    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();
 }