diff --git a/IO/Writer.h b/IO/Writer.h
index 08ac0775..710fa0d8 100644
--- a/IO/Writer.h
+++ b/IO/Writer.h
@@ -24,7 +24,7 @@ namespace IO {
  *                          silo - Silo
  * @param[in] append        Append any existing data (default is false)
  */
-void initialize( const std::string& path="", const std::string& format="new", bool append=false );
+void initialize( const std::string& path="", const std::string& format="silo", bool append=false );
 
 
 /*!
diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp
index af8a3915..14a4fddb 100644
--- a/analysis/runAnalysis.cpp
+++ b/analysis/runAnalysis.cpp
@@ -15,23 +15,23 @@
 
 AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs)  
 {
-	lhs = static_cast<AnalysisType>(
-			static_cast<std::underlying_type<AnalysisType>::type>(lhs) |
-			static_cast<std::underlying_type<AnalysisType>::type>(rhs)
-	);
-	return lhs;
+    lhs = static_cast<AnalysisType>(
+          static_cast<std::underlying_type<AnalysisType>::type>(lhs) |
+          static_cast<std::underlying_type<AnalysisType>::type>(rhs)
+    );
+    return lhs;
 }
 bool matches( AnalysisType x, AnalysisType y )
 {
-	return static_cast<std::underlying_type<AnalysisType>::type>(x) &
-			static_cast<std::underlying_type<AnalysisType>::type>(y) != 0;
+    return ( static_cast<std::underlying_type<AnalysisType>::type>(x) &
+             static_cast<std::underlying_type<AnalysisType>::type>(y) ) != 0;
 }
 
 
 template<class TYPE>
 void DeleteArray( const TYPE *p )
 {
-	delete [] p;
+    delete [] p;
 }
 
 
@@ -39,36 +39,36 @@ void DeleteArray( const TYPE *p )
 class WriteRestartWorkItem: public ThreadPool::WorkItemRet<void>
 {
 public:
-	WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_, std::shared_ptr<double> cfq_, int N_ ):
-		filename(filename_), cDen(cDen_), cfq(cfq_), N(N_) {}
-	virtual void run() {
-		PROFILE_START("Save Checkpoint",1);
-		double value;
-		ofstream File(filename,ios::binary);
-		for (int n=0; n<N; n++){
-			// Write the two density values
-			value = cDen.get()[n];
-			File.write((char*) &value, sizeof(value));
-			value = cDen.get()[N+n];
-			File.write((char*) &value, sizeof(value));
-			
-		}
-		for (int n=0; n<N; n++){
-			// Write the distributions
-			for (int q=0; q<19; q++){
-				value = cfq.get()[q*N+n];
-				File.write((char*) &value, sizeof(value));
-			}
-		}
-		File.close();
-		PROFILE_STOP("Save Checkpoint",1);
-	};
+    WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_, std::shared_ptr<double> cfq_, int N_ ):
+        filename(filename_), cfq(cfq_), cDen(cDen_), N(N_) {}
+    virtual void run() {
+        PROFILE_START("Save Checkpoint",1);
+        double value;
+        ofstream File(filename,ios::binary);
+        for (int n=0; n<N; n++){
+            // Write the two density values
+            value = cDen.get()[n];
+            File.write((char*) &value, sizeof(value));
+            value = cDen.get()[N+n];
+            File.write((char*) &value, sizeof(value));
+            
+        }
+        for (int n=0; n<N; n++){
+            // Write the distributions
+            for (int q=0; q<19; q++){
+                value = cfq.get()[q*N+n];
+                File.write((char*) &value, sizeof(value));
+            }
+        }
+        File.close();
+        PROFILE_STOP("Save Checkpoint",1);
+    };
 private:
     WriteRestartWorkItem();
     const char* filename;
     std::shared_ptr<double> cfq,cDen;
-  // const DoubleArray& phase;
-  //const DoubleArray& dist;
+    // const DoubleArray& phase;
+    //const DoubleArray& dist;
     const int N;
 };
 
@@ -78,78 +78,78 @@ static const std::string id_map_filename = "lbpm_id_map.txt";
 class BlobIdentificationWorkItem1: public ThreadPool::WorkItemRet<void>
 {
 public:
-	BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, 
-			std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
-			BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ):
-				timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
-				phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
+    BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, 
+            std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
+            BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ):
+                timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
+                phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
 {
 }
-	~BlobIdentificationWorkItem1() { }
-	virtual void run() {
-		// Compute the global blob id and compare to the previous version
-		PROFILE_START("Identify blobs",1);
-		double vF = 0.0;
-		double vS = -1.0; // one voxel buffer region around solid
-		IntArray& ids = new_index->second;
-		new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm);
-		PROFILE_STOP("Identify blobs",1);
-	}
+    ~BlobIdentificationWorkItem1() { }
+    virtual void run() {
+        // Compute the global blob id and compare to the previous version
+        PROFILE_START("Identify blobs",1);
+        double vF = 0.0;
+        double vS = -1.0; // one voxel buffer region around solid
+        IntArray& ids = new_index->second;
+        new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm);
+        PROFILE_STOP("Identify blobs",1);
+    }
 private:
-	BlobIdentificationWorkItem1();
-	int timestep;
-	int Nx, Ny, Nz;
-	const RankInfoStruct& rank_info;
-	std::shared_ptr<const DoubleArray> phase;
-	const DoubleArray& dist;
-	BlobIDstruct last_id, new_index, new_id;
-	BlobIDList new_list;
-	runAnalysis::commWrapper comm;
+    BlobIdentificationWorkItem1();
+    int timestep;
+    int Nx, Ny, Nz;
+    const RankInfoStruct& rank_info;
+    std::shared_ptr<const DoubleArray> phase;
+    const DoubleArray& dist;
+    BlobIDstruct last_id, new_index, new_id;
+    BlobIDList new_list;
+    runAnalysis::commWrapper comm;
 };
 class BlobIdentificationWorkItem2: public ThreadPool::WorkItemRet<void>
 {
 public:
-	BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, 
-			std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
-			BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ):
-				timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
-				phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
+    BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, 
+            std::shared_ptr<const DoubleArray> phase_, const DoubleArray& dist_,
+            BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ):
+                timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_),
+                phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_))
 {
 }
-	~BlobIdentificationWorkItem2() { }
-	virtual void run() {
-		// Compute the global blob id and compare to the previous version
-		PROFILE_START("Identify blobs maps",1);
-		const IntArray& ids = new_index->second;
-		static int max_id = -1;
-		new_id->first = new_index->first;
-		new_id->second = new_index->second;
-		if ( last_id.get()!=NULL ) {
-			// Compute the timestep-timestep map
-			const IntArray& old_ids = last_id->second;
-			ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm);
-			// Renumber the current timestep's ids
-			getNewIDs(map,max_id,*new_list);
-			renumberIDs(*new_list,new_id->second);
-			writeIDMap(map,timestep,id_map_filename);
-		} else {
-			max_id = -1;
-			ID_map_struct map(new_id->first);
-			getNewIDs(map,max_id,*new_list);
-			writeIDMap(map,timestep,id_map_filename);
-		}
-		PROFILE_STOP("Identify blobs maps",1);
-	}
+    ~BlobIdentificationWorkItem2() { }
+    virtual void run() {
+        // Compute the global blob id and compare to the previous version
+        PROFILE_START("Identify blobs maps",1);
+        const IntArray& ids = new_index->second;
+        static int max_id = -1;
+        new_id->first = new_index->first;
+        new_id->second = new_index->second;
+        if ( last_id.get()!=NULL ) {
+            // Compute the timestep-timestep map
+            const IntArray& old_ids = last_id->second;
+            ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm);
+            // Renumber the current timestep's ids
+            getNewIDs(map,max_id,*new_list);
+            renumberIDs(*new_list,new_id->second);
+            writeIDMap(map,timestep,id_map_filename);
+        } else {
+            max_id = -1;
+            ID_map_struct map(new_id->first);
+            getNewIDs(map,max_id,*new_list);
+            writeIDMap(map,timestep,id_map_filename);
+        }
+        PROFILE_STOP("Identify blobs maps",1);
+    }
 private:
-	BlobIdentificationWorkItem2();
-	int timestep;
-	int Nx, Ny, Nz;
-	const RankInfoStruct& rank_info;
-	std::shared_ptr<const DoubleArray> phase;
-	const DoubleArray& dist;
-	BlobIDstruct last_id, new_index, new_id;
-	BlobIDList new_list;
-	runAnalysis::commWrapper comm;
+    BlobIdentificationWorkItem2();
+    int timestep;
+    int Nx, Ny, Nz;
+    const RankInfoStruct& rank_info;
+    std::shared_ptr<const DoubleArray> phase;
+    const DoubleArray& dist;
+    BlobIDstruct last_id, new_index, new_id;
+    BlobIDList new_list;
+    runAnalysis::commWrapper comm;
 };
 
 
@@ -190,12 +190,12 @@ public:
         PROFILE_STOP("Save Vis",1);
     };
 private:
-	WriteVisWorkItem();
-	int timestep;
-	std::vector<IO::MeshDataStruct>& visData;
-	TwoPhase& Averages;
-	fillHalo<double>& fillData;
-	runAnalysis::commWrapper comm;
+    WriteVisWorkItem();
+    int timestep;
+    std::vector<IO::MeshDataStruct>& visData;
+    TwoPhase& Averages;
+    fillHalo<double>& fillData;
+    runAnalysis::commWrapper comm;
 };
 
 
@@ -204,46 +204,46 @@ private:
 class AnalysisWorkItem: public ThreadPool::WorkItemRet<void>
 {
 public:
-	AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, 
-			BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
-				type(type_), timestep(timestep_), Averages(Averages_), 
-				blob_ids(ids), id_list(id_list_), beta(beta_) { }
-	~AnalysisWorkItem() { }
-	virtual void run() {
-		Averages.NumberComponents_NWP = blob_ids->first;
-		Averages.Label_NWP = blob_ids->second;
-		Averages.Label_NWP_map = *id_list;
-		Averages.NumberComponents_WP = 1;
-		Averages.Label_WP.fill(0.0);
-		if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
-			// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
-		}
-		if ( matches(type,AnalysisType::ComputeAverages) ) {
-			PROFILE_START("Compute dist",1);
-			Averages.Initialize();
-			Averages.ComputeDelPhi();
-			Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
-			Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
-			Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus);
-			Averages.UpdateMeshValues();
-			Averages.ComputeLocal();
-			Averages.Reduce();
-			Averages.PrintAll(timestep);
-			Averages.Initialize();
-			Averages.ComponentAverages();
-			Averages.SortBlobs();
-			Averages.PrintComponents(timestep);
-			PROFILE_STOP("Compute dist",1);
-		}
-	}
+    AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, 
+            BlobIDstruct ids, BlobIDList id_list_, double beta_ ):
+                type(type_), timestep(timestep_), Averages(Averages_), 
+                blob_ids(ids), id_list(id_list_), beta(beta_) { }
+    ~AnalysisWorkItem() { }
+    virtual void run() {
+        Averages.NumberComponents_NWP = blob_ids->first;
+        Averages.Label_NWP = blob_ids->second;
+        Averages.Label_NWP_map = *id_list;
+        Averages.NumberComponents_WP = 1;
+        Averages.Label_WP.fill(0.0);
+        if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
+            // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
+        }
+        if ( matches(type,AnalysisType::ComputeAverages) ) {
+            PROFILE_START("Compute dist",1);
+            Averages.Initialize();
+            Averages.ComputeDelPhi();
+            Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn);
+            Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus);
+            Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus);
+            Averages.UpdateMeshValues();
+            Averages.ComputeLocal();
+            Averages.Reduce();
+            Averages.PrintAll(timestep);
+            Averages.Initialize();
+            Averages.ComponentAverages();
+            Averages.SortBlobs();
+            Averages.PrintComponents(timestep);
+            PROFILE_STOP("Compute dist",1);
+        }
+    }
 private:
-	AnalysisWorkItem();
-	AnalysisType type;
-	int timestep;
-	TwoPhase& Averages;
-	BlobIDstruct blob_ids;
-	BlobIDList id_list;
-	double beta;
+    AnalysisWorkItem();
+    AnalysisType type;
+    int timestep;
+    TwoPhase& Averages;
+    BlobIDstruct blob_ids;
+    BlobIDList id_list;
+    double beta;
 };
 
 
@@ -251,44 +251,44 @@ private:
  *  MPI comm wrapper for use with analysis                         *
  ******************************************************************/
 runAnalysis::commWrapper::commWrapper( int tag_, MPI_Comm comm_, runAnalysis* analysis_ ):
-    		comm(comm_),
-    		tag(tag_),
-    		analysis(analysis_)
+            comm(comm_),
+            tag(tag_),
+            analysis(analysis_)
 {
 }
 runAnalysis::commWrapper::commWrapper( commWrapper &&rhs ):
-    		comm(rhs.comm),
-    		tag(rhs.tag),
-    		analysis(rhs.analysis)
+            comm(rhs.comm),
+            tag(rhs.tag),
+            analysis(rhs.analysis)
 {
-	rhs.tag = -1;
+    rhs.tag = -1;
 }
 runAnalysis::commWrapper::~commWrapper()
 {
-	if ( tag == -1 )
-		return;
-	MPI_Barrier( comm );
-	analysis->d_comm_used[tag] = false;
+    if ( tag == -1 )
+        return;
+    MPI_Barrier( comm );
+    analysis->d_comm_used[tag] = false;
 }
 runAnalysis::commWrapper runAnalysis::getComm( )
 {
-	// Get a tag from root
-	int tag = -1;
-	if ( d_rank == 0 ) {
-		for (int i=0; i<1024; i++) {
-			if ( !d_comm_used[i] ) {
-				tag = i;
-				break;
-			}
-		}
-		if ( tag == -1 )
-			ERROR("Unable to get comm");
-	}
-	MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm );
-	d_comm_used[tag] = true;
-	if ( d_comms[tag] == MPI_COMM_NULL )
-		MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] );
-	return commWrapper(tag,d_comms[tag],this);
+    // Get a tag from root
+    int tag = -1;
+    if ( d_rank == 0 ) {
+        for (int i=0; i<1024; i++) {
+            if ( !d_comm_used[i] ) {
+                tag = i;
+                break;
+            }
+        }
+        if ( tag == -1 )
+            ERROR("Unable to get comm");
+    }
+    MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm );
+    d_comm_used[tag] = true;
+    if ( d_comms[tag] == MPI_COMM_NULL )
+        MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] );
+    return commWrapper(tag,d_comms[tag],this);
 }
 
 
@@ -296,116 +296,122 @@ runAnalysis::commWrapper runAnalysis::getComm( )
  *  Constructor/Destructors                                        *
  ******************************************************************/
 runAnalysis::runAnalysis( std::shared_ptr<Database> db,
-		const RankInfoStruct& rank_info, std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> Dm,
-		int Np, bool Regular, double beta, IntArray Map ):
-    		d_Np( Np ),
-    		d_beta( beta ),
-    		d_regular ( Regular),
-    		d_rank_info( rank_info ),
-    		d_Map( Map ),
-    		d_ScaLBL_Comm( ScaLBL_Comm),
-    		d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1)
+        const RankInfoStruct& rank_info, std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> Dm,
+        int Np, bool Regular, double beta, IntArray Map ):
+            d_Np( Np ),
+            d_beta( beta ),
+            d_regular ( Regular),
+            d_rank_info( rank_info ),
+            d_Map( Map ),
+            d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1),
+            d_ScaLBL_Comm( ScaLBL_Comm)
 {
 
-	char rankString[20];
-	sprintf(rankString,"%05d",Dm->rank());
-	d_N[0] = Dm->Nx;
-	d_N[1] = Dm->Ny;
-	d_N[2] = Dm->Nz;
-	d_restart_interval = db->getScalar<int>( "restart_interval" );
-	d_analysis_interval = db->getScalar<int>( "analysis_interval" );
-	d_blobid_interval = db->getScalar<int>( "blobid_interval" );
-	d_visualization_interval = db->getScalar<int>( "visualization_interval" );
-	auto restart_file = db->getScalar<std::string>( "restart_file" );
-	d_restartFile = restart_file + "." + rankString;
-	d_rank = MPI_WORLD_RANK();
-	writeIDMap(ID_map_struct(),0,id_map_filename);
-	// Initialize IO for silo
-	IO::initialize("","silo","false");
-	// Create the MeshDataStruct	
-	d_meshData.resize(1);
-	d_meshData[0].meshName = "domain";
-	d_meshData[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 PhaseVar = std::make_shared<IO::Variable>();
-	auto PressVar = 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>();
-	auto SignDistVar = std::make_shared<IO::Variable>();
-	auto BlobIDVar = std::make_shared<IO::Variable>();
-	
-	PhaseVar->name = "phase";
-	PhaseVar->type = IO::VariableType::VolumeVariable;
-	PhaseVar->dim = 1;
-	PhaseVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
-	d_meshData[0].vars.push_back(PhaseVar);
-	PressVar->name = "Pressure";
-	PressVar->type = IO::VariableType::VolumeVariable;
-	PressVar->dim = 1;
-	PressVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
-	d_meshData[0].vars.push_back(PressVar);
-	
-	VxVar->name = "Velocity_x";
-	VxVar->type = IO::VariableType::VolumeVariable;
-	VxVar->dim = 1;
-	VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
-	d_meshData[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);
-	d_meshData[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);
-	d_meshData[0].vars.push_back(VzVar);
-	
-	SignDistVar->name = "SignDist";
-	SignDistVar->type = IO::VariableType::VolumeVariable;
-	SignDistVar->dim = 1;
-	SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
-	d_meshData[0].vars.push_back(SignDistVar);
-	BlobIDVar->name = "BlobID";
-	BlobIDVar->type = IO::VariableType::VolumeVariable;
-	BlobIDVar->dim = 1;
-	BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
-	d_meshData[0].vars.push_back(BlobIDVar);
-	// Initialize the comms
-	MPI_Comm_dup(MPI_COMM_WORLD,&d_comm);
-	for (int i=0; i<1024; i++) {
-		d_comms[i] = MPI_COMM_NULL;
-		d_comm_used[i] = false;
-	}
-	// Initialize the threads
-	int N_threads = db->getWithDefault<int>( "N_threads", 4 );
-	auto method = db->getWithDefault<std::string>( "load_balance", "default" );
-	createThreads( method, N_threads );
+    // Ids of work items to use for dependencies
+    ThreadPool::thread_id_t d_wait_blobID;
+    ThreadPool::thread_id_t d_wait_analysis;
+    ThreadPool::thread_id_t d_wait_vis;
+    ThreadPool::thread_id_t d_wait_restart;
+
+    char rankString[20];
+    sprintf(rankString,"%05d",Dm->rank());
+    d_N[0] = Dm->Nx;
+    d_N[1] = Dm->Ny;
+    d_N[2] = Dm->Nz;
+    d_restart_interval = db->getScalar<int>( "restart_interval" );
+    d_analysis_interval = db->getScalar<int>( "analysis_interval" );
+    d_blobid_interval = db->getScalar<int>( "blobid_interval" );
+    d_visualization_interval = db->getScalar<int>( "visualization_interval" );
+    auto restart_file = db->getScalar<std::string>( "restart_file" );
+    d_restartFile = restart_file + "." + rankString;
+    d_rank = MPI_WORLD_RANK();
+    writeIDMap(ID_map_struct(),0,id_map_filename);
+    // Initialize IO for silo
+    IO::initialize("","silo","false");
+    // Create the MeshDataStruct    
+    d_meshData.resize(1);
+    d_meshData[0].meshName = "domain";
+    d_meshData[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 PhaseVar = std::make_shared<IO::Variable>();
+    auto PressVar = 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>();
+    auto SignDistVar = std::make_shared<IO::Variable>();
+    auto BlobIDVar = std::make_shared<IO::Variable>();
+    
+    PhaseVar->name = "phase";
+    PhaseVar->type = IO::VariableType::VolumeVariable;
+    PhaseVar->dim = 1;
+    PhaseVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    d_meshData[0].vars.push_back(PhaseVar);
+    PressVar->name = "Pressure";
+    PressVar->type = IO::VariableType::VolumeVariable;
+    PressVar->dim = 1;
+    PressVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    d_meshData[0].vars.push_back(PressVar);
+    
+    VxVar->name = "Velocity_x";
+    VxVar->type = IO::VariableType::VolumeVariable;
+    VxVar->dim = 1;
+    VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    d_meshData[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);
+    d_meshData[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);
+    d_meshData[0].vars.push_back(VzVar);
+    
+    SignDistVar->name = "SignDist";
+    SignDistVar->type = IO::VariableType::VolumeVariable;
+    SignDistVar->dim = 1;
+    SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    d_meshData[0].vars.push_back(SignDistVar);
+    BlobIDVar->name = "BlobID";
+    BlobIDVar->type = IO::VariableType::VolumeVariable;
+    BlobIDVar->dim = 1;
+    BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+    d_meshData[0].vars.push_back(BlobIDVar);
+    // Initialize the comms
+    MPI_Comm_dup(MPI_COMM_WORLD,&d_comm);
+    for (int i=0; i<1024; i++) {
+        d_comms[i] = MPI_COMM_NULL;
+        d_comm_used[i] = false;
+    }
+    // Initialize the threads
+    int N_threads = db->getWithDefault<int>( "N_threads", 4 );
+    auto method = db->getWithDefault<std::string>( "load_balance", "default" );
+    createThreads( method, N_threads );
 }
 runAnalysis::~runAnalysis( )
 {
-	// Finish processing analysis
-	finish();
-	// Clear internal data
-	MPI_Comm_free( &d_comm );
-	for (int i=0; i<1024; i++) {
-		if ( d_comms[i] != MPI_COMM_NULL )
-			MPI_Comm_free(&d_comms[i]);
-	}
+    // Finish processing analysis
+    finish();
+    // Clear internal data
+    MPI_Comm_free( &d_comm );
+    for (int i=0; i<1024; i++) {
+        if ( d_comms[i] != MPI_COMM_NULL )
+            MPI_Comm_free(&d_comms[i]);
+    }
 }
 void runAnalysis::finish( )
 {
-	PROFILE_START("finish");
-	// Wait for the work items to finish
-	d_tpool.wait_pool_finished();
-	// Clear the wait ids
-	d_wait_blobID.reset();
-	d_wait_analysis.reset();
-	d_wait_vis.reset();
-	d_wait_restart.reset();
-	// Syncronize
-	MPI_Barrier( d_comm );
-	PROFILE_STOP("finish");
+    PROFILE_START("finish");
+    // Wait for the work items to finish
+    d_tpool.wait_pool_finished();
+    // Clear the wait ids
+    d_wait_blobID.reset();
+    d_wait_analysis.reset();
+    d_wait_vis.reset();
+    d_wait_restart.reset();
+    // Syncronize
+    MPI_Barrier( d_comm );
+    PROFILE_STOP("finish");
 }
 
 
@@ -414,50 +420,50 @@ void runAnalysis::finish( )
  ******************************************************************/
 void print( const std::vector<int>& ids )
 {
-	if ( ids.empty() )
-		return;
-	printf("%i",ids[0]);
-	for (size_t i=1; i<ids.size(); i++)
-		printf(", %i",ids[i]);
-	printf("\n");
+    if ( ids.empty() )
+        return;
+    printf("%i",ids[0]);
+    for (size_t i=1; i<ids.size(); i++)
+        printf(", %i",ids[i]);
+    printf("\n");
 }
 void runAnalysis::createThreads( const std::string& method, int N_threads )
 {
-	// Check if we are not using analysis threads
-	if ( method == "none" )
-		return;
-	// Check if we have thread support
-	int thread_support;
-	MPI_Query_thread( &thread_support );
-	if ( thread_support < MPI_THREAD_MULTIPLE ) {
-		std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
-		return;
-	}
-	// Create the threads
-	const auto cores = d_tpool.getProcessAffinity();
-	if ( cores.empty() ) {
-		// We were not able to get the cores for the process
-		d_tpool.setNumThreads( N_threads );
-	} else if ( method == "default" ) {
-		// Create the given number of threads, but let the OS manage affinities
-		d_tpool.setNumThreads( N_threads );
-	} else if ( method == "independent" ) {
-		int N = cores.size() - 1;
-		d_tpool.setNumThreads( N );
-		d_tpool.setThreadAffinity( { cores[0] } );
-		for ( int i=0; i<N; i++)
-			d_tpool.setThreadAffinity( i, { cores[i+1] } );
-	}
-	// Print the current affinities
-	if ( d_rank == 0 ) {
-		printf("Affinities - rank 0:\n");
-		printf("Main: ");
-		print(d_tpool.getProcessAffinity());
-		for (int i=0; i<d_tpool.getNumThreads(); i++) {
-			printf("Thread %i: ",i+1);
-			print(d_tpool.getThreadAffinity(i));
-		}
-	}
+    // Check if we are not using analysis threads
+    if ( method == "none" )
+        return;
+    // Check if we have thread support
+    int thread_support;
+    MPI_Query_thread( &thread_support );
+    if ( thread_support < MPI_THREAD_MULTIPLE ) {
+        std::cerr << "Warning: Failed to start MPI with necessary thread support, thread support will be disabled" << std::endl;
+        return;
+    }
+    // Create the threads
+    const auto cores = d_tpool.getProcessAffinity();
+    if ( cores.empty() ) {
+        // We were not able to get the cores for the process
+        d_tpool.setNumThreads( N_threads );
+    } else if ( method == "default" ) {
+        // Create the given number of threads, but let the OS manage affinities
+        d_tpool.setNumThreads( N_threads );
+    } else if ( method == "independent" ) {
+        int N = cores.size() - 1;
+        d_tpool.setNumThreads( N );
+        d_tpool.setThreadAffinity( { cores[0] } );
+        for ( int i=0; i<N; i++)
+            d_tpool.setThreadAffinity( i, { cores[i+1] } );
+    }
+    // Print the current affinities
+    if ( d_rank == 0 ) {
+        printf("Affinities - rank 0:\n");
+        printf("Main: ");
+        print(d_tpool.getProcessAffinity());
+        for (int i=0; i<d_tpool.getNumThreads(); i++) {
+            printf("Thread %i: ",i+1);
+            print(d_tpool.getThreadAffinity(i));
+        }
+    }
 }
 
 
@@ -466,17 +472,17 @@ void runAnalysis::createThreads( const std::string& method, int N_threads )
  ******************************************************************/
 AnalysisType runAnalysis::computeAnalysisType( int timestep )
 {
-	AnalysisType type = AnalysisType::AnalyzeNone;
-	if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
-		// Copy the phase indicator field for the earlier timestep
-		// printf("Copy phase indicator,timestep=%i\n",timestep);
-		type |= AnalysisType::CopyPhaseIndicator;
-	}
-	if ( timestep%d_blobid_interval == 0 ) {
-		// Identify blobs and update global ids in time
-		type |= AnalysisType::IdentifyBlobs;
-	}
-	/*#ifdef USE_CUDA
+    AnalysisType type = AnalysisType::AnalyzeNone;
+    if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
+        // Copy the phase indicator field for the earlier timestep
+        // printf("Copy phase indicator,timestep=%i\n",timestep);
+        type |= AnalysisType::CopyPhaseIndicator;
+    }
+    if ( timestep%d_blobid_interval == 0 ) {
+        // Identify blobs and update global ids in time
+        type |= AnalysisType::IdentifyBlobs;
+    }
+    /*#ifdef USE_CUDA
         if ( tpool.getQueueSize()<=3 && tpool.getNumThreads()>0 && timestep%50==0 ) {
             // Keep a few blob identifications queued up to keep the processors busy,
             // allowing us to track the blobs as fast as possible
@@ -484,28 +490,28 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep )
             type |= AnalysisType::IdentifyBlobs;
         }
     #endif */
-	if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
-		// Copy the averages to the CPU (and identify blobs)
-		//printf("Copy sim state, timestep=%i \n",timestep);
-		type |= AnalysisType::CopySimState;
-		type |= AnalysisType::IdentifyBlobs;
-	}
-	if ( timestep%d_analysis_interval == 0 ) {
-		// Run the analysis
-		//printf("Compute averages, timestep=%i \n",timestep);
-		type |= AnalysisType::ComputeAverages;
-	}
-	if (timestep%d_restart_interval == 0) {
-		// Write the restart file
-		type |= AnalysisType::CreateRestart;
-	}
-	if (timestep%d_visualization_interval == 0) {
-		// Write the visualization data
-		type |= AnalysisType::WriteVis;
-		type |= AnalysisType::CopySimState;
-		type |= AnalysisType::IdentifyBlobs;
-	}
-	return type;
+    if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
+        // Copy the averages to the CPU (and identify blobs)
+        //printf("Copy sim state, timestep=%i \n",timestep);
+        type |= AnalysisType::CopySimState;
+        type |= AnalysisType::IdentifyBlobs;
+    }
+    if ( timestep%d_analysis_interval == 0 ) {
+        // Run the analysis
+        //printf("Compute averages, timestep=%i \n",timestep);
+        type |= AnalysisType::ComputeAverages;
+    }
+    if (timestep%d_restart_interval == 0) {
+        // Write the restart file
+        type |= AnalysisType::CreateRestart;
+    }
+    if (timestep%d_visualization_interval == 0) {
+        // Write the visualization data
+        type |= AnalysisType::WriteVis;
+        type |= AnalysisType::CopySimState;
+        type |= AnalysisType::IdentifyBlobs;
+    }
+    return type;
 }
 
 
@@ -514,56 +520,56 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep )
  *  Run the analysis                                               *
  ******************************************************************/
 void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
-		double *Pressure, double *Velocity, double *fq, double *Den)
+        double *Pressure, double *Velocity, double *fq, double *Den)
 {
-	int N = d_N[0]*d_N[1]*d_N[2];
+    int N = d_N[0]*d_N[1]*d_N[2];
 
-	// Check which analysis steps we need to perform
-	auto type = computeAnalysisType( timestep );
-	if ( type == AnalysisType::AnalyzeNone )
-		return;
+    // Check which analysis steps we need to perform
+    auto type = computeAnalysisType( timestep );
+    if ( type == AnalysisType::AnalyzeNone )
+        return;
 
-	// Check how may queued items we have
-	if ( d_tpool.N_queued() > 20 ) {
-		std::cerr << "Analysis queue is getting behind, waiting ...\n";
-		finish();
-	}
+    // Check how may queued items we have
+    if ( d_tpool.N_queued() > 20 ) {
+        std::cerr << "Analysis queue is getting behind, waiting ...\n";
+        finish();
+    }
 
-	PROFILE_START("run");
+    PROFILE_START("run");
 
-	// Copy the appropriate variables to the host (so we can spawn new threads)
-	ScaLBL_DeviceBarrier();
-	PROFILE_START("Copy data to host",1);
-	std::shared_ptr<DoubleArray> phase;
-	/*   if ( matches(type,AnalysisType::CopyPhaseIndicator) ||
-    		matches(type,AnalysisType::ComputeAverages) ||
-    		matches(type,AnalysisType::CopySimState) || 
-    		matches(type,AnalysisType::IdentifyBlobs) )
+    // Copy the appropriate variables to the host (so we can spawn new threads)
+    ScaLBL_DeviceBarrier();
+    PROFILE_START("Copy data to host",1);
+    std::shared_ptr<DoubleArray> phase;
+    /*   if ( matches(type,AnalysisType::CopyPhaseIndicator) ||
+            matches(type,AnalysisType::ComputeAverages) ||
+            matches(type,AnalysisType::CopySimState) || 
+            matches(type,AnalysisType::IdentifyBlobs) )
     {
-    	phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
-    	//ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
-    	//   try 2 d_ScaLBL_Comm.RegulLayout(d_Map,Phi,Averages.Phase);
-    	//   memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
-    	int Nx = d_N[0];
-    	int Ny = d_N[1];
-    	int Nz = d_N[2];
-    	double *TmpDat;
-    	TmpDat = new double [d_Np];
-    	ScaLBL_CopyToHost(&TmpDat[0],&Phi[0], d_Np*sizeof(double));
-    	for (int k=0; k<Nz; k++){
-    		for (int j=0; j<Ny; j++){
-    			for (int i=0; i<Nx; i++){
-    				int n=k*Nx*Ny+j*Nx+i;
-    				int idx=d_Map(i,j,k);
-    				if (!(idx<0)){
-    					double value=TmpDat[idx];
-    					//regdata(i,j,k)=value;
-    					phase->data()[n]=value;
-    				}
-    			}
-    		}
-    	}
-    	delete [] TmpDat;
+        phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
+        //ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
+        //   try 2 d_ScaLBL_Comm.RegulLayout(d_Map,Phi,Averages.Phase);
+        //   memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
+        int Nx = d_N[0];
+        int Ny = d_N[1];
+        int Nz = d_N[2];
+        double *TmpDat;
+        TmpDat = new double [d_Np];
+        ScaLBL_CopyToHost(&TmpDat[0],&Phi[0], d_Np*sizeof(double));
+        for (int k=0; k<Nz; k++){
+            for (int j=0; j<Ny; j++){
+                for (int i=0; i<Nx; i++){
+                    int n=k*Nx*Ny+j*Nx+i;
+                    int idx=d_Map(i,j,k);
+                    if (!(idx<0)){
+                        double value=TmpDat[idx];
+                        //regdata(i,j,k)=value;
+                        phase->data()[n]=value;
+                    }
+                }
+            }
+        }
+        delete [] TmpDat;
     }
     */
     //if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
@@ -571,57 +577,57 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
       if (d_regular)
         d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus);
       else 
-	ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
+    ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
         //memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
     }
     if ( timestep%d_analysis_interval == 0 ) {
       if (d_regular)
         d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus);
       else 
-	ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double));
+    ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double));
         //memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
     }
     //if ( matches(type,AnalysisType::CopySimState) ) {
     if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
-    	// Copy the members of Averages to the cpu (phase was copied above)
-    	PROFILE_START("Copy-Pressure",1);
-    	ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
-    	//ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
-    	ScaLBL_DeviceBarrier();
-    	PROFILE_STOP("Copy-Pressure",1);
-    	PROFILE_START("Copy-Wait",1);
-    	PROFILE_STOP("Copy-Wait",1);
-    	PROFILE_START("Copy-State",1);
-    	//memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
-    	if (d_regular)
-    		d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
-    	else
-    		ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
-    	// copy other variables
-    	d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
-    	d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
-    	d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
-    	d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
-    	PROFILE_STOP("Copy-State",1);
+        // Copy the members of Averages to the cpu (phase was copied above)
+        PROFILE_START("Copy-Pressure",1);
+        ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
+        //ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
+        ScaLBL_DeviceBarrier();
+        PROFILE_STOP("Copy-Pressure",1);
+        PROFILE_START("Copy-Wait",1);
+        PROFILE_STOP("Copy-Wait",1);
+        PROFILE_START("Copy-State",1);
+        //memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
+        if (d_regular)
+            d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
+        else
+            ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
+        // copy other variables
+        d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
+        d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
+        d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
+        d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
+        PROFILE_STOP("Copy-State",1);
     }
     std::shared_ptr<double> cfq,cDen;
     //if ( matches(type,AnalysisType::CreateRestart) ) {
     if (timestep%d_restart_interval==0){
-    	// Copy restart data to the CPU
-    	cDen = std::shared_ptr<double>(new double[2*d_Np],DeleteArray<double>);
-    	cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
-    	ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
-    	ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double));
+        // Copy restart data to the CPU
+        cDen = std::shared_ptr<double>(new double[2*d_Np],DeleteArray<double>);
+        cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
+        ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
+        ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double));
     }
     PROFILE_STOP("Copy data to host",1);
 
     // Spawn threads to do blob identification work
     if ( matches(type,AnalysisType::IdentifyBlobs) ) {
-    	phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
-    	if (d_regular)
+        phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
+        if (d_regular)
             d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase);
-    	else
-	  ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
+        else
+      ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
 
         BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
         BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
@@ -642,11 +648,11 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
     //if (timestep%d_restart_interval==0){
     // if ( matches(type,AnalysisType::ComputeAverages) ) {
     if ( timestep%d_analysis_interval == 0 ) {
-    	auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta);
-    	work->add_dependency(d_wait_blobID);
-    	work->add_dependency(d_wait_analysis);
-    	work->add_dependency(d_wait_vis);     // Make sure we are done using analysis before modifying
-    	d_wait_analysis = d_tpool.add_work(work);
+        auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta);
+        work->add_dependency(d_wait_blobID);
+        work->add_dependency(d_wait_analysis);
+        work->add_dependency(d_wait_vis);     // Make sure we are done using analysis before modifying
+        d_wait_analysis = d_tpool.add_work(work);
     }
 
     // Spawn a thread to write the restart file
diff --git a/common/Array.hpp b/common/Array.hpp
index f7d4398c..cdf58125 100644
--- a/common/Array.hpp
+++ b/common/Array.hpp
@@ -1002,7 +1002,7 @@ Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::coarsen(
             throw std::invalid_argument( "Array must be multiple of filter size" );
     }
     Array<TYPE, FUN, Allocator> y( S2 );
-    if ( d_size.ndim() <= 3 )
+    if ( d_size.ndim() > 3 )
         throw std::logic_error( "Function programmed for more than 3 dimensions" );
     const auto &Nh = filter.d_size;
     for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) {
diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp
index ccb0f5af..80dbb52b 100644
--- a/tests/lbpm_uCT_pp.cpp
+++ b/tests/lbpm_uCT_pp.cpp
@@ -30,463 +30,415 @@
 int main(int argc, char **argv)
 {
 
-	// Initialize MPI
-	int rank, nprocs;
-	MPI_Init(&argc,&argv);
-	MPI_Comm comm = MPI_COMM_WORLD;
-	MPI_Comm_rank(comm,&rank);
-	MPI_Comm_size(comm,&nprocs);
-	{
-		Utilities::setErrorHandlers();
-		PROFILE_START("Main");
-
-		//std::vector<std::string> filenames;
-		if ( argc<2 ) {
-			if ( rank == 0 ){
-				printf("At least one filename must be specified\n");
-			}
-			return 1;
-		}
-		std::string filename = std::string(argv[1]);
-		if ( rank == 0 ){
-			printf("Input data file: %s\n",filename.c_str());
-		}
-
-		auto db = std::make_shared<Database>( filename );
-		auto domain_db = db->getDatabase( "Domain" );
-		auto uct_db = db->getDatabase( "uCT" );
-		auto analysis_db = db->getDatabase( "Analysis" );
-
-		// Read domain values
-		auto L = domain_db->getVector<double>( "L" );
-		auto size = domain_db->getVector<int>( "n" );
-		auto nproc = domain_db->getVector<int>( "nproc" );
-		int BoundaryCondition = domain_db->getScalar<int>( "BC" );
-		int nx = size[0];
-		int ny = size[1];
-		int nz = size[2];
-		double Lx = L[0];
-		double Ly = L[1];
-		double Lz = L[2];
-		int nprocx = nproc[0];
-		int nprocy = nproc[1];
-		int nprocz = nproc[2];
-
-		auto InputFile=uct_db->getScalar<std::string>( "InputFile" );
-		auto target=uct_db->getScalar<int>("target");
-		auto background=uct_db->getScalar<int>("background");
-		auto rough_cutoff=uct_db->getScalar<float>( "rough_cutoff" );    
-		auto lamda=uct_db->getScalar<float>( "lamda" );    
-		auto nlm_sigsq=uct_db->getScalar<float>( "nlm_sigsq" );    
-		auto nlm_depth=uct_db->getScalar<int>( "nlm_depth" );    
-		auto cx=uct_db->getScalar<int>( "center_x" );
-		auto cy=uct_db->getScalar<int>( "center_y" );
-		auto cz=uct_db->getScalar<int>( "center_z" );
-		auto CylRad=uct_db->getScalar<float>( "cylinder_radius" );
-
-		//.......................................................................
-		// Reading the domain information file
-		//.......................................................................
-		//    std::shared_ptr<Domain> Dm ();
-		//for (int i=0; i<Dm->Nx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = 1;
-		//Dm->CommInit();
-
-		// Check that the number of processors >= the number of ranks
-		if ( rank==0 ) {
-			printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
-			printf("Number of MPI ranks used: %i \n", nprocs);
-			printf("Full domain size: %i x %i x %i  \n",nx*nprocx,ny*nprocy,nz*nprocz);
-		}
-		if ( nprocs < nprocx*nprocy*nprocz ){
-			ERROR("Insufficient number of processors");
-		}
-
-		// Determine the maximum number of levels for the desired coarsen ratio
-		int ratio[3] = {2,2,2};
-		//std::vector<size_t> ratio = {4,4,4};
-		// need to set up databases for each level of the mesh
-		std:vector<Database> multidomain_db;
-
-		std::vector<int> Nx(1,nx), Ny(1,ny), Nz(1,nz);
-		while ( Nx.back()%ratio[0]==0 && Nx.back()>8 &&
-				Ny.back()%ratio[1]==0 && Ny.back()>8 &&
-				Nz.back()%ratio[2]==0 && Nz.back()>8 )
-		{
-			Nx.push_back( Nx.back()/ratio[0] );
-			Ny.push_back( Ny.back()/ratio[1] );
-			Nz.push_back( Nz.back()/ratio[2] );
-			// clone the domain and create coarse version based on Nx,Ny,Nz
-			//multidomain_db.push_back();
-		}
-		int N_levels = Nx.size();
-
-		// Initialize the domain
-		std::vector<std::shared_ptr<Domain>> Dm(N_levels);
-		for (int i=0; i<N_levels; i++) {
-			// This line is no good -- will create identical Domain structures instead of
-			// Need a way to define a coarse structure for the coarse domain (see above)
-			Dm[i].reset( new Domain(domain_db, comm) );
-			int N = (Nx[i]+2)*(Ny[i]+2)*(Nz[i]+2);
-			for (int n=0; n<N; n++){
-				Dm[i]->id[n] = 1;
-			}
-			Dm[i]->CommInit();
-		}
-
-		// array containing a distance mask
-		Array<float> MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2);
-
-		// Create the level data
-		std::vector<Array<char>>  ID(N_levels);
-		std::vector<Array<float>> LOCVOL(N_levels);
-		std::vector<Array<float>> Dist(N_levels);
-		std::vector<Array<float>> MultiScaleSmooth(N_levels);
-		std::vector<Array<float>> Mean(N_levels);
-		std::vector<Array<float>> NonLocalMean(N_levels);
-		std::vector<std::shared_ptr<fillHalo<double>>> fillDouble(N_levels);
-		std::vector<std::shared_ptr<fillHalo<float>>>  fillFloat(N_levels);
-		std::vector<std::shared_ptr<fillHalo<char>>>   fillChar(N_levels);
-		for (int i=0; i<N_levels; i++) {
-			ID[i] = Array<char>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
-			LOCVOL[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
-			Dist[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
-			MultiScaleSmooth[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
-			Mean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
-			NonLocalMean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
-			ID[i].fill(0);
-			LOCVOL[i].fill(0);
-			Dist[i].fill(0);
-			MultiScaleSmooth[i].fill(0);
-			Mean[i].fill(0);
-			NonLocalMean[i].fill(0);
-			fillDouble[i].reset(new fillHalo<double>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
-			fillFloat[i].reset(new fillHalo<float>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
-			fillChar[i].reset(new fillHalo<char>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
-		}
-
-		// Read the subvolume of interest on each processor
-		PROFILE_START("ReadVolume");
-		int fid = netcdf::open(InputFile,netcdf::READ);
-		std::string varname("VOLUME");
-		netcdf::VariableType type = netcdf::getVarType( fid, varname );
-		std::vector<size_t> dim = netcdf::getVarDim( fid, varname );
-		if ( rank == 0 ) {
-			printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str());
-			printf("   dims =  %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2]));
-		}
-		{
-			RankInfoStruct info( rank, nprocx, nprocy, nprocz );
-			int x = info.ix*nx;
-			int y = info.jy*ny;
-			int z = info.kz*nz;
-			// Read the local data
-			Array<short> VOLUME = netcdf::getVar<short>( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} );
-			// Copy the data and fill the halos
-			LOCVOL[0].fill(0);
-			fillFloat[0]->copy( VOLUME, LOCVOL[0] );
-			fillFloat[0]->fill( LOCVOL[0] );
-		}
-		netcdf::close( fid );
-		MPI_Barrier(comm);
-		PROFILE_STOP("ReadVolume");
-		if (rank==0) printf("Read complete\n");
-
-
-		// Filter the original data
-		filter_src( *Dm[0], LOCVOL[0] );
-
-		// Set up the mask to be distance to cylinder (crop outside cylinder)
-		//	float CylRad=900;
-		for (int k=0;k<Nz[0]+2;k++) {
-			for (int j=0;j<Ny[0]+2;j++) {
-				for (int i=0;i<Nx[0]+2;i++) {
-					int iproc = Dm[0]->iproc();
-					int jproc = Dm[0]->jproc();
-					int kproc = Dm[0]->kproc();
-
-					int x=iproc*Nx[0]+i-1;
-					int y=jproc*Ny[0]+j-1;
-					int z=kproc*Nz[0]+k-1;
-
-					//int cx = 0.5*nprocx*Nx[0];
-					//int cy = 0.5*nprocy*Ny[0];
-					//int cz = 0.5*nprocz*Nz[0];
-
-					// distance from the center line 
-					MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) );
-
-				}
-			}
-		}
-
-		// Compute the means for the high/low regions
-		// (should use automated mixture model to approximate histograms)
-		//float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) );
-		float THRESHOLD=0.5*float(target+background);
-		float mean_plus=0;
-		float mean_minus=0;
-		int count_plus=0;
-		int count_minus=0;
-		for (int k=1;k<Nz[0]+1;k++) {
-			for (int j=1;j<Ny[0]+1;j++) {
-				for (int i=1;i<Nx[0]+1;i++) {
-					if (MASK(i,j,k) > 0.f ){
-						auto tmp = LOCVOL[0](i,j,k);
-						/*						if ((tmp-background)*(tmp-target) > 0){
-							// direction to background / target is the same
-							if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background
-							else                             			 tmp=target;     // tmp closer to target
-						}
-						 */
-						if ( tmp > THRESHOLD ) {
-							mean_plus += tmp;
-							count_plus++;
-						} else if ( tmp < -THRESHOLD ) {
-							mean_minus += tmp;
-							count_minus++;
-						}
-					}
-				}
-			}
-		}
-		mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus );
-		mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus );
-		if (rank==0) printf("	Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus);
-
-
-		MPI_Barrier(comm);
-		// Scale the source data to +-1.0
-		for (size_t i=0; i<LOCVOL[0].length(); i++) {
-			if (MASK(i) < 0.f){
-				LOCVOL[0](i) = 1.0;
-			}
-			else if ( LOCVOL[0](i) >= 0 ) {
-				LOCVOL[0](i) /= mean_plus;
-				LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f );
-			} else {
-				LOCVOL[0](i) /= -mean_minus;
-				LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f );
-			}
-		}
-
-		// Fill the source data for the coarse meshes
-		PROFILE_START("CoarsenMesh");
-		for (int i=1; i<N_levels; i++) {
-			Array<float> filter(ratio[0],ratio[1],ratio[2]);
-			filter.fill(1.0f/filter.length());
-			Array<float> tmp(Nx[i-1],Ny[i-1],Nz[i-1]);
-			fillFloat[i-1]->copy( LOCVOL[i-1], tmp );
-			Array<float> coarse = tmp.coarsen( filter );
-			fillFloat[i]->copy( coarse, LOCVOL[i] );
-			fillFloat[i]->fill( LOCVOL[i] );
-			if (rank==0){
-				printf("Coarsen level %i \n",i);
-				printf("   Nx=%i, Ny=%i, Nz=%i \n",int(tmp.size(0)),int(tmp.size(1)),int(tmp.size(2))  );
-				printf("   filter_x=%i, filter_y=%i, filter_z=%i \n",int(filter.size(0)),int(filter.size(1)),int(filter.size(2))  );
-				printf("   ratio= %i,%i,%i \n",int(ratio[0]),int(ratio[1]),int(ratio[2])  );
-			}
-			MPI_Barrier(comm);
-		}
-		PROFILE_STOP("CoarsenMesh");
-
-		// Initialize the coarse level
-		PROFILE_START("Solve full mesh");
-		if (rank==0)
-			printf("Initialize full mesh\n");
-		solve( LOCVOL[0], Mean[0], ID[0], Dist[0], MultiScaleSmooth[0],
-				NonLocalMean[0], *fillFloat[0], *Dm[0], nprocx, 
-				rough_cutoff, lamda, nlm_sigsq, nlm_depth);
-		PROFILE_STOP("Solve full mesh");
-		MPI_Barrier(comm);
-
-		/*	// Initialize the coarse level
-	PROFILE_START("Solve coarse mesh");
-	if (rank==0)
-		printf("Initialize coarse mesh\n");
-	solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(),
-			NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx );
-	PROFILE_STOP("Solve coarse mesh");
-	MPI_Barrier(comm);
-
-	// Refine the solution
-	PROFILE_START("Refine distance");
-	if (rank==0)
-		printf("Refine mesh\n");
-	for (int i=int(Nx.size())-2; i>=0; i--) {
-		if (rank==0)
-			printf("   Refining to level %i\n",int(i));
-		refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i],
-				NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i );
-	}
-	PROFILE_STOP("Refine distance");
-	MPI_Barrier(comm);    
-
-	// Perform a final filter
-	PROFILE_START("Filtering final domains");
-	if (rank==0)
-		printf("Filtering final domains\n");
-	Array<float> filter_Mean, filter_Dist1, filter_Dist2;
-	filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 );
-	PROFILE_STOP("Filtering final domains");
-		 */
-		//removeDisconnected( ID[0], *Dm[0] );
-
-		/*
-        // Write the distance function to a netcdf file
-    const char* netcdf_filename = "Distance.nc";
+    // Initialize MPI
+    int rank, nprocs;
+    MPI_Init(&argc,&argv);
+    MPI_Comm comm = MPI_COMM_WORLD;
+    MPI_Comm_rank(comm,&rank);
+    MPI_Comm_size(comm,&nprocs);
     {
-        RankInfoStruct info( rank, nprocx, nprocy, nprocz );
-        std::vector<int> dim = { Nx[0]*nprocx, Ny[0]*nprocy, Nz[0]*nprocz };
-        int fid = netcdf::open( netcdf_filename, netcdf::CREATE, MPI_COMM_WORLD );
-        auto dims =  netcdf::defDim( fid, {"X", "Y", "Z"}, dim );
-        Array<float> data(Nx[0],Ny[0],Nz[0]);
-        fillFloat[0]->copy( Dist[0], data );
-        netcdf::write( fid, "Distance", dims, data, info );
+        Utilities::setErrorHandlers();
+        PROFILE_START("Main");
+
+        //std::vector<std::string> filenames;
+        if ( argc<2 ) {
+            if ( rank == 0 ){
+                printf("At least one filename must be specified\n");
+            }
+            return 1;
+        }
+        std::string filename = std::string(argv[1]);
+        if ( rank == 0 ){
+            printf("Input data file: %s\n",filename.c_str());
+        }
+
+        auto db = std::make_shared<Database>( filename );
+        auto domain_db = db->getDatabase( "Domain" );
+        auto uct_db = db->getDatabase( "uCT" );
+        auto analysis_db = db->getDatabase( "Analysis" );
+
+        // Read domain values
+        auto L = domain_db->getVector<double>( "L" );
+        auto size = domain_db->getVector<int>( "n" );
+        auto nproc = domain_db->getVector<int>( "nproc" );
+        int BoundaryCondition = domain_db->getScalar<int>( "BC" );
+        int nx = size[0];
+        int ny = size[1];
+        int nz = size[2];
+        double Lx = L[0];
+        double Ly = L[1];
+        double Lz = L[2];
+        int nprocx = nproc[0];
+        int nprocy = nproc[1];
+        int nprocz = nproc[2];
+
+        auto InputFile    = uct_db->getScalar<std::string>( "InputFile" );
+        auto target       = uct_db->getScalar<int>("target");
+        auto background   = uct_db->getScalar<int>("background");
+        auto rough_cutoff = uct_db->getScalar<float>( "rough_cutoff" );    
+        auto lamda        = uct_db->getScalar<float>( "lamda" );    
+        auto nlm_sigsq    = uct_db->getScalar<float>( "nlm_sigsq" );    
+        auto nlm_depth    = uct_db->getScalar<int>( "nlm_depth" );    
+        auto center       = uct_db->getVector<int>( "center" );
+        auto CylRad       = uct_db->getScalar<float>( "cylinder_radius" );
+        auto maxLevels    = uct_db->getScalar<int>( "max_levels" );
+        std::vector<int> offset( 3, 0 );
+        if ( uct_db->keyExists( "offset" ) )
+            offset = uct_db->getVector<int>( "offset" );
+
+        // Check that the number of processors >= the number of ranks
+        if ( rank==0 ) {
+            printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
+            printf("Number of MPI ranks used: %i \n", nprocs);
+            printf("Full domain size: %i x %i x %i  \n",nx*nprocx,ny*nprocy,nz*nprocz);
+        }
+        if ( nprocs < nprocx*nprocy*nprocz ){
+            ERROR("Insufficient number of processors");
+        }
+
+        // Determine the maximum number of levels for the desired coarsen ratio
+        int ratio[3] = {2,2,2};
+        //std::vector<size_t> ratio = {4,4,4};
+        // need to set up databases for each level of the mesh
+        std::vector<std::shared_ptr<Database>> multidomain_db(1,domain_db);
+        std::vector<int> Nx(1,nx), Ny(1,ny), Nz(1,nz);
+        while ( Nx.back()%ratio[0]==0 && Nx.back()>8 &&
+                Ny.back()%ratio[1]==0 && Ny.back()>8 &&
+                Nz.back()%ratio[2]==0 && Nz.back()>8 &&
+                (int) Nx.size() < maxLevels )
+        {
+            Nx.push_back( Nx.back()/ratio[0] );
+            Ny.push_back( Ny.back()/ratio[1] );
+            Nz.push_back( Nz.back()/ratio[2] );
+            // clone the domain and create coarse version based on Nx,Ny,Nz
+            auto db2 = domain_db->cloneDatabase();
+            db2->putVector<int>( "n", { Nx.back(), Ny.back(), Nz.back() } );
+            multidomain_db.push_back(db2);
+        }
+        int N_levels = Nx.size();
+
+        // Initialize the domain
+        std::vector<std::shared_ptr<Domain>> Dm(N_levels);
+        for (int i=0; i<N_levels; i++) {
+            // This line is no good -- will create identical Domain structures instead of
+            // Need a way to define a coarse structure for the coarse domain (see above)
+            Dm[i].reset( new Domain(multidomain_db[i], comm) );
+            int N = (Nx[i]+2)*(Ny[i]+2)*(Nz[i]+2);
+            for (int n=0; n<N; n++){
+                Dm[i]->id[n] = 1;
+            }
+            Dm[i]->CommInit();
+        }
+
+        // array containing a distance mask
+        Array<float> MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2);
+        MASK.fill(0);
+
+        // Create the level data
+        std::vector<Array<char>>  ID(N_levels);
+        std::vector<Array<float>> LOCVOL(N_levels);
+        std::vector<Array<float>> Dist(N_levels);
+        std::vector<Array<float>> MultiScaleSmooth(N_levels);
+        std::vector<Array<float>> Mean(N_levels);
+        std::vector<Array<float>> NonLocalMean(N_levels);
+        std::vector<std::shared_ptr<fillHalo<double>>> fillDouble(N_levels);
+        std::vector<std::shared_ptr<fillHalo<float>>>  fillFloat(N_levels);
+        std::vector<std::shared_ptr<fillHalo<char>>>   fillChar(N_levels);
+        for (int i=0; i<N_levels; i++) {
+            ID[i] = Array<char>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
+            LOCVOL[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
+            Dist[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
+            MultiScaleSmooth[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
+            Mean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
+            NonLocalMean[i] = Array<float>(Nx[i]+2,Ny[i]+2,Nz[i]+2);
+            ID[i].fill(0);
+            LOCVOL[i].fill(0);
+            Dist[i].fill(0);
+            MultiScaleSmooth[i].fill(0);
+            Mean[i].fill(0);
+            NonLocalMean[i].fill(0);
+            fillDouble[i].reset(new fillHalo<double>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
+            fillFloat[i].reset(new fillHalo<float>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
+            fillChar[i].reset(new fillHalo<char>(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) );
+        }
+
+        // Read the subvolume of interest on each processor
+        PROFILE_START("ReadVolume");
+        int fid = netcdf::open(InputFile,netcdf::READ);
+        std::string varname("VOLUME");
+        auto type = netcdf::getVarType( fid, varname );
+        auto dim = netcdf::getVarDim( fid, varname );
+        if ( rank == 0 ) {
+            printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str());
+            printf("   dims =  %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2]));
+        }
+        {
+            // Read the local data
+            int x = Dm[0]->iproc()*nx + offset[0];
+            int y = Dm[0]->jproc()*ny + offset[1];
+            int z = Dm[0]->kproc()*nz + offset[2];
+            Array<short> VOLUME = netcdf::getVar<short>( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} );
+            // Copy the data and fill the halos
+            LOCVOL[0].fill(0);
+            fillFloat[0]->copy( VOLUME, LOCVOL[0] );
+            fillFloat[0]->fill( LOCVOL[0] );
+        }
         netcdf::close( fid );
+        MPI_Barrier(comm);
+        PROFILE_STOP("ReadVolume");
+        if (rank==0) printf("Read complete\n");
+
+
+        // Filter the original data
+        filter_src( *Dm[0], LOCVOL[0] );
+
+        // Set up the mask to be distance to cylinder (crop outside cylinder)
+        for (int k=0;k<Nz[0]+2;k++) {
+            for (int j=0;j<Ny[0]+2;j++) {
+                for (int i=0;i<Nx[0]+2;i++) {
+                    int x=Dm[0]->iproc()*Nx[0]+i-1;
+                    int y=Dm[0]->jproc()*Ny[0]+j-1;
+                    int z=Dm[0]->kproc()*Nz[0]+k-1;
+                    int cx = center[0] - offset[0];
+                    int cy = center[1] - offset[1];
+                    int cz = center[2] - offset[2];
+                    // distance from the center line 
+                    MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) );
+                }
+            }
+        }
+
+        // Compute the means for the high/low regions
+        // (should use automated mixture model to approximate histograms)
+        //float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) );
+        double THRESHOLD=0.5*(target+background);
+        double mean_plus=0;
+        double mean_minus=0;
+        int count_plus=0;
+        int count_minus=0;
+        for (int k=1;k<Nz[0]+1;k++) {
+            for (int j=1;j<Ny[0]+1;j++) {
+                for (int i=1;i<Nx[0]+1;i++) {
+                    if (MASK(i,j,k) > 0.f ){
+                        auto tmp = LOCVOL[0](i,j,k);
+                        /*                        if ((tmp-background)*(tmp-target) > 0){
+                            // direction to background / target is the same
+                            if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background
+                            else                                          tmp=target;     // tmp closer to target
+                        }
+                         */
+                        if ( tmp > THRESHOLD ) {
+                            mean_plus += tmp;
+                            count_plus++;
+                        } else if ( tmp < -THRESHOLD ) {
+                            mean_minus += tmp;
+                            count_minus++;
+                        }
+                    }
+                }
+            }
+        }
+        ASSERT( count_plus > 0 && count_minus > 0 );
+        mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus );
+        mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus );
+        if (rank==0) printf("    Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus);
+
+
+        MPI_Barrier(comm);
+        // Scale the source data to +-1.0
+        for (size_t i=0; i<LOCVOL[0].length(); i++) {
+            if (MASK(i) < 0.f){
+                LOCVOL[0](i) = 1.0;
+            }
+            else if ( LOCVOL[0](i) >= 0 ) {
+                LOCVOL[0](i) /= mean_plus;
+                LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f );
+            } else {
+                LOCVOL[0](i) /= -mean_minus;
+                LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f );
+            }
+        }
+
+
+        // Fill the source data for the coarse meshes
+        PROFILE_START("CoarsenMesh");
+        for (int i=1; i<N_levels; i++) {
+            Array<float> filter(ratio[0],ratio[1],ratio[2]);
+            filter.fill(1.0f/filter.length());
+            Array<float> tmp(Nx[i-1],Ny[i-1],Nz[i-1]);
+            fillFloat[i-1]->copy( LOCVOL[i-1], tmp );
+            Array<float> coarse = tmp.coarsen( filter );
+            fillFloat[i]->copy( coarse, LOCVOL[i] );
+            fillFloat[i]->fill( LOCVOL[i] );
+            if (rank==0){
+                printf("Coarsen level %i \n",i);
+                printf("   Nx=%i, Ny=%i, Nz=%i \n",int(tmp.size(0)),int(tmp.size(1)),int(tmp.size(2))  );
+                printf("   filter_x=%i, filter_y=%i, filter_z=%i \n",int(filter.size(0)),int(filter.size(1)),int(filter.size(2))  );
+                printf("   ratio= %i,%i,%i \n",int(ratio[0]),int(ratio[1]),int(ratio[2])  );
+            }
+            MPI_Barrier(comm);
+        }
+        PROFILE_STOP("CoarsenMesh");
+
+        // Initialize the coarse level
+        PROFILE_START("Solve coarse mesh");
+        if (rank==0)
+            printf("Initialize full mesh\n");
+        solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(),
+                NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx, 
+                rough_cutoff, lamda, nlm_sigsq, nlm_depth);
+        PROFILE_STOP("Solve coarse mesh");
+        MPI_Barrier(comm);
+
+        // Refine the solution
+        PROFILE_START("Refine distance");
+        if (rank==0)
+            printf("Refine mesh\n");
+        for (int i=N_levels-2; i>=0; i--) {
+            if (rank==0)
+                printf("   Refining to level %i\n",i);
+            refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i],
+                    NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i, 
+                rough_cutoff, lamda, nlm_sigsq, nlm_depth);
+        }
+        PROFILE_STOP("Refine distance");
+        MPI_Barrier(comm);    
+
+        // Perform a final filter
+        PROFILE_START("Filtering final domains");
+        if (rank==0)
+            printf("Filtering final domains\n");
+        Array<float> filter_Mean, filter_Dist1, filter_Dist2;
+        filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 );
+        PROFILE_STOP("Filtering final domains");
+        //removeDisconnected( ID[0], *Dm[0] );
+        
+        // Write the distance function to a netcdf file
+        /* const char* netcdf_filename = "Distance.nc";
+        {
+            RankInfoStruct info( rank, nprocx, nprocy, nprocz );
+            std::vector<int> dim = { Nx[0]*nprocx, Ny[0]*nprocy, Nz[0]*nprocz };
+            int fid = netcdf::open( netcdf_filename, netcdf::CREATE, MPI_COMM_WORLD );
+            auto dims =  netcdf::defDim( fid, {"X", "Y", "Z"}, dim );
+            Array<float> data(Nx[0],Ny[0],Nz[0]);
+            fillFloat[0]->copy( Dist[0], data );
+            netcdf::write( fid, "Distance", dims, data, info );
+            netcdf::close( fid );
+        } */
+
+        // Write the results
+        if (rank==0) printf("Setting up visualization structure \n");
+        std::vector<IO::MeshDataStruct> meshData(N_levels);
+        for (size_t i=0; i<Nx.size(); i++) {
+            // Mesh
+            meshData[i].meshName = "image_" + std::to_string( i );
+            meshData[i].mesh = std::make_shared<IO::DomainMesh>(Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],Lx,Ly,Lz);
+            // Source data
+            auto OrigData = std::make_shared<IO::Variable>();
+            OrigData->name = "Source_Data_" + std::to_string( i );
+            OrigData->type = IO::VariableType::VolumeVariable;
+            OrigData->dim = 1;
+            OrigData->data.resize(Nx[i],Ny[i],Nz[i]);
+            meshData[i].vars.push_back(OrigData);
+            fillDouble[i]->copy( LOCVOL[i], OrigData->data );
+            // Non-Local Mean
+            auto NonLocMean = std::make_shared<IO::Variable>();
+            NonLocMean->name = "NonLocal_Mean_" + std::to_string( i );
+            NonLocMean->type = IO::VariableType::VolumeVariable;
+            NonLocMean->dim = 1;
+            NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]);
+            meshData[i].vars.push_back(NonLocMean);
+            fillDouble[i]->copy( NonLocalMean[i], NonLocMean->data );
+            // Segmented Data
+            auto SegData = std::make_shared<IO::Variable>();
+            SegData->name = "Segmented_Data_" + std::to_string( i );
+            SegData->type = IO::VariableType::VolumeVariable;
+            SegData->dim = 1;
+            SegData->data.resize(Nx[i],Ny[i],Nz[i]);
+            meshData[i].vars.push_back(SegData);
+            fillDouble[i]->copy( ID[i], SegData->data );
+            // Signed Distance
+            auto DistData = std::make_shared<IO::Variable>();
+            DistData->name = "Signed_Distance_" + std::to_string( i );
+            DistData->type = IO::VariableType::VolumeVariable;
+            DistData->dim = 1;
+            DistData->data.resize(Nx[i],Ny[i],Nz[i]);
+            meshData[i].vars.push_back(DistData);
+            fillDouble[i]->copy( Dist[i], DistData->data );
+            // Smoothed Data
+            auto SmoothData = std::make_shared<IO::Variable>();
+            SmoothData->name = "Smoothed_Data_" + std::to_string( i );
+            SmoothData->type = IO::VariableType::VolumeVariable;
+            SmoothData->dim = 1;
+            SmoothData->data.resize(Nx[i],Ny[i],Nz[i]);
+            meshData[i].vars.push_back(SmoothData);
+            fillDouble[i]->copy( MultiScaleSmooth[i], SmoothData->data );
+
+        }
+        #if 0
+            std::shared_ptr<IO::Variable> filter_Mean_var( new IO::Variable() );
+            filter_Mean_var->name = "Mean";
+            filter_Mean_var->type = IO::VariableType::VolumeVariable;
+            filter_Mean_var->dim = 1;
+            filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]);
+            meshData[0].vars.push_back(filter_Mean_var);
+            fillDouble[0]->copy( filter_Mean, filter_Mean_var->data );
+            std::shared_ptr<IO::Variable> filter_Dist1_var( new IO::Variable() );
+            filter_Dist1_var->name = "Dist1";
+            filter_Dist1_var->type = IO::VariableType::VolumeVariable;
+            filter_Dist1_var->dim = 1;
+            filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]);
+            meshData[0].vars.push_back(filter_Dist1_var);
+            fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data );
+            std::shared_ptr<IO::Variable> filter_Dist2_var( new IO::Variable() );
+            filter_Dist2_var->name = "Dist2";
+            filter_Dist2_var->type = IO::VariableType::VolumeVariable;
+            filter_Dist2_var->dim = 1;
+            filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]);
+            meshData[0].vars.push_back(filter_Dist2_var);
+            fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data );
+        #endif
+        MPI_Barrier(comm);
+        if (rank==0) printf("Writing output \n");
+        // Write visulization data
+        IO::writeData( 0, meshData, comm );
+        if (rank==0) printf("Finished. \n");
+    
+        // Compute the Minkowski functionals
+        MPI_Barrier(comm);
+        auto Averages = std::make_shared<Minkowski>(Dm[0]);
+        
+        Array <char> phase_label(Nx[0]+2,Ny[0]+2,Nz[0]+2);
+        Array <double> phase_distance(Nx[0]+2,Ny[0]+2,Nz[0]+2);
+        // Analyze the wetting fluid
+        for (int k=1;k<Nz[0]+1;k++) {
+            for (int j=1;j<Ny[0]+1;j++) {
+                for (int i=1;i<Nx[0]+1;i++) {
+                    int n = k*Nx[0]*Ny[0]+j*Nx[0]+i;
+                    if (!(Dm[0]->id[n] > 0)){
+                        // Solid phase
+                        phase_label(i,j,k) = 0;
+                    }
+                    else if (Dist[0](i,j,k) < 0.0){
+                        // wetting phase
+                        phase_label(i,j,k) = 1;
+                    }
+                    else {
+                        // non-wetting phase
+                        phase_label(i,j,k) = 0;
+                    }
+                    phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0;
+                }
+            }
+        }    
+        CalcDist(phase_distance,phase_label,*Dm[0]);
+        Averages->ComputeScalar(phase_distance,0.f);
+        Averages->PrintAll();
     }
-		 */
-
-		{
-			// Write the results
-			if (rank==0) printf("Setting up visualization structure \n");
-			//	std::vector<IO::MeshDataStruct> meshData(N_levels);
-			std::vector<IO::MeshDataStruct> meshData(1);
-			//    for (size_t i=0; i<Nx.size(); i++) {
-			// Mesh
-			meshData[0].meshName = "image";
-			meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm[0]->rank_info,Nx[0],Ny[0],Nz[0],Lx,Ly,Lz) );
-			// Source data
-			std::shared_ptr<IO::Variable> OrigData( new IO::Variable() );
-			OrigData->name = "Source Data";
-			OrigData->type = IO::VariableType::VolumeVariable;
-			OrigData->dim = 1;
-			OrigData->data.resize(Nx[0],Ny[0],Nz[0]);
-			meshData[0].vars.push_back(OrigData);
-			fillDouble[0]->copy( LOCVOL[0], OrigData->data );
-			// Non-Local Mean
-			std::shared_ptr<IO::Variable> NonLocMean( new IO::Variable() );
-			NonLocMean->name = "Non-Local Mean";
-			NonLocMean->type = IO::VariableType::VolumeVariable;
-			NonLocMean->dim = 1;
-			NonLocMean->data.resize(Nx[0],Ny[0],Nz[0]);
-			meshData[0].vars.push_back(NonLocMean);
-			fillDouble[0]->copy( NonLocalMean[0], NonLocMean->data );
-			std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
-			SegData->name = "Segmented Data";
-			SegData->type = IO::VariableType::VolumeVariable;
-			SegData->dim = 1;
-			SegData->data.resize(Nx[0],Ny[0],Nz[0]);
-			meshData[0].vars.push_back(SegData);
-			fillDouble[0]->copy( ID[0], SegData->data );
-			// Signed Distance
-			std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
-			DistData->name = "Signed Distance";
-			DistData->type = IO::VariableType::VolumeVariable;
-			DistData->dim = 1;
-			DistData->data.resize(Nx[0],Ny[0],Nz[0]);
-			meshData[0].vars.push_back(DistData);
-			fillDouble[0]->copy( Dist[0], DistData->data );
-			// Smoothed Data
-			std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
-			SmoothData->name = "Smoothed Data";
-			SmoothData->type = IO::VariableType::VolumeVariable;
-			SmoothData->dim = 1;
-			SmoothData->data.resize(Nx[0],Ny[0],Nz[0]);
-			meshData[0].vars.push_back(SmoothData);
-			fillDouble[0]->copy( MultiScaleSmooth[0], SmoothData->data );
-
-			/*// Segmented Data
-	       std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
-	    SegData->name = "Segmented Data";
-	    SegData->type = IO::VariableType::VolumeVariable;
-	    SegData->dim = 1;
-	    SegData->data.resize(Nx[i],Ny[i],Nz[i]);
-	    meshData[i].vars.push_back(SegData);
-        fillDouble[i]->copy( ID[i], SegData->data );
-        // Signed Distance
-        std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
-	    DistData->name = "Signed Distance";
-	    DistData->type = IO::VariableType::VolumeVariable;
-	    DistData->dim = 1;
-	    DistData->data.resize(Nx[i],Ny[i],Nz[i]);
-	    meshData[i].vars.push_back(DistData);
-        fillDouble[i]->copy( Dist[i], DistData->data );
-        // Smoothed Data
-        std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
-	    SmoothData->name = "Smoothed Data";
-	    SmoothData->type = IO::VariableType::VolumeVariable;
-	    SmoothData->dim = 1;
-	    SmoothData->data.resize(Nx[i],Ny[i],Nz[i]);
-	    meshData[i].vars.push_back(SmoothData);
-        fillDouble[i]->copy( MultiScaleSmooth[i], SmoothData->data );
-
-    }
-      #if 0
-        std::shared_ptr<IO::Variable> filter_Mean_var( new IO::Variable() );
-	    filter_Mean_var->name = "Mean";
-	    filter_Mean_var->type = IO::VariableType::VolumeVariable;
-	    filter_Mean_var->dim = 1;
-	    filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]);
-	    meshData[0].vars.push_back(filter_Mean_var);
-        fillDouble[0]->copy( filter_Mean, filter_Mean_var->data );
-        std::shared_ptr<IO::Variable> filter_Dist1_var( new IO::Variable() );
-	    filter_Dist1_var->name = "Dist1";
-	    filter_Dist1_var->type = IO::VariableType::VolumeVariable;
-	    filter_Dist1_var->dim = 1;
-	    filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]);
-	    meshData[0].vars.push_back(filter_Dist1_var);
-        fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data );
-        std::shared_ptr<IO::Variable> filter_Dist2_var( new IO::Variable() );
-	    filter_Dist2_var->name = "Dist2";
-	    filter_Dist2_var->type = IO::VariableType::VolumeVariable;
-	    filter_Dist2_var->dim = 1;
-	    filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]);
-	    meshData[0].vars.push_back(filter_Dist2_var);
-        fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data );
-    #endif
-			 */
-			MPI_Barrier(comm);
-			if (rank==0) printf("Writing output \n");
-			// Write visulization data
-			IO::writeData( 0, meshData, comm );
-			if (rank==0) printf("Finished. \n");
-		}
-	
-		// Compute the Minkowski functionals
-		MPI_Barrier(comm);
-		std::shared_ptr<Minkowski> Averages(new Minkowski(Dm[0]));
-		
-		Array <char> phase_label(Nx[0],Ny[0],Nz[0]);
-		Array <double> phase_distance(Nx[0],Ny[0],Nz[0]);
-		// Analyze the wetting fluid
-		for (int k=1;k<Nz[0]+1;k++) {
-			for (int j=1;j<Ny[0]+1;j++) {
-				for (int i=1;i<Nx[0]+1;i++) {
-					int n = k*Nx[0]*Ny[0]+j*Nx[0]+i;
-					if (!(Dm[0]->id[n] > 0)){
-						// Solid phase
-						phase_label(i,j,k) = 0;
-					}
-					else if (Dist[0](i,j,k) < 0.0){
-						// wetting phase
-						phase_label(i,j,k) = 1;
-					}
-					else {
-						// non-wetting phase
-						phase_label(i,j,k) = 0;
-					}
-					phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0;
-				}
-			}
-		}	
-		CalcDist(phase_distance,phase_label,*Dm[0]);
-		Averages->ComputeScalar(phase_distance,0.f);
-		Averages->PrintAll();
-	}
-	PROFILE_STOP("Main");
-	PROFILE_SAVE("lbpm_uCT_pp",true);
-	MPI_Barrier(comm);
-	MPI_Finalize();
-	return 0;
+    PROFILE_STOP("Main");
+    PROFILE_SAVE("lbpm_uCT_pp",true);
+    MPI_Barrier(comm);
+    MPI_Finalize();
+    return 0;
 }