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

IO/netcdf.cpp

@@ -204,6 +204,7 @@ std::vector<size_t> getAttDim( int fid, const std::string& att )
 {
     std::vector<size_t> dim(1,0);
     int err = nc_inq_attlen( fid, NC_GLOBAL, att.c_str(), dim.data() );
+    CHECK_NC_ERR( err );
     return dim;
 }
 std::vector<std::string> getVarNames( int fid )

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    = WEXITSTATUS( status );
+    int code = pclose( pipe );
+    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

@@ -261,6 +261,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();
 
@@ -304,6 +308,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 {
 
 
@@ -97,12 +109,11 @@ inline size_t findfirst( const std::vector<TYPE> &X, TYPE Y )
  *  Function to terminate the program                                        *
  ****************************************************************************/
 static bool abort_throwException = false;
-static printStackType abort_stackType = printStackType::global;
 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/SubPhase.cpp

@@ -161,15 +161,14 @@ void SubPhase::Basic(){
 	// If external boundary conditions are set, do not average over the inlet
 	kmin=1; kmax=Nz-1;
 	imin=jmin=1;
-	// If inlet/outlet layers exist use these as default
+	/*// If inlet/outlet layers exist use these as default
 	if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
 	if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
 	if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z; 
 	if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z; 
-	
+	*/
 	nb.reset(); wb.reset();
 
-	double nA,nB;
 	double count_w = 0.0;
 	double count_n = 0.0;
 	
@@ -281,7 +280,7 @@ void SubPhase::Basic(){
 			dir_y = 0.0;
 			dir_z = 1.0;
 		}
-		if (Dm->BoundaryCondition > 0 ){
+		if (Dm->BoundaryCondition == 1 || Dm->BoundaryCondition == 2 || Dm->BoundaryCondition == 3 || Dm->BoundaryCondition == 4 ){
 			// compute the pressure drop
 			double pressure_drop = (Pressure(Nx*Ny + Nx + 1) - 1.0) / 3.0;
 			double length = ((Nz-2)*Dm->nprocz());
@@ -297,8 +296,8 @@ void SubPhase::Basic(){
 		double saturation=gwb.V/(gwb.V + gnb.V);
 		double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume;
 		double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume;
-		double total_flow_rate = water_flow_rate + not_water_flow_rate;
-		double fractional_flow= water_flow_rate / total_flow_rate;
+		//double total_flow_rate = water_flow_rate + not_water_flow_rate;
+		//double fractional_flow = water_flow_rate / total_flow_rate;
 
 		double h = Dm->voxel_length;		
 		double krn = h*h*nu_n*not_water_flow_rate / force_mag ;
@@ -377,16 +376,17 @@ void SubPhase::Full(){
 
 	// If external boundary conditions are set, do not average over the inlet
 	kmin=1; kmax=Nz-1;
-	if (Dm->BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
-	if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
-
+	/*if (Dm->BoundaryCondition > 0 && Dm->BoundaryCondition != 5 && Dm->kproc() == 0) kmin=4;
+	if (Dm->BoundaryCondition > 0 && Dm->BoundaryCondition != 5 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
+	*/
 	imin=jmin=1;
-	// If inlet layers exist use these as default
+	/*// If inlet layers exist use these as default
+	 * NOTE -- excluding inlet / outlet will screw up topological averages!!!
 	if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
 	if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
 	if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z; 
 	if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z; 
-		
+	*/
 	nd.reset();	nc.reset(); wd.reset();	wc.reset();	iwn.reset();	iwnc.reset();
 
  	Dm->CommunicateMeshHalo(Phi);
@@ -697,7 +697,8 @@ void SubPhase::Full(){
 }
 
 
-void SubPhase::AggregateLabels(char *FILENAME){
+void SubPhase::AggregateLabels( const std::string& filename )
+{
 	
 	int nx = Dm->Nx;
 	int ny = Dm->Ny;
@@ -721,7 +722,7 @@ void SubPhase::AggregateLabels(char *FILENAME){
 	}
 	MPI_Barrier(Dm->Comm);
 
-	Dm->AggregateLabels(FILENAME);
+	Dm->AggregateLabels( filename );
 
 }
 

analysis/SubPhase.h

@@ -101,7 +101,7 @@ public:
 	void Basic();
 	void Full();
 	void Write(int time);
-    void AggregateLabels(char *FILENAME);
+    void AggregateLabels( const std::string& filename );
 
 private:
 	FILE *TIMELOG;

analysis/TwoPhase.cpp

@@ -234,6 +234,7 @@ TwoPhase::~TwoPhase()
 
 void TwoPhase::ColorToSignedDistance(double Beta, DoubleArray &ColorData, DoubleArray &DistData)
 {
+    NULL_USE( Beta );
   /*double factor,temp,value;
 		factor=0.5/Beta;
 	// Initialize to -1,1 (segmentation)
@@ -657,8 +658,8 @@ void TwoPhase::ComputeLocal()
 
 void TwoPhase::AssignComponentLabels()
 {
-	int LabelNWP=1;
-	int LabelWP=2;
+	//int LabelNWP=1;
+	//int LabelWP=2;
 	// NOTE: labeling the wetting phase components is tricky! One sandstone media had over 800,000 components
 	// NumberComponents_WP = ComputeGlobalPhaseComponent(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->rank_info,PhaseID,LabelWP,Label_WP);
 	// treat all wetting phase is connected
@@ -1202,6 +1203,8 @@ void TwoPhase::Reduce()
 
 void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT)
 {
+    NULL_USE( viscosity );
+    NULL_USE( IFT );
 	awn_global *= D;
 	ans_global *= D;
 	ans_global *= D;

analysis/dcel.cpp

@@ -352,6 +352,8 @@ double DECL::EdgeAngle(int edge)
 
 void Isosurface(DoubleArray &A, const double &v)
 {
+    NULL_USE( v );
+
 	Point P,Q;
 	Point PlaceHolder;
 	Point C0,C1,C2,C3,C4,C5,C6,C7;
@@ -562,7 +564,7 @@ void Isosurface(DoubleArray &A, const double &v)
 					if (P.z == 1.0 && Q.z == 1.0) HalfEdge[idx_edge][3] = -6;  // ghost twin for z=1 face
 				}
 				// Find all the angles
-				for (int idx=0; idx<EdgeCount; idx++){
+				/*for (int idx=0; idx<EdgeCount; idx++){
 					int V1=HalfEdge[idx][0];
 					int V2=HalfEdge[idx][1];
 					int T1= HalfEdge[idx_edge][2];
@@ -570,7 +572,7 @@ void Isosurface(DoubleArray &A, const double &v)
 					if (twin == -1){
 						
 					}
-				}
+				}*/
 
 				// Map vertices to global coordinates
 				for (int idx=0;idx<VertexCount;idx++) {

analysis/morphology.cpp

@@ -34,9 +34,6 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 	int nx = Dm->Nx;
 	int ny = Dm->Ny;
 	int nz = Dm->Nz;
-	int iproc = Dm->iproc();
-	int jproc = Dm->jproc();
-	int kproc = Dm->kproc();
 	int nprocx = Dm->nprocx();
 	int nprocy = Dm->nprocy();
 	int nprocz = Dm->nprocz();
@@ -122,7 +119,6 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 	int sendtag,recvtag;
 	sendtag = recvtag = 7;
 
-	int x,y,z;
 	int ii,jj,kk;
 	int Nx = nx;
 	int Ny = ny;
@@ -135,7 +131,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 
 	// Increase the critical radius until the target saturation is met
 	double deltaR=0.05; // amount to change the radius in voxel units
-	double Rcrit_old;
+	double Rcrit_old=0.0;
 
 	double GlobalNumber = 1.f;
 	int imin,jmin,kmin,imax,jmax,kmax;
@@ -336,9 +332,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 	int nx = Dm->Nx;
 	int ny = Dm->Ny;
 	int nz = Dm->Nz;
-	int iproc = Dm->iproc();
-	int jproc = Dm->jproc();
-	int kproc = Dm->kproc();
 	int nprocx = Dm->nprocx();
 	int nprocy = Dm->nprocy();
 	int nprocz = Dm->nprocz();
@@ -427,7 +420,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 	int sendtag,recvtag;
 	sendtag = recvtag = 7;
 
-	int x,y,z;
 	int ii,jj,kk;
 	int Nx = nx;
 	int Ny = ny;
@@ -693,17 +685,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain
 	return final_void_fraction;
 }
 
-double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetGrowth){
-	
+double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetGrowth, double WallFactor)
+{
 	int Nx = Dm->Nx;
 	int Ny = Dm->Ny;
 	int Nz = Dm->Nz;
-	int iproc = Dm->iproc();
-	int jproc = Dm->jproc();
-	int kproc = Dm->kproc();
-	int nprocx = Dm->nprocx();
-	int nprocy = Dm->nprocy();
-	int nprocz = Dm->nprocz();
 	int rank = Dm->rank();
 		
 	double count=0.0;
@@ -736,8 +722,7 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
 			for (int j=1; j<Ny-1; j++){
 				for (int i=1; i<Nx-1; i++){
 					double walldist=BoundaryDist(i,j,k);
-					double wallweight = 1.0 / (1+exp(-5.f*(walldist-1.f))); 
-					//wallweight = 1.0;
+					double wallweight = WallFactor/ (1+exp(-5.f*(walldist-1.f))); 
 					if (fabs(wallweight*morph_delta) > MAX_DISPLACEMENT) MAX_DISPLACEMENT= fabs(wallweight*morph_delta);
 					
 					if (Dist(i,j,k) - wallweight*morph_delta < 0.0){
@@ -783,7 +768,7 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
 		for (int j=1; j<Ny-1; j++){
 			for (int i=1; i<Nx-1; i++){
 				double walldist=BoundaryDist(i,j,k);
-				double wallweight = 1.0 / (1+exp(-5.f*(walldist-1.f))); 
+				double wallweight = WallFactor / (1+exp(-5.f*(walldist-1.f))); 
 				//wallweight = 1.0;
 				Dist(i,j,k) -= wallweight*morph_delta;
 				if (Dist(i,j,k) < 0.0)	count+=1.0;

analysis/morphology.h

@@ -5,4 +5,4 @@
 
 double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction, signed char ErodeLabel, signed char ReplaceLabel);
 double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
-double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetVol);
+double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetVol, double WallFactor);

analysis/runAnalysis.cpp

@@ -782,6 +782,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 );
@@ -920,12 +922,12 @@ void runAnalysis::run(int timestep, std::shared_ptr<Database> input_db, TwoPhase
     // Spawn a thread to write the restart file
     //    if ( matches(type,AnalysisType::CreateRestart) ) {
     if (timestep%d_restart_interval==0){
-
+    	auto Restart_db = input_db->cloneDatabase();
+	//	Restart_db->putScalar<bool>( "Restart", true );
     	if (d_rank==0) {
-    		input_db->putScalar<bool>( "Restart", true );
-    		std::ofstream OutStream("Restart.db");
-    		input_db->print(OutStream, "");
-    		OutStream.close();
+	  //	std::ofstream OutStream("Restart.db");
+	  //	Restart_db->print(OutStream, "");
+	  //	OutStream.close();
     	}
     	// Write the restart file (using a seperate thread)
         auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cDen,cfq,d_Np);
@@ -952,8 +954,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" );
@@ -969,7 +969,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();
@@ -998,7 +998,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) ) {
@@ -1025,21 +1024,21 @@ void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPha
     	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));
-
+    	// clone the input database to avoid modifying shared data
+    	auto Restart_db = input_db->cloneDatabase();
+    	auto tmp_color_db =  Restart_db->getDatabase( "Color" );
+    	tmp_color_db->putScalar<int>("timestep",timestep);    		
+    	tmp_color_db->putScalar<bool>( "Restart", true );
+    	Restart_db->putDatabase("Color", tmp_color_db);
     	if (d_rank==0) {
-    		color_db->putScalar<int>("timestep",timestep);    		
-    		color_db->putScalar<bool>( "Restart", true );
-    		input_db->putDatabase("Color", color_db);
     		std::ofstream OutStream("Restart.db");
-    		input_db->print(OutStream, "");
+    		Restart_db->print(OutStream, "");
     		OutStream.close();
-  
     	}
     	// Write the restart file (using a seperate thread)
     	auto work1 = new WriteRestartWorkItem(d_restartFile.c_str(),cDen,cfq,d_Np);
     	work1->add_dependency(d_wait_restart);
     	d_wait_restart = d_tpool.add_work(work1);
-
     }
     
     if (timestep%d_visualization_interval==0){
@@ -1051,12 +1050,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 );
@@ -1083,7 +1081,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");
 }

analysis/uCT.cpp

@@ -172,6 +172,7 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
     //	int depth = 5;
     //	float sigsq=0.1;
 	int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq);
+    NULL_USE( nlm_count );
 	fillFloat.fill(NonLocalMean);
 }
 
@@ -216,6 +217,7 @@ void refine( const Array<float>& Dist_coarse,
 //	int depth = 3;
 //	float sigsq = 0.1;
 	int nlm_count = NLM3D( MultiScaleSmooth, Mean, Dist, NonLocalMean, depth, sigsq);
+    NULL_USE( nlm_count );
 	fillFloat.fill(NonLocalMean);
     segment( NonLocalMean, ID, 0.001 );
     for (size_t i=0; i<ID.length(); i++) {

cmake/libraries.cmake

@@ -305,4 +305,9 @@ MACRO ( CONFIGURE_LBPM )
         SET( CMAKE_BUILD_WITH_INSTALL_RPATH TRUE )
         SET( CMAKE_INSTALL_RPATH ${CMAKE_INSTALL_RPATH} "${TIMER_DIRECTORY}" "${LBPM_INSTALL_DIR}/lib" )
     ENDIF()
+    # Suppress some common warnings
+    IF ( USING_GCC )
+        SET( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-reorder -Wno-unused-parameter")
+        SET( CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} --compiler-options -Wno-reorder,-Wno-unused-parameter")
+    ENDIF()
 ENDMACRO ()

cmake/macros.cmake

@@ -848,7 +848,7 @@ FUNCTION( ADD_${PROJ}_TEST EXEFILE ${ARGN} )
     ADD_PROJ_PROVISIONAL_TEST( ${EXEFILE} )
     CREATE_TEST_NAME( ${EXEFILE} ${ARGN} )
     IF ( USE_MPI_FOR_SERIAL_TESTS )
-        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} "${MPIEXEC_NUMPROC_FLAG}" 1 $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
+        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} ${MPIFLAGS} "${MPIEXEC_NUMPROC_FLAG}" 1 $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
         SET_PROPERTY( TEST ${TESTNAME} APPEND PROPERTY ENVIRONMENT OMPI_MCA_hwloc_base_binding_policy=none )
     ELSE()
         ADD_TEST( NAME ${TESTNAME} COMMAND $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
@@ -877,7 +877,7 @@ FUNCTION( ADD_${PROJ}_WEEKLY_TEST EXEFILE PROCS ${ARGN} )
     ELSEIF( ${PROCS} STREQUAL "1" )
         CREATE_TEST_NAME( "${EXEFILE}_WEEKLY" ${ARGN} )
         IF ( USE_MPI_FOR_SERIAL_TESTS )
-            ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} "${MPIEXEC_NUMPROC_FLAG}" 1 $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
+            ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} ${MPIFLAGS} "${MPIEXEC_NUMPROC_FLAG}" 1 $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
             SET_PROPERTY( TEST ${TESTNAME} APPEND PROPERTY ENVIRONMENT OMPI_MCA_hwloc_base_binding_policy=none )
         ELSE()
             ADD_TEST( NAME ${TESTNAME} COMMAND $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
@@ -909,7 +909,7 @@ FUNCTION( ADD_${PROJ}_TEST_PARALLEL EXEFILE PROCS ${ARGN} )
     ELSEIF ( ${PROCS} GREATER ${TEST_MAX_PROCS} )
         MESSAGE("Disabling test ${TESTNAME} (exceeds maximum number of processors ${TEST_MAX_PROCS})")
     ELSE()
-        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} "${MPIEXEC_NUMPROC_FLAG}" ${PROCS} $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
+        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} ${MPIFLAGS} "${MPIEXEC_NUMPROC_FLAG}" ${PROCS} $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
         SET_PROPERTY( TEST ${TESTNAME} APPEND PROPERTY ENVIRONMENT OMPI_MCA_hwloc_base_binding_policy=none )
         SET_TESTS_PROPERTIES( ${TESTNAME} PROPERTIES FAIL_REGULAR_EXPRESSION "${TEST_FAIL_REGULAR_EXPRESSION}" PROCESSORS ${PROCS} )
         ADD_RESOURCE_LOCK( ${TESTNAME} ${EXEFILE} ${ARGN} )
@@ -930,7 +930,7 @@ MACRO( ADD_${PROJ}_TEST_THREAD_MPI EXEFILE PROCS THREADS ${ARGN} )
         SET_TESTS_PROPERTIES ( ${TESTNAME} PROPERTIES FAIL_REGULAR_EXPRESSION "${TEST_FAIL_REGULAR_EXPRESSION}" PROCESSORS ${TOT_PROCS} )
         ADD_RESOURCE_LOCK( ${TESTNAME} ${EXEFILE} ${ARGN} )
     ELSEIF ( USE_MPI OR USE_EXT_MPI )
-        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} "${MPIEXEC_NUMPROC_FLAG}" ${PROCS} $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
+        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} ${MPIFLAGS} "${MPIEXEC_NUMPROC_FLAG}" ${PROCS} $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
         SET_PROPERTY( TEST ${TESTNAME} APPEND PROPERTY ENVIRONMENT OMPI_MCA_hwloc_base_binding_policy=none )
         SET_TESTS_PROPERTIES ( ${TESTNAME} PROPERTIES FAIL_REGULAR_EXPRESSION "${TEST_FAIL_REGULAR_EXPRESSION}" PROCESSORS ${TOT_PROCS} )
         ADD_RESOURCE_LOCK( ${TESTNAME} ${EXEFILE} ${ARGN} )
@@ -966,7 +966,7 @@ FUNCTION( ADD_${PROJ}_EXAMPLE EXEFILE PROCS ${ARGN} )
         ADD_TEST( NAME ${TESTNAME} COMMAND $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
     ELSEIF ( USE_EXT_MPI AND NOT (${PROCS} GREATER ${TEST_MAX_PROCS}) )
         CREATE_TEST_NAME( "example--${EXEFILE}_${PROCS}procs" ${ARGN} )
-        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} "${MPIEXEC_NUMPROC_FLAG}" ${PROCS} $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
+        ADD_TEST( NAME ${TESTNAME} COMMAND ${MPIEXEC} ${MPIFLAGS} "${MPIEXEC_NUMPROC_FLAG}" ${PROCS} $<TARGET_FILE:${LAST_TEST}> ${ARGN} )
         SET_PROPERTY( TEST ${TESTNAME} APPEND PROPERTY ENVIRONMENT OMPI_MCA_hwloc_base_binding_policy=none )
     ENDIF()
     SET_TESTS_PROPERTIES( ${TESTNAME} PROPERTIES FAIL_REGULAR_EXPRESSION "${TEST_FAIL_REGULAR_EXPRESSION}" PROCESSORS ${PROCS} )

common/Communication.hpp

@@ -74,9 +74,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 ( size_t j=0; j<dst_rank.ny; j++ ) {
-                for ( size_t k=0; k<dst_rank.nz; k++ ) {
+        for ( int i=0; i<dst_rank.nx; i++ ) {
+            for ( int j=0; j<dst_rank.ny; j++ ) {
+                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 );
@@ -95,9 +95,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 ( size_t j=0; j<src_rank.ny; j++ ) {
-            for ( size_t k=0; k<src_rank.nz; k++ ) {
+    for ( int i=0; i<src_rank.nx; i++ ) {
+        for ( int j=0; j<src_rank.ny; j++ ) {
+            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/Domain.cpp

@@ -73,6 +73,9 @@ Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
 	recvData_xY(NULL), recvData_yZ(NULL), recvData_Xz(NULL), recvData_XY(NULL), recvData_YZ(NULL), recvData_XZ(NULL),
 	id(NULL)
 {	
+    NULL_USE( rnk );
+    NULL_USE( npy );
+    NULL_USE( npz );
 	// set up the neighbor ranks
     int myrank;
     MPI_Comm_rank( Comm, &myrank );
@@ -256,7 +259,7 @@ void Domain::initialize( std::shared_ptr<Database> db )
 	INSIST(nprocs == nproc[0]*nproc[1]*nproc[2],"Fatal error in processor count!");
 }
 
-void Domain::Decomp(std::string Filename)
+void Domain::Decomp( const std::string& Filename )
 {
 	//.......................................................................
 	// Reading the domain information file
@@ -266,9 +269,9 @@ void Domain::Decomp(std::string Filename)
 	int nprocs, nprocx, nprocy, nprocz, nx, ny, nz;
 	int64_t global_Nx,global_Ny,global_Nz;
 	int64_t i,j,k,n;
-	int BC=0;
 	int64_t xStart,yStart,zStart;
 	int checkerSize;
+	bool USE_CHECKER = false;
 	//int inlet_layers_x, inlet_layers_y, inlet_layers_z;
 	//int outlet_layers_x, outlet_layers_y, outlet_layers_z;
 	xStart=yStart=zStart=0;
@@ -308,6 +311,7 @@ void Domain::Decomp(std::string Filename)
 	}
 	if (database->keyExists( "checkerSize" )){
 		checkerSize = database->getScalar<int>( "checkerSize" );
+		USE_CHECKER = true;
 	}
 	else {
 		checkerSize = SIZE[0];
@@ -346,7 +350,7 @@ void Domain::Decomp(std::string Filename)
 	if (RANK==0){
 		printf("Input media: %s\n",Filename.c_str());
 		printf("Relabeling %lu values\n",ReadValues.size());
-		for (int idx=0; idx<ReadValues.size(); idx++){
+		for (size_t idx=0; idx<ReadValues.size(); idx++){
 			int oldvalue=ReadValues[idx];
 			int newvalue=WriteValues[idx];
 			printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
@@ -389,7 +393,7 @@ void Domain::Decomp(std::string Filename)
 					n = k*global_Nx*global_Ny+j*global_Nx+i;
 					//char locval = loc_id[n];
 					char locval = SegData[n];
-					for (int idx=0; idx<ReadValues.size(); idx++){
+					for (size_t idx=0; idx<ReadValues.size(); idx++){
 						signed char oldvalue=ReadValues[idx];
 						signed char newvalue=WriteValues[idx];
 						if (locval == oldvalue){
@@ -401,14 +405,12 @@ void Domain::Decomp(std::string Filename)
 				}
 			}
 		}
-		if (RANK==0){
-			for (int idx=0; idx<ReadValues.size(); idx++){
+		for (size_t idx=0; idx<ReadValues.size(); idx++){
 			long int label=ReadValues[idx];
 			long int count=LabelCount[idx];
-				printf("Label=%d, Count=%d \n",label,count);
+			printf("Label=%ld, Count=%ld \n",label,count);
 		}
-		}
-
+		if (USE_CHECKER) {
 			if (inlet_layers_x > 0){
 				// use checkerboard pattern
 				printf("Checkerboard pattern at x inlet for %i layers \n",inlet_layers_x);
@@ -448,7 +450,7 @@ void Domain::Decomp(std::string Filename)
 			}
 
 			if (inlet_layers_z > 0){
-			printf("Checkerboard pattern at z inlet for %i layers \n",inlet_layers_z);
+				printf("Checkerboard pattern at z inlet for %i layers, saturated with phase label=%i \n",inlet_layers_z,inlet_layers_phase);
 				// use checkerboard pattern
 				for (int k = zStart; k < zStart+inlet_layers_z; k++){
 					for (int j = 0; j<global_Ny; j++){
@@ -490,7 +492,7 @@ void Domain::Decomp(std::string Filename)
 				printf("Checkerboard pattern at y outlet for %i layers \n",outlet_layers_y);
 				// use checkerboard pattern
 				for (int k = 0; k<global_Nz; k++){
-				for (int j = yStart + ny*nprocy - outlet_layers_y; i < yStart + ny*nprocy; j++){
+					for (int j = yStart + ny*nprocy - outlet_layers_y; j < yStart + ny*nprocy; j++){
 						for (int i = 0; i<global_Nx; i++){
 							if ( (i/checkerSize + k/checkerSize)%2 == 0){
 								// void checkers
@@ -507,7 +509,7 @@ void Domain::Decomp(std::string Filename)
 			}
 
 			if (outlet_layers_z > 0){
-			printf("Checkerboard pattern at z outlet for %i layers \n",outlet_layers_z);
+				printf("Checkerboard pattern at z outlet for %i layers, saturated with phase label=%i \n",outlet_layers_z,outlet_layers_phase);
 				// use checkerboard pattern
 				for (int k = zStart + nz*nprocz - outlet_layers_z; k < zStart + nz*nprocz; k++){
 					for (int j = 0; j<global_Ny; j++){
@@ -526,6 +528,37 @@ void Domain::Decomp(std::string Filename)
 				}
 			}
 		}
+		else {
+			if (inlet_layers_z > 0){
+				printf("Mixed reflection pattern at z inlet for %i layers, saturated with phase label=%i \n",inlet_layers_z,inlet_layers_phase);
+				for (int k = zStart; k < zStart+inlet_layers_z; k++){
+					for (int j = 0; j<global_Ny; j++){
+						for (int i = 0; i<global_Nx; i++){
+							signed char local_id = SegData[k*global_Nx*global_Ny+j*global_Nx+i];
+							signed char reflection_id = SegData[(zStart + nz*nprocz - 1)*global_Nx*global_Ny+j*global_Nx+i];
+							if ( local_id < 1 && reflection_id > 0){
+								SegData[k*global_Nx*global_Ny+j*global_Nx+i] = reflection_id;
+							}
+						}
+					}
+				}
+			}
+			if (outlet_layers_z > 0){
+				printf("Mixed reflection pattern at z outlet for %i layers, saturated with phase label=%i \n",outlet_layers_z,outlet_layers_phase);
+				for (int k = zStart + nz*nprocz - outlet_layers_z; k < zStart + nz*nprocz; k++){
+					for (int j = 0; j<global_Ny; j++){
+						for (int i = 0; i<global_Nx; i++){
+							signed char local_id = SegData[k*global_Nx*global_Ny+j*global_Nx+i];
+							signed char reflection_id = SegData[zStart*global_Nx*global_Ny+j*global_Nx+i];
+							if ( local_id < 1 && reflection_id > 0){
+								SegData[k*global_Nx*global_Ny+j*global_Nx+i] = reflection_id;
+							}
+						}
+					}
+				}
+			}
+		}
+	}
 	
 	// Get the rank info
 	int64_t N = (nx+2)*(ny+2)*(nz+2);
@@ -599,10 +632,65 @@ void Domain::Decomp(std::string Filename)
 		//printf("Ready to recieve data %i at process %i \n", N,rank);
 		MPI_Recv(id,N,MPI_CHAR,0,15,Comm,MPI_STATUS_IGNORE);
 	}
+	//Comm.barrier();
 	MPI_Barrier(Comm);
+	// Compute the porosity
+	double sum;
+	double sum_local=0.0;
+	double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+	if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx*(Ny-2)*nprocy*((Nz-2)*nprocz-6));
+	//.........................................................
+	// If external boundary conditions are applied remove solid
+	if (BoundaryCondition >  0 && BoundaryCondition !=5 && kproc() == 0){
+    	if (inlet_layers_z < 4){
+            inlet_layers_z=4;
+            if(RANK==0){
+                printf("NOTE:Non-periodic BC is applied, but the number of Z-inlet layers is not specified (or is smaller than 3 voxels) \n     the number of Z-inlet layer is reset to %i voxels, saturated with phase label=%i \n",inlet_layers_z-1,inlet_layers_phase);
+            } 
+        }	
+		for (int k=0; k<inlet_layers_z; k++){
+			for (int j=0;j<Ny;j++){
+				for (int i=0;i<Nx;i++){
+					int n = k*Nx*Ny+j*Nx+i;
+					id[n] = inlet_layers_phase;
+				}                    
+			}
+ 		}
+ 	}
+    if (BoundaryCondition >  0 && BoundaryCondition !=5 && kproc() == nprocz-1){
+    	if (outlet_layers_z < 4){
+            outlet_layers_z=4;
+            if(RANK==nprocs-1){
+                printf("NOTE:Non-periodic BC is applied, but the number of Z-outlet layers is not specified (or is smaller than 3 voxels) \n     the number of Z-outlet layer is reset to %i voxels, saturated with phase label=%i \n",outlet_layers_z-1,outlet_layers_phase);
+            } 
+        }	
+ 		for (int k=Nz-outlet_layers_z; 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;
+ 					id[n] = outlet_layers_phase;
+ 				}                    
+ 			}
+ 		}
+ 	}
+    for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
+        for (int j=1;j<Ny-1;j++){
+            for (int i=1;i<Nx-1;i++){
+                int n = k*Nx*Ny+j*Nx+i;
+                if (id[n] > 0){
+                    sum_local+=1.0;
+                }
+            }
+        }
+    }
+    MPI_Allreduce(&sum_local,&sum,1,MPI_DOUBLE,MPI_SUM,Comm);
+    //sum = Comm.sumReduce(sum_local);
+    porosity = sum*iVol_global;
+    if (rank()==0) printf("Media porosity = %f \n",porosity);
+ 	//.........................................................
 }
 
-void Domain::AggregateLabels(char *FILENAME){
+void Domain::AggregateLabels( const std::string& filename ){
 	
 	int nx = Nx;
 	int ny = Ny;
@@ -679,8 +767,7 @@ void Domain::AggregateLabels(char *FILENAME){
 			}
 		}
 		// write the output
-		FILE *OUTFILE;
-		OUTFILE = fopen(FILENAME,"wb");
+		FILE *OUTFILE = fopen(filename.c_str(),"wb");
 		fwrite(FullID,1,full_size,OUTFILE);
 		fclose(OUTFILE);
 	}
@@ -1020,10 +1107,10 @@ void Domain::ReadIDs(){
 	double sum;
 	double sum_local=0.0;
 	double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
-	if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
+	if (BoundaryCondition > 0 && BoundaryCondition !=5) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
 	//.........................................................
 	// If external boundary conditions are applied remove solid
-	if (BoundaryCondition >  0  && kproc() == 0){
+	if (BoundaryCondition >  0 && BoundaryCondition !=5 && kproc() == 0){
     	if (inlet_layers_z < 4)	inlet_layers_z=4;
 		for (int k=0; k<inlet_layers_z; k++){
 			for (int j=0;j<Ny;j++){
@@ -1034,7 +1121,7 @@ void Domain::ReadIDs(){
 			}
  		}
  	}
-    if (BoundaryCondition >  0  && kproc() == nprocz()-1){
+    if (BoundaryCondition >  0 && BoundaryCondition !=5 && kproc() == nprocz()-1){
     	if (outlet_layers_z < 4)	outlet_layers_z=4;
  		for (int k=Nz-outlet_layers_z; k<Nz; k++){
  			for (int j=0;j<Ny;j++){
@@ -1160,19 +1247,18 @@ void Domain::CommunicateMeshHalo(DoubleArray &Mesh)
 }
 
 // Ideally stuff below here should be moved somewhere else -- doesn't really belong here
-void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np)
+void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, size_t Np)
 {
-    int q,n;
     double value;
     ofstream File(FILENAME,ios::binary);
-    for (n=0; n<Np; n++){
+    for (size_t n=0; n<Np; n++){
         // Write the two density values
         value = cDen[n];
         File.write((char*) &value, sizeof(value));
         value = cDen[Np+n];
         File.write((char*) &value, sizeof(value));
         // Write the even distributions
-        for (q=0; q<19; q++){
+        for (size_t q=0; q<19; q++){
             value = cfq[q*Np+n];
             File.write((char*) &value, sizeof(value));
         }
@@ -1181,16 +1267,15 @@ void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq
 
 }
 
-void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, int Np)
+void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, size_t Np)
 {
-    int q=0, n=0;
     double value=0;
     ifstream File(FILENAME,ios::binary);
-    for (n=0; n<Np; n++){
+    for (size_t n=0; n<Np; n++){
         File.read((char*) &value, sizeof(value));
         cPhi[n] = value;
         // Read the distributions
-        for (q=0; q<19; q++){
+        for (size_t q=0; q<19; q++){
             File.read((char*) &value, sizeof(value));
             cfq[q*Np+n] = value;
         }
@@ -1198,13 +1283,12 @@ void ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, int Np)
     File.close();
 }
 
-void ReadBinaryFile(char *FILENAME, double *Data, int N)
+void ReadBinaryFile(char *FILENAME, double *Data, size_t N)
 {
-  int n;
   double value;
   ifstream File(FILENAME,ios::binary);
   if (File.good()){
-    for (n=0; n<N; n++){
+    for (size_t n=0; n<N; n++){
       // Write the two density values                                                                                
       File.read((char*) &value, sizeof(value));
       Data[n] = value;
@@ -1212,7 +1296,7 @@ void ReadBinaryFile(char *FILENAME, double *Data, int N)
     }
   }
   else {
-    for (n=0; n<N; n++) Data[n] = 1.2e-34;
+    for (size_t n=0; n<N; n++) Data[n] = 1.2e-34;
   }
   File.close();
 }

common/Domain.h

@@ -192,11 +192,11 @@ public: // Public variables (need to create accessors instead)
     signed char *id;
 
     void ReadIDs();
-    void Decomp(std::string Filename);
+    void Decomp( const std::string& filename );
     void CommunicateMeshHalo(DoubleArray &Mesh);
     void CommInit(); 
     int PoreCount();
-    void AggregateLabels(char *FILENAME);
+    void AggregateLabels( const std::string& filename );
 
 private:
 
@@ -259,10 +259,10 @@ private:
 
 };
 
-void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, int Np);
+void WriteCheckpoint(const char *FILENAME, const double *cDen, const double *cfq, size_t Np);
 
-void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, int Np);
+void ReadCheckpoint(char *FILENAME, double *cDen, double *cfq, size_t Np);
 
-void ReadBinaryFile(char *FILENAME, double *Data, int N);
+void ReadBinaryFile(char *FILENAME, double *Data, size_t N);
 
 #endif

common/MPI_Helpers.h

@@ -203,6 +203,12 @@ inline int sumReduce( MPI_Comm comm, int x )
 	MPI_Allreduce(&x,&y,1,MPI_INT,MPI_SUM,comm);
     return y;
 }
+inline long long sumReduce( MPI_Comm comm, long long x )
+{
+    long long y = 0;
+	MPI_Allreduce(&x,&y,1,MPI_LONG_LONG,MPI_SUM,comm);
+    return y;
+}
 inline bool sumReduce( MPI_Comm comm, bool x )
 {
     int y = sumReduce( comm, x?1:0 );

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/ScaLBL.cpp

@@ -1026,6 +1026,52 @@ void ScaLBL_Communicator::RecvD3Q19AA(double *dist){
 	ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Y,3*recvCount_Y,recvCount_Y,recvbuf_Y,dist,N);
 	ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Y,4*recvCount_Y,recvCount_Y,recvbuf_Y,dist,N);
 	//...................................................................................
+	//...Pack the xy edge (8)................................
+	ScaLBL_D3Q19_Unpack(8,dvcRecvDist_xy,0,recvCount_xy,recvbuf_xy,dist,N);
+	//...Pack the Xy edge (9)................................
+	ScaLBL_D3Q19_Unpack(9,dvcRecvDist_Xy,0,recvCount_Xy,recvbuf_Xy,dist,N);
+	//...Pack the xY edge (10)................................
+	ScaLBL_D3Q19_Unpack(10,dvcRecvDist_xY,0,recvCount_xY,recvbuf_xY,dist,N);
+	//...Pack the XY edge (7)................................
+	ScaLBL_D3Q19_Unpack(7,dvcRecvDist_XY,0,recvCount_XY,recvbuf_XY,dist,N);
+
+	if (BoundaryCondition > 0){
+		if (kproc != 0){
+			//...Packing for z face(6,12,13,16,17)................................
+			ScaLBL_D3Q19_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,dist,N);
+			ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,dist,N);
+			ScaLBL_D3Q19_Unpack(13,dvcRecvDist_z,2*recvCount_z,recvCount_z,recvbuf_z,dist,N);
+			ScaLBL_D3Q19_Unpack(16,dvcRecvDist_z,3*recvCount_z,recvCount_z,recvbuf_z,dist,N);
+			ScaLBL_D3Q19_Unpack(17,dvcRecvDist_z,4*recvCount_z,recvCount_z,recvbuf_z,dist,N);
+			//...Pack the xz edge (12)................................
+			ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,recvbuf_xz,dist,N);
+			//...Pack the Xz edge (13)................................
+			ScaLBL_D3Q19_Unpack(13,dvcRecvDist_Xz,0,recvCount_Xz,recvbuf_Xz,dist,N);
+			//...Pack the yz edge (16)................................
+			ScaLBL_D3Q19_Unpack(16,dvcRecvDist_yz,0,recvCount_yz,recvbuf_yz,dist,N);
+			//...Pack the Yz edge (17)................................
+			ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,recvbuf_Yz,dist,N);
+			//..................................................................................
+		}
+		if (kproc != nprocz-1){
+			//...Packing for Z face(5,11,14,15,18)................................
+			ScaLBL_D3Q19_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,dist,N);
+			ScaLBL_D3Q19_Unpack(11,dvcRecvDist_Z,recvCount_Z,recvCount_Z,recvbuf_Z,dist,N);
+			ScaLBL_D3Q19_Unpack(14,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N);
+			ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N);
+			ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N);
+			//...Pack the xZ edge (14)................................
+			ScaLBL_D3Q19_Unpack(14,dvcRecvDist_xZ,0,recvCount_xZ,recvbuf_xZ,dist,N);
+			//...Pack the XZ edge (11)................................
+			ScaLBL_D3Q19_Unpack(11,dvcRecvDist_XZ,0,recvCount_XZ,recvbuf_XZ,dist,N);
+			//...Pack the yZ edge (18)................................
+			ScaLBL_D3Q19_Unpack(18,dvcRecvDist_yZ,0,recvCount_yZ,recvbuf_yZ,dist,N);
+			//...Pack the YZ edge (15)................................
+			ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,recvbuf_YZ,dist,N);
+			//..................................................................................
+		}
+	}
+	else {
 		//...Packing for z face(6,12,13,16,17)................................
 		ScaLBL_D3Q19_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,dist,N);
 		ScaLBL_D3Q19_Unpack(12,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,dist,N);
@@ -1039,14 +1085,6 @@ void ScaLBL_Communicator::RecvD3Q19AA(double *dist){
 		ScaLBL_D3Q19_Unpack(15,dvcRecvDist_Z,3*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N);
 		ScaLBL_D3Q19_Unpack(18,dvcRecvDist_Z,4*recvCount_Z,recvCount_Z,recvbuf_Z,dist,N);
 		//..................................................................................
-	//...Pack the xy edge (8)................................
-	ScaLBL_D3Q19_Unpack(8,dvcRecvDist_xy,0,recvCount_xy,recvbuf_xy,dist,N);
-	//...Pack the Xy edge (9)................................
-	ScaLBL_D3Q19_Unpack(9,dvcRecvDist_Xy,0,recvCount_Xy,recvbuf_Xy,dist,N);
-	//...Pack the xY edge (10)................................
-	ScaLBL_D3Q19_Unpack(10,dvcRecvDist_xY,0,recvCount_xY,recvbuf_xY,dist,N);
-	//...Pack the XY edge (7)................................
-	ScaLBL_D3Q19_Unpack(7,dvcRecvDist_XY,0,recvCount_XY,recvbuf_XY,dist,N);
 		//...Pack the xz edge (12)................................
 		ScaLBL_D3Q19_Unpack(12,dvcRecvDist_xz,0,recvCount_xz,recvbuf_xz,dist,N);
 		//...Pack the xZ edge (14)................................
@@ -1063,6 +1101,8 @@ void ScaLBL_Communicator::RecvD3Q19AA(double *dist){
 		ScaLBL_D3Q19_Unpack(17,dvcRecvDist_Yz,0,recvCount_Yz,recvbuf_Yz,dist,N);
 		//...Pack the YZ edge (15)................................
 		ScaLBL_D3Q19_Unpack(15,dvcRecvDist_YZ,0,recvCount_YZ,recvbuf_YZ,dist,N);
+	}
+
 	//...................................................................................
 	Lock=false; // unlock the communicator after communications complete
 	//...................................................................................
@@ -1241,17 +1281,17 @@ void ScaLBL_Communicator::BiRecvD3Q7AA(double *Aq, double *Bq){
 	ScaLBL_D3Q7_Unpack(3,dvcRecvDist_Y,recvCount_Y,recvCount_Y,recvbuf_Y,Bq,N);
 	//...................................................................................
 
-	if (BoundaryCondition > 0 && kproc == 0){
-		// don't unpack little z
+	if (BoundaryCondition > 0){
+		if (kproc != 0){
+			//...Packing for z face(6,12,13,16,17)................................
+			ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,Aq,N);
+			ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,Bq,N);
+		}
+		if (kproc != nprocz-1){
 			//...Packing for Z face(5,11,14,15,18)................................
 			ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,Aq,N);
 			ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,recvCount_Z,recvCount_Z,recvbuf_Z,Bq,N);
 		}
-	else if (BoundaryCondition > 0 && kproc == nprocz-1){
-		// don't unpack big z
-		//...Packing for z face(6,12,13,16,17)................................
-		ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,Aq,N);
-		ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,Bq,N);
 	}
 	else {
 		//...Packing for z face(6,12,13,16,17)................................
@@ -1261,7 +1301,17 @@ void ScaLBL_Communicator::BiRecvD3Q7AA(double *Aq, double *Bq){
 		ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,Aq,N);
 		ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,recvCount_Z,recvCount_Z,recvbuf_Z,Bq,N);
 	}
-	
+/*	if (BoundaryCondition == 5){
+		if (kproc == 0){
+			ScaLBL_D3Q7_Reflection_BC_z(dvcSendList_z, Aq, sendCount_z, N);
+			ScaLBL_D3Q7_Reflection_BC_z(dvcSendList_z, Bq, sendCount_z, N);
+		}
+		if (kproc == nprocz-1){
+			ScaLBL_D3Q7_Reflection_BC_Z(dvcSendList_Z, Aq, sendCount_Z, N);
+			ScaLBL_D3Q7_Reflection_BC_Z(dvcSendList_Z, Bq, sendCount_Z, N);
+		}
+	}
+	*/
 	//...................................................................................
 	Lock=false; // unlock the communicator after communications complete
 	//...................................................................................
@@ -1358,19 +1408,19 @@ void ScaLBL_Communicator::TriRecvD3Q7AA(double *Aq, double *Bq, double *Cq){
 	ScaLBL_D3Q7_Unpack(3,dvcRecvDist_Y,2*recvCount_Y,recvCount_Y,recvbuf_Y,Cq,N);
 	//...................................................................................
 	
-	if (BoundaryCondition > 0 && kproc == 0){
-		// don't unpack little z
+	if (BoundaryCondition > 0){
+		if (kproc != 0){
+			//...Packing for z face(6,12,13,16,17)................................
+			ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,Aq,N);
+			ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,Bq,N);
+			ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,2*recvCount_z,recvCount_z,recvbuf_z,Cq,N);
+		}
+		if (kproc != nprocz-1){
 			//...Packing for Z face(5,11,14,15,18)................................
 			ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,0,recvCount_Z,recvbuf_Z,Aq,N);
 			ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,recvCount_Z,recvCount_Z,recvbuf_Z,Bq,N);
 			ScaLBL_D3Q7_Unpack(5,dvcRecvDist_Z,2*recvCount_Z,recvCount_Z,recvbuf_Z,Cq,N);
 		}
-	else if (BoundaryCondition > 0 && kproc == nprocz-1){
-		// don't unpack big z
-		//...Packing for z face(6,12,13,16,17)................................
-		ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,0,recvCount_z,recvbuf_z,Aq,N);
-		ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,recvCount_z,recvCount_z,recvbuf_z,Bq,N);
-		ScaLBL_D3Q7_Unpack(6,dvcRecvDist_z,2*recvCount_z,recvCount_z,recvbuf_z,Cq,N);
 	}
 	else {
 		//...Packing for z face(6,12,13,16,17)................................
@@ -1472,30 +1522,57 @@ void ScaLBL_Communicator::RecvHalo(double *data){
 	//...................................................................................
 	ScaLBL_Scalar_Unpack(dvcRecvList_x, recvCount_x,recvbuf_x, data, N);
 	ScaLBL_Scalar_Unpack(dvcRecvList_y, recvCount_y,recvbuf_y, data, N);
-	ScaLBL_Scalar_Unpack(dvcRecvList_z, recvCount_z,recvbuf_z, data, N);
 	ScaLBL_Scalar_Unpack(dvcRecvList_X, recvCount_X,recvbuf_X, data, N);
 	ScaLBL_Scalar_Unpack(dvcRecvList_Y, recvCount_Y,recvbuf_Y, data, N);
-	ScaLBL_Scalar_Unpack(dvcRecvList_Z, recvCount_Z,recvbuf_Z, data, N);
 	ScaLBL_Scalar_Unpack(dvcRecvList_xy, recvCount_xy,recvbuf_xy, data, N);
 	ScaLBL_Scalar_Unpack(dvcRecvList_xY, recvCount_xY,recvbuf_xY, data, N);
 	ScaLBL_Scalar_Unpack(dvcRecvList_Xy, recvCount_Xy,recvbuf_Xy, data, N);
 	ScaLBL_Scalar_Unpack(dvcRecvList_XY, recvCount_XY,recvbuf_XY, data, N);
+	
+	if (BoundaryCondition > 0){
+		if (kproc != 0){
+			//...Packing for z face(6,12,13,16,17)................................
+			ScaLBL_Scalar_Unpack(dvcRecvList_z, recvCount_z,recvbuf_z, data, N);
 			ScaLBL_Scalar_Unpack(dvcRecvList_xz, recvCount_xz,recvbuf_xz, data, N);
-	ScaLBL_Scalar_Unpack(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, data, N);
 			ScaLBL_Scalar_Unpack(dvcRecvList_Xz, recvCount_Xz,recvbuf_Xz, data, N);
-	ScaLBL_Scalar_Unpack(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, data, N);
 			ScaLBL_Scalar_Unpack(dvcRecvList_yz, recvCount_yz,recvbuf_yz, data, N);
-	ScaLBL_Scalar_Unpack(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, data, N);
 			ScaLBL_Scalar_Unpack(dvcRecvList_Yz, recvCount_Yz,recvbuf_Yz, data, N);
+		}
+		if (kproc != nprocz-1){
+			//...Packing for Z face(5,11,14,15,18)................................
+			ScaLBL_Scalar_Unpack(dvcRecvList_Z, recvCount_Z,recvbuf_Z, data, N);
+			ScaLBL_Scalar_Unpack(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, data, N);
+			ScaLBL_Scalar_Unpack(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, data, N);
+			ScaLBL_Scalar_Unpack(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, data, N);
 			ScaLBL_Scalar_Unpack(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, data, N);
+		}
+	}
+	else {
+		ScaLBL_Scalar_Unpack(dvcRecvList_z, recvCount_z,recvbuf_z, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_xz, recvCount_xz,recvbuf_xz, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_Xz, recvCount_Xz,recvbuf_Xz, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_yz, recvCount_yz,recvbuf_yz, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_Yz, recvCount_Yz,recvbuf_Yz, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_Z, recvCount_Z,recvbuf_Z, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_xZ, recvCount_xZ,recvbuf_xZ, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_XZ, recvCount_XZ,recvbuf_XZ, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_yZ, recvCount_yZ,recvbuf_yZ, data, N);
+		ScaLBL_Scalar_Unpack(dvcRecvList_YZ, recvCount_YZ,recvbuf_YZ, data, N);
+	}
 	//...................................................................................
 	Lock=false; // unlock the communicator after communications complete
 	//...................................................................................
+	if (BoundaryCondition == 5 && kproc == 0){
+		ScaLBL_CopySlice_z(data,Nx,Ny,Nz,1,0);
+	}
+	if (BoundaryCondition == 5 && kproc == nprocz-1){
+		ScaLBL_CopySlice_z(data,Nx,Ny,Nz,Nz-2,Nz-1);
+	}	
 }
 
 void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, DoubleArray &regdata){
 	// Gets data from the device and stores in regular layout
-	int i,j,k,n,idx;
+	int i,j,k,idx;
 	int Nx = map.size(0);
 	int Ny = map.size(1);
 	int Nz = map.size(2);
@@ -1524,25 +1601,32 @@ void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, Double
 	delete [] TmpDat;
 }
 
-
 void ScaLBL_Communicator::Color_BC_z(int *Map, double *Phi, double *Den, double vA, double vB){
-	double Value=(vA-vB)/(vA+vB);
 	if (kproc == 0) {
+		if (BoundaryCondition == 5){
+			//ScaLBL_CopySlice_z(Phi,Nx,Ny,Nz,1,0);
+		}
+		else {
 			// Set the phase indicator field and density on the z inlet
 			ScaLBL_Color_BC_z(dvcSendList_z, Map, Phi, Den, vA, vB, sendCount_z, N);
+		}
 		//ScaLBL_SetSlice_z(Phi,Value,Nx,Ny,Nz,0);
 	}
 }
 
 void ScaLBL_Communicator::Color_BC_Z(int *Map, double *Phi, double *Den, double vA, double vB){
-	double Value=(vA-vB)/(vA+vB);
 	if (kproc == nprocz-1){
+		if (BoundaryCondition == 5){
+			//ScaLBL_CopySlice_z(Phi,Nx,Ny,Nz,Nz-2,Nz-1);
+		}
+		else {
 		// Set the phase indicator field and density on the Z outlet
 			ScaLBL_Color_BC_Z(dvcSendList_Z, Map, Phi, Den, vA, vB, sendCount_Z, N);
-		//ScaLBL_SetSlice_z(Phi,Value,Nx,Ny,Nz,Nz-1);
 		}
 	}
 
+}
+
 void ScaLBL_Communicator::D3Q19_Pressure_BC_z(int *neighborList, double *fq, double din, int time){
     //ScaLBL_D3Q19_Pressure_BC_z(int *LIST,fq,din,Nx,Ny,Nz);
 	if (kproc == 0) {
@@ -1608,6 +1692,17 @@ double ScaLBL_Communicator::D3Q19_Flux_BC_z(int *neighborList, double *fq, doubl
 	return din;
 }
 
+void ScaLBL_Communicator::D3Q19_Reflection_BC_z(double *fq){
+	if (kproc == 0)
+		ScaLBL_D3Q19_Reflection_BC_z(dvcSendList_z, fq, sendCount_z, N);
+	
+}
+
+void ScaLBL_Communicator::D3Q19_Reflection_BC_Z(double *fq){
+	if (kproc == nprocz-1)
+		ScaLBL_D3Q19_Reflection_BC_Z(dvcSendList_Z, fq, sendCount_Z, N);
+}
+
 void ScaLBL_Communicator::PrintD3Q19(){
 	printf("Printing D3Q19 communication buffer contents \n");
 

common/ScaLBL.h

@@ -61,6 +61,7 @@ extern "C" void ScaLBL_UnpackDenD3Q7(int *list, int count, double *recvbuf, int
 
 extern "C" void ScaLBL_D3Q19_Init(double *Dist, int Np);
 
+
 extern "C" void ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np);
 
 extern "C" void ScaLBL_D3Q19_Pressure(double *dist, double *press, int Np);
@@ -70,6 +71,23 @@ extern "C" void ScaLBL_D3Q19_AAeven_BGK(double *dist, int start, int finish, int
 
 extern "C" void ScaLBL_D3Q19_AAodd_BGK(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double Fx, double Fy, double Fz);
 
+// GREYSCALE MODEL
+
+extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *Dist, int Np, double Den);
+
+extern "C" void ScaLBL_D3Q19_AAeven_Greyscale(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz,
+                                              double *Poros,double *Perm, double *Velocity,double *Pressure);
+
+extern "C" void ScaLBL_D3Q19_AAodd_Greyscale(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz, 
+                                             double *Poros,double *Perm, double *Velocity,double *Pressure);
+
+extern "C" void ScaLBL_D3Q19_AAeven_Greyscale_IMRT(double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz,
+                                              double *Poros,double *Perm, double *Velocity,double Den,double *Pressure);
+
+extern "C" void ScaLBL_D3Q19_AAodd_Greyscale_IMRT(int *neighborList, double *dist, int start, int finish, int Np, double rlx, double rlx_eff, double Fx, double Fy, double Fz, 
+                                             double *Poros,double *Perm, double *Velocity,double Den,double *Pressure);
+
+
 // MRT MODEL
 extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx,
 		double Fy, double Fz);
@@ -116,11 +134,6 @@ extern "C" void ScaLBL_D3Q19_Gradient_DFH(int *NeighborList, double *Phi, double
 
 // BOUNDARY CONDITION ROUTINES
 
-//extern "C" void ScaLBL_D3Q19_Pressure_BC_z(double *disteven, double *distodd, double din,
-//		int Nx, int Ny, int Nz);
-//extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(double *disteven, double *distodd, double dout,
-//		int Nx, int Ny, int Nz, int outlet);
-
 extern "C" void ScaLBL_D3Q19_AAodd_Pressure_BC_z(int *neighborList, int *list, double *dist, double din, int count, int Np);
 
 extern "C" void ScaLBL_D3Q19_AAodd_Pressure_BC_Z(int *neighborList, int *list, double *dist, double dout, int count, int Np);
@@ -139,8 +152,18 @@ extern "C" void ScaLBL_Color_BC_z(int *list, int *Map, double *Phi, double *Den,
 
 extern "C" void ScaLBL_Color_BC_Z(int *list, int *Map, double *Phi, double *Den, double vA, double vB, int count, int Np);
 
+extern "C" void ScaLBL_D3Q19_Reflection_BC_z(int *list, double *dist, int count, int Np);
+
+extern "C" void ScaLBL_D3Q19_Reflection_BC_Z(int *list, double *dist, int count, int Np);
+
+extern "C" void ScaLBL_D3Q7_Reflection_BC_z(int *list, double *dist, int count, int Np);
+
+extern "C" void ScaLBL_D3Q7_Reflection_BC_Z(int *list, double *dist, int count, int Np);
+
 extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, int Nz, int Slice);
 
+extern "C" void ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int Source, int Destination);
+
 class ScaLBL_Communicator{
 public:
 	//......................................................................................
@@ -149,6 +172,7 @@ public:
 	//ScaLBL_Communicator(Domain &Dm, IntArray &Map);
 	~ScaLBL_Communicator();
 	//......................................................................................
+	MPI_Comm MPI_COMM_SCALBL;		// MPI Communicator
 	unsigned long int CommunicationCount,SendCount,RecvCount;
 	int Nx,Ny,Nz,N;
 	int BoundaryCondition;
@@ -174,18 +198,8 @@ public:
 	int LastInterior();
 	
 	int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np);
-//	void MemoryOptimizedLayout(IntArray &Map, int *neighborList, char *id, int Np);
-//	void MemoryOptimizedLayoutFull(IntArray &Map, int *neighborList, char *id, int Np);
-//	void MemoryDenseLayout(IntArray &Map, int *neighborList, char *id, int Np);
-//	void MemoryDenseLayoutFull(IntArray &Map, int *neighborList, char *id, int Np);
-//	void SendD3Q19(double *f_even, double *f_odd);
-//	void RecvD3Q19(double *f_even, double *f_odd);
-//	void SendD3Q19AA(double *f_even, double *f_odd);
-//	void RecvD3Q19AA(double *f_even, double *f_odd);
 	void SendD3Q19AA(double *dist);
 	void RecvD3Q19AA(double *dist);
-//	void BiSendD3Q7(double *A_even, double *A_odd, double *B_even, double *B_odd);
-//	void BiRecvD3Q7(double *A_even, double *A_odd, double *B_even, double *B_odd);
 	void BiSendD3Q7AA(double *Aq, double *Bq);
 	void BiRecvD3Q7AA(double *Aq, double *Bq);
 	void TriSendD3Q7AA(double *Aq, double *Bq, double *Cq);
@@ -200,11 +214,10 @@ public:
 	void Color_BC_Z(int *Map, double *Phi, double *Den, double vA, double vB);
 	void D3Q19_Pressure_BC_z(int *neighborList, double *fq, double din, int time);
 	void D3Q19_Pressure_BC_Z(int *neighborList, double *fq, double dout, int time);
+	void D3Q19_Reflection_BC_z(double *fq);
+	void D3Q19_Reflection_BC_Z(double *fq);
 	double D3Q19_Flux_BC_z(int *neighborList, double *fq, double flux, int time);
 
-//	void TestSendD3Q19(double *f_even, double *f_odd);
-//	void TestRecvD3Q19(double *f_even, double *f_odd);
-
 	// Debugging and unit testing functions
 	void PrintD3Q19();
 
@@ -222,7 +235,6 @@ private:
 	// Give the object it's own MPI communicator
 	RankInfoStruct rank_info;
 	MPI_Group Group;	// Group of processors associated with this domain
-	MPI_Comm MPI_COMM_SCALBL;		// MPI Communicator for this domain
 	MPI_Request req1[18],req2[18];
 	MPI_Status stat1[18],stat2[18];
 	//......................................................................................

common/Utilities.cpp

@@ -14,12 +14,132 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 #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
+// OS specific includes / definitions
+// clang-format off
+#if defined( WIN32 ) || defined( _WIN32 ) || defined( WIN64 ) || defined( _WIN64 )
+    #define USE_WINDOWS
+#elif defined( __APPLE__ )
+    #define USE_MAC
+#elif defined( __linux ) || defined( __linux__ ) || defined( __unix ) || defined( __posix )
+    #define USE_LINUX
+#else
+    #error Unknown OS
+#endif
+// clang-format on
+
+
+// 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 ) 
@@ -69,9 +189,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

@@ -40,6 +40,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, ... );
 

cpu/Color.cpp

@@ -2825,3 +2825,11 @@ extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, doubl
 	}
 }
 
+extern "C" void ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int Source, int Dest){
+	int n; double value;
+	for (n=0; n<Nx*Ny; n++){
+		value = Phi[Source*Nx*Ny+n];
+		Phi[Dest*Nx*Ny+n] = value;
+	}
+}
+

cpu/D3Q7.cpp

@@ -87,6 +87,23 @@ extern  "C" void ScaLBL_UnpackDenD3Q7(int *list, int count, double *recvbuf, int
 	}
 }
 
+extern "C" void ScaLBL_D3Q7_Reflection_BC_z(int *list, double *dist, int count, int Np){
+  int n;
+	for (int idx=0; idx<count; idx++){
+		n = list[idx];
+		double f5 = 0.222222222222222222222222 - dist[6*Np+n];
+		dist[6*Np+n] = f5;
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_Reflection_BC_Z(int *list, double *dist, int count, int Np){
+  int n;
+  for (int idx=0; idx<count; idx++){
+		n = list[idx];
+		double f6 = 0.222222222222222222222222 - dist[5*Np+n];
+		dist[5*Np+n] = f6;
+	}
+}
 extern "C" void ScaLBL_D3Q7_Init(char *ID, double *f_even, double *f_odd, double *Den, int Nx, int Ny, int Nz)
 {
 	int n,N;

cpu/dfh.cpp

@@ -53,11 +53,10 @@ extern "C" void ScaLBL_Gradient_Unpack(double weight, double Cqx, double Cqy, do
 	}
 }
 
-extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np){
-	int idx,n;
+extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq, int start, int finish, int Np)
+{
+	for (int idx=start; idx<finish; idx++){
 	    double phi,nA,nB;
-
-	for (idx=start; idx<finish; idx++){
 		phi = Phi[idx];
 		if (phi > 0.f){
 			nA = 1.0; nB = 0.f;
@@ -90,15 +89,13 @@ extern "C" void ScaLBL_DFH_Init(double *Phi, double *Den, double *Aq, double *Bq
 // LBM based on density functional hydrodynamics
 extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
 		double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
-		double Fx, double Fy, double Fz, int start, int finish, int Np){
-
-	int ijk,nn,n;
+		double Fx, double Fy, double Fz, int start, int finish, int Np)
+{
 	double fq;
 	// conserved momemnts
 	double rho,jx,jy,jz;
 	// non-conserved moments
 	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
-	double m3,m5,m7;
 	double nA,nB; // number density
 	double a1,b1,a2,b2,nAB,delta;
 	double C,nx,ny,nz; //color gradient magnitude and direction
@@ -616,7 +613,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *
 		double *Phi, double *Gradient, double *SolidForce, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int start, int finish, int Np){
 	
-	int n,nn,ijk,nread;
+	int nread;
 	int nr1,nr2,nr3,nr4,nr5,nr6;
 	int nr7,nr8,nr9,nr10;
 	int nr11,nr12,nr13,nr14;
@@ -626,7 +623,6 @@ extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *
 	double rho,jx,jy,jz;
 	// non-conserved moments
 	double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
-	double m3,m5,m7;
 	double nA,nB; // number density
 	double a1,b1,a2,b2,nAB,delta;
 	double C,nx,ny,nz; //color gradient magnitude and direction
@@ -1212,12 +1208,12 @@ extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *
 }
 
 extern "C" void ScaLBL_D3Q7_AAodd_DFH(int *neighborList, double *Aq, double *Bq, 
-			double *Den, double *Phi, int start, int finish, int Np){
-
-	int idx,n,nread;
-	double fq,nA,nB;
+			double *Den, double *Phi, int start, int finish, int Np)
+{
 
 	for (int n=start; n<finish; n++){
+		int nread;
+		double fq,nA,nB;
 		
 		//..........Compute the number density for component A............
 		// q=0
@@ -1300,11 +1296,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_DFH(int *neighborList, double *Aq, double *Bq,
 }
 
 extern "C" void ScaLBL_D3Q7_AAeven_DFH(double *Aq, double *Bq, double *Den, double *Phi, 
-			int start, int finish, int Np){
-	int idx,n,nread;
-	double fq,nA,nB;
+			int start, int finish, int Np)
+{
 	for (int n=start; n<finish; n++){
-		
+		double fq,nA,nB;
 		// compute number density for component A
 		// q=0
 		fq = Aq[n];

example/Workflow/ComputeSaturation.py

@@ -0,0 +1,37 @@
+import sys
+import numpy as np
+import matplotlib.pylab as plt
+
+FILENAME=sys.argv[1]
+Nx=int(sys.argv[2])
+Ny=int(sys.argv[3])
+Nz=int(sys.argv[4])
+
+# read the input image
+Output = np.fromfile(FILENAME,dtype = np.uint8)
+Output.shape = (Nz,Ny,Nx)
+
+Oil=np.count_nonzero(Output==1)
+Water=np.count_nonzero(Output==2)
+Sw=Water/(Oil+Water)
+
+Porosity=1.0-(Oil+Water)/(Nx*Ny*Nz)
+
+print(FILENAME,"Porosity=", Porosity)
+
+SaturationProfile=np.zeros(Nz)
+PorosityProfile=np.zeros(Nz)
+# Compute saturation slice by slice 
+for idx in range(0, Nz):
+   Slice = Output[idx,:,:]
+   Oil=np.count_nonzero(Slice==1)
+   Water=np.count_nonzero(Slice==2)
+   SaturationProfile[idx]=Water/(Oil+Water)
+   PorosityProfile[idx]=(Oil+Water)/(Nx*Ny)
+   
+
+plt.figure()
+plt.plot(SaturationProfile)
+plt.xlabel('Position (z)')
+plt.ylabel('Water Saturation')
+plt.show()

example/Workflow/HelperFunctions.R

@@ -2,16 +2,57 @@ require("ggplot2")
 
 GG_THEME=theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
 
+ReadDatabase<-function(FILE){
+
+    INPUT<-gsub(';','',readLines(FILE))
+
+    S<-gsub('tauA = ','',gsub("\\s+"," ",(grep("tauA",INPUT,value=TRUE))))
+    TAU_A = as.numeric(gsub("/.*","",S))
+    S<-gsub('tauB = ','',gsub("\\s+"," ",(grep("tauB",INPUT,value=TRUE))))
+    TAU_B = as.numeric(gsub("/.*","",S))
+    S<-gsub('rhoA = ','',gsub("\\s+"," ",(grep("rhoA",INPUT,value=TRUE))))
+    RHO_A = as.numeric(gsub("/.*","",S))
+    S<-gsub('rhoB = ','',gsub("\\s+"," ",(grep("rhoB",INPUT,value=TRUE))))
+    RHO_B = as.numeric(gsub("/.*","",S))
+
+    S<-gsub('alpha = ','',gsub("\\s+"," ",(grep("alpha",INPUT,value=TRUE))))
+    ALPHA = as.numeric(gsub("/.*","",S))
+
+    # Read the affinity
+    S<-gsub('ComponentAffinity = ','',gsub("\\s+"," ",(grep("ComponentAffinity",INPUT,value=TRUE))))
+    S<-gsub("/.*","",S)
+    AFFINITY<-as.numeric(unlist(strsplit(S,", ")))
+
+    PARAMETERS<-c(TAU_A,TAU_B,RHO_A,RHO_B,ALPHA,AFFINITY)
+    
+    return(PARAMETERS)
+}
+
 ReadSubphase<-function(PATH){
   FILE=paste0(PATH,"/subphase.csv")
   S<-read.csv(FILE,head=TRUE,sep=" ")
   S$Vw<-S$Vwc+S$Vwd
   S$Vn<-S$Vnc+S$Vnd
+  S$Aw<-S$Awc+S$Awd
+  S$An<-S$Anc+S$And
+  S$Hw<-S$Hwc+S$Hwd
+  S$Hn<-S$Hnc+S$Hnd
+  S$Xw<-S$Xwc+S$Xwd
+  S$Xn<-S$Xnc+S$Xnd
+
   S$Sw<-S$Vw/(S$Vn+S$Vw)
+  S$pw<-(S$pwc*S$Vwc+S$pwd*S$Vwd) / (S$Vwc+S$Vwd)
+  S$pn<-(S$pnc*S$Vnc+S$pnd*S$Vnd) / (S$Vnc+S$Vnd)
+
   S$Qwx<-S$Vw*(S$Pwc_x+S$Pwd_x)/(S$Mwc+S$Mwd)
   S$Qnx<-S$Vn*(S$Pnc_x+S$Pnd_x)/(S$Mnc+S$Mnd)
-  S$Krn<-S$nun*S$Qnx/S$Fx
-  S$Krw<-S$nuw*S$Qwx/S$Fx
+  S$Qwy<-S$Vw*(S$Pwc_y+S$Pwd_y)/(S$Mwc+S$Mwd)
+  S$Qny<-S$Vn*(S$Pnc_y+S$Pnd_y)/(S$Mnc+S$Mnd)
+  S$Qwz<-S$Vw*(S$Pwc_z+S$Pwd_z)/(S$Mwc+S$Mwd)
+  S$Qnz<-S$Vn*(S$Pnc_z+S$Pnd_z)/(S$Mnc+S$Mnd)
+
+  S$Krn<-S$nun*S$Qnz/S$Fz
+  S$Krw<-S$nuw*S$Qwz/S$Fz
   S$Case<-PATH
   return(S)
 }

gpu/Color.cu

@@ -143,7 +143,7 @@ __global__  void dvc_ScaLBL_Color_InitDistance(char *ID, double *Den, double *Ph
 
 __global__  void dvc_ScaLBL_Color_BC(int *list, int *Map, double *Phi, double *Den, double vA, double vB, int count, int Np)
 {
-	int idx,n,nm;
+	int idx,n;
 	// Fill the outlet with component b
 	idx = blockIdx.x*blockDim.x + threadIdx.x;
 	if (idx < count){
@@ -1255,6 +1255,15 @@ __global__  void dvc_ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny
 }
 
 
+__global__  void dvc_ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int Source, int Dest){
+	double value;
+	int n =  blockIdx.x*blockDim.x + threadIdx.x;
+	if (n < Nx*Ny){
+		value = Phi[Source*Nx*Ny+n];
+		Phi[Dest*Nx*Ny+n] = value;
+	}
+}
+
 
 __global__  void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
 		double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
@@ -3486,13 +3495,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_ColorMass(double *Aq, double *Bq, double
 		double *Velocity, double *ColorGrad, double beta, int start, int finish, int Np){
 
 	int n;
-	double fq;
 	// non-conserved moments
 	double nA,nB; // number density
 	double a1,b1,a2,b2,nAB,delta;
 	double C,nx,ny,nz; //color gradient magnitude and direction
 	double ux,uy,uz;
-	double phi,tau,rho0,rlx_setA,rlx_setB;
 
 	int S = Np/NBLOCKS/NTHREADS + 1;
 	for (int s=0; s<S; s++){
@@ -3580,13 +3587,11 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_ColorMass(int *neighborList, double *Aq,
 		double *Velocity, double *ColorGrad, double beta, int start, int finish, int Np){
 
 	int n,nread;
-	double fq;
 	// non-conserved moments
 	double nA,nB; // number density
 	double a1,b1,a2,b2,nAB,delta;
 	double C,nx,ny,nz; //color gradient magnitude and direction
 	double ux,uy,uz;
-	double phi,tau,rho0,rlx_setA,rlx_setB;
 
 	int S = Np/NBLOCKS/NTHREADS + 1;
 	for (int s=0; s<S; s++){
@@ -4153,5 +4158,9 @@ extern "C" void ScaLBL_Color_BC_Z(int *list, int *Map, double *Phi, double *Den,
 	}
 }
 
+extern "C" void ScaLBL_CopySlice_z(double *Phi, int Nx, int Ny, int Nz, int Source, int Dest){
+	int GRID = Nx*Ny / 512 + 1;
+	dvc_ScaLBL_CopySlice_z<<<GRID,512>>>(Phi,Nx,Ny,Nz,Source,Dest);
+}
 
 

gpu/D3Q19.cu

@@ -282,6 +282,36 @@ __global__ void dvc_ScaLBL_D3Q19_Init(double *dist, int Np)
 	}
 }
 
+__global__ void dvc_ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den)
+{
+	int n;
+	int S = Np/NBLOCKS/NTHREADS + 1;
+	for (int s=0; s<S; s++){
+		//........Get 1-D index for this thread....................
+		n = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
+		if (n<Np ){
+			dist[n] = Den - 0.6666666666666667;
+			dist[Np+n] = 0.055555555555555555;		//double(100*n)+1.f;
+			dist[2*Np+n] = 0.055555555555555555;	//double(100*n)+2.f;
+			dist[3*Np+n] = 0.055555555555555555;	//double(100*n)+3.f;
+			dist[4*Np+n] = 0.055555555555555555;	//double(100*n)+4.f;
+			dist[5*Np+n] = 0.055555555555555555;	//double(100*n)+5.f;
+			dist[6*Np+n] = 0.055555555555555555;	//double(100*n)+6.f;
+			dist[7*Np+n] = 0.0277777777777778;   //double(100*n)+7.f;
+			dist[8*Np+n] = 0.0277777777777778;   //double(100*n)+8.f;
+			dist[9*Np+n] = 0.0277777777777778;   //double(100*n)+9.f;
+			dist[10*Np+n] = 0.0277777777777778;  //double(100*n)+10.f;
+			dist[11*Np+n] = 0.0277777777777778;  //double(100*n)+11.f;
+			dist[12*Np+n] = 0.0277777777777778;  //double(100*n)+12.f;
+			dist[13*Np+n] = 0.0277777777777778;  //double(100*n)+13.f;
+			dist[14*Np+n] = 0.0277777777777778;  //double(100*n)+14.f;
+			dist[15*Np+n] = 0.0277777777777778;  //double(100*n)+15.f;
+			dist[16*Np+n] = 0.0277777777777778;  //double(100*n)+16.f;
+			dist[17*Np+n] = 0.0277777777777778;  //double(100*n)+17.f;
+			dist[18*Np+n] = 0.0277777777777778;  //double(100*n)+18.f;
+		}
+	}
+}
 
 //*************************************************************************
 __global__  void dvc_ScaLBL_D3Q19_Swap_Compact(int *neighborList, double *disteven, double *distodd, int Np, int q){
@@ -1553,7 +1583,6 @@ __global__  void dvc_ScaLBL_D3Q19_Momentum(double *dist, double *vel, int N)
 	double f1,f2,f3,f4,f5,f6,f7,f8,f9;
 	double f10,f11,f12,f13,f14,f15,f16,f17,f18;
 	double vx,vy,vz;
-	char id;
 
 	int S = N/NBLOCKS/NTHREADS + 1;
 	for (int s=0; s<S; s++){
@@ -1744,6 +1773,43 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_Pressure_BC_Z(int *list, double *dist,
 		//...................................................
 	}
 }
+__global__  void dvc_ScaLBL_D3Q19_Reflection_BC_z(int *list, double *dist, int count, int Np){
+	int idx, n;
+	idx = blockIdx.x*blockDim.x + threadIdx.x;
+	if (idx < count){
+		n = list[idx];
+		double f5 = 0.111111111111111111111111 - dist[6*Np+n];
+		double f11 = 0.05555555555555555555556 - dist[12*Np+n];
+		double f14 = 0.05555555555555555555556 - dist[13*Np+n];
+		double f15 = 0.05555555555555555555556 - dist[16*Np+n];
+		double f18 = 0.05555555555555555555556 - dist[17*Np+n];
+		
+		dist[6*Np+n] = f5;
+		dist[12*Np+n] = f11;
+		dist[13*Np+n] = f14;
+		dist[16*Np+n] = f15;
+		dist[17*Np+n] = f18;
+	}
+}
+
+__global__  void dvc_ScaLBL_D3Q19_Reflection_BC_Z(int *list, double *dist, int count, int Np){
+	int idx, n;
+	idx = blockIdx.x*blockDim.x + threadIdx.x;
+	if (idx < count){
+		n = list[idx];
+		double f6 = 0.111111111111111111111111 - dist[5*Np+n];
+		double f12 = 0.05555555555555555555556 - dist[11*Np+n];
+		double f13 = 0.05555555555555555555556 - dist[14*Np+n] ;
+		double f16 = 0.05555555555555555555556 - dist[15*Np+n];
+		double f17 = 0.05555555555555555555556 - dist[18*Np+n];
+		
+		dist[5*Np+n] = f6;
+		dist[11*Np+n] = f12;
+		dist[14*Np+n] = f13;
+		dist[15*Np+n] = f16;
+		dist[18*Np+n] = f17;
+	}
+}
 
 __global__  void dvc_ScaLBL_D3Q19_AAodd_Pressure_BC_z(int *d_neighborList, int *list, double *dist, double din, int count, int Np)
 {
@@ -2340,7 +2406,15 @@ extern "C" void ScaLBL_D3Q19_Init(double *dist, int Np){
 	dvc_ScaLBL_D3Q19_Init<<<NBLOCKS,NTHREADS >>>(dist, Np);
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){
-		printf("CUDA error in ScaLBL_D3Q19_AA_Init: %s \n",cudaGetErrorString(err));
+		printf("CUDA error in ScaLBL_D3Q19_Init: %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q19_GreyIMRT_Init(double *dist, int Np, double Den){
+	dvc_ScaLBL_D3Q19_GreyIMRT_Init<<<NBLOCKS,NTHREADS >>>(dist, Np, Den);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_GreyIMRT_Init: %s \n",cudaGetErrorString(err));
 	}
 }
 
@@ -2630,11 +2704,23 @@ extern "C" double deviceReduce(double *in, double* out, int N) {
 	return sum;
 }
 
-//
-//extern "C" void ScaLBL_D3Q19_Pressure_BC_Z(int *list, double *dist, double dout, int count, int Np){
-//	int GRID = count / 512 + 1;
-//	dvc_ScaLBL_D3Q19_Pressure_BC_Z<<<GRID,512>>>(disteven, distodd, dout, Nx, Ny, Nz, outlet);
-//}
+extern "C" void ScaLBL_D3Q19_Reflection_BC_z(int *list, double *dist, int count, int Np){
+	int GRID = count / 512 + 1;
+	dvc_ScaLBL_D3Q19_Reflection_BC_z<<<GRID,512>>>(list, dist, count, Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_Reflection_BC_z (kernel): %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q19_Reflection_BC_Z(int *list, double *dist, int count, int Np){
+	int GRID = count / 512 + 1;
+	dvc_ScaLBL_D3Q19_Reflection_BC_Z<<<GRID,512>>>(list, dist, count, Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q19_Reflection_BC_Z (kernel): %s \n",cudaGetErrorString(err));
+	}
+}
 
 extern "C" void ScaLBL_D3Q19_AAeven_MRT(double *dist, int start, int finish, int Np, double rlx_setA, double rlx_setB, double Fx,
        double Fy, double Fz){

gpu/D3Q7.cu

@@ -14,6 +14,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 // GPU Functions for D3Q7 Lattice Boltzmann Methods
+#include <stdio.h>
 
 #define NBLOCKS 560
 #define NTHREADS 128
@@ -94,6 +95,25 @@ __global__ void dvc_ScaLBL_D3Q7_Unpack(int q,  int *list,  int start, int count,
 	}
 }
 
+__global__  void dvc_ScaLBL_D3Q7_Reflection_BC_z(int *list, double *dist, int count, int Np){
+	int idx, n;
+	idx = blockIdx.x*blockDim.x + threadIdx.x;
+	if (idx < count){
+		n = list[idx];
+		double f5 = 0.222222222222222222222222 - dist[6*Np+n];
+		dist[6*Np+n] = f5;
+	}
+}
+
+__global__  void dvc_ScaLBL_D3Q7_Reflection_BC_Z(int *list, double *dist, int count, int Np){
+	int idx, n;
+	idx = blockIdx.x*blockDim.x + threadIdx.x;
+	if (idx < count){
+		n = list[idx];
+		double f6 = 0.222222222222222222222222 - dist[5*Np+n];
+		dist[5*Np+n] = f6;
+	}
+}
 __global__ void dvc_ScaLBL_D3Q7_Init(char *ID, double *f_even, double *f_odd, double *Den, int Nx, int Ny, int Nz)
 {
 	int n,N;
@@ -222,6 +242,24 @@ __global__  void dvc_ScaLBL_D3Q7_Density(char *ID, double *disteven, double *dis
 	}
 }
 
+extern "C" void ScaLBL_D3Q7_Reflection_BC_z(int *list, double *dist, int count, int Np){
+	int GRID = count / 512 + 1;
+	dvc_ScaLBL_D3Q7_Reflection_BC_z<<<GRID,512>>>(list, dist, count, Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q7_Reflection_BC_z (kernel): %s \n",cudaGetErrorString(err));
+	}
+}
+
+extern "C" void ScaLBL_D3Q7_Reflection_BC_Z(int *list, double *dist, int count, int Np){
+	int GRID = count / 512 + 1;
+	dvc_ScaLBL_D3Q7_Reflection_BC_Z<<<GRID,512>>>(list, dist, count, Np);
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err){
+		printf("CUDA error in ScaLBL_D3Q7_Reflection_BC_Z (kernel): %s \n",cudaGetErrorString(err));
+	}
+}
+
 extern "C" void ScaLBL_D3Q7_Unpack(int q, int *list,  int start, int count, double *recvbuf, double *dist, int N){
 	int GRID = count / 512 + 1;
 	dvc_ScaLBL_D3Q7_Unpack <<<GRID,512 >>>(q, list, start, count, recvbuf, dist, N);

models/ColorModel.cpp

@@ -19,6 +19,7 @@ color lattice boltzmann model
 #include "models/ColorModel.h"
 #include "analysis/distance.h"
 #include "analysis/morphology.h"
+#include "common/Communication.h"
 #include "common/ReadMicroCT.h"
 #include <stdlib.h>
 #include <time.h>
@@ -70,8 +71,6 @@ void ScaLBL_ColorModel::ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq
     File.close();
 }
  */
-
-
 void ScaLBL_ColorModel::ReadParams(string filename){
 	// read the input database 
 	db = std::make_shared<Database>( filename );
@@ -144,21 +143,21 @@ void ScaLBL_ColorModel::ReadParams(string filename){
 	// Override user-specified boundary condition for specific protocols
 	auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
 	if (protocol == "seed water"){
-		if (BoundaryCondition != 0 ){
+		if (BoundaryCondition != 0 && BoundaryCondition != 5){
 			BoundaryCondition = 0;
 			if (rank==0) printf("WARNING: protocol (seed water) supports only full periodic boundary condition \n");
 		}
 		domain_db->putScalar<int>( "BC", BoundaryCondition );
 	}
 	else if (protocol == "open connected oil"){
-		if (BoundaryCondition != 0 ){
+		if (BoundaryCondition != 0 && BoundaryCondition != 5){
 			BoundaryCondition = 0;
 			if (rank==0) printf("WARNING: protocol (open connected oil) supports only full periodic boundary condition \n");
 		}
 		domain_db->putScalar<int>( "BC", BoundaryCondition );
 	}
 	else if (protocol == "shell aggregation"){
-		if (BoundaryCondition != 0 ){
+		if (BoundaryCondition != 0 && BoundaryCondition != 5){
 			BoundaryCondition = 0;
 			if (rank==0) printf("WARNING: protocol (shell aggregation) supports only full periodic boundary condition \n");
 		}
@@ -200,14 +199,22 @@ void ScaLBL_ColorModel::ReadInput(){
 	if (color_db->keyExists( "image_sequence" )){
 		auto ImageList = color_db->getVector<std::string>( "image_sequence");
 		int IMAGE_INDEX = color_db->getWithDefault<int>( "image_index", 0 );
-		int IMAGE_COUNT = ImageList.size();
 		std::string first_image = ImageList[IMAGE_INDEX];
 		Mask->Decomp(first_image);
 		IMAGE_INDEX++;
 	}
 	else if (domain_db->keyExists( "GridFile" )){
+        // Read the local domain data
 	    auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD );
-	    for (int i=0; i<Nx*Ny*Nz; i++) Mask->id[i] = input_id(i);
+        // Fill the halo (assuming GCW of 1)
+        array<int,3> size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) };
+        ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz };
+        ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
+        fillHalo<signed char> fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
+        Array<signed char> id_view;
+        id_view.viewRaw( size1, Mask->id );
+        fill.copy( input_id, id_view );
+        fill.fill( id_view );
 	}
 	else if (domain_db->keyExists( "Filename" )){
 		auto Filename = domain_db->getScalar<std::string>( "Filename" );
@@ -221,7 +228,6 @@ void ScaLBL_ColorModel::ReadInput(){
 	// Generate the signed distance map
 	// Initialize the domain and communication
 	Array<char> id_solid(Nx,Ny,Nz);
-	int count = 0;
 	// Solve for the position of the solid phase
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
@@ -238,7 +244,6 @@ void ScaLBL_ColorModel::ReadInput(){
 	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;
 				// Initialize distance to +/- 1
 				Averages->SDs(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
 			}
@@ -271,7 +276,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 	double label_count_global[NLABELS];
 	// Assign the labels
 
-	for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
+	for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
 
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
@@ -299,7 +304,8 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 	// Set Dm to match Mask
 	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; 
 	
-	for (int idx=0; idx<NLABELS; idx++)		label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
+	for (size_t idx=0; idx<NLABELS; idx++)
+		label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
 
 	if (rank==0){
 		printf("Component labels: %lu \n",NLABELS);
@@ -378,16 +384,16 @@ void ScaLBL_ColorModel::Create(){
 	}
 	// check that TmpMap is valid
 	for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
-		int n = TmpMap[idx];
+		auto n = TmpMap[idx];
 		if (n > Nx*Ny*Nz){
-			printf("Bad value! idx=%i \n");
+			printf("Bad value! idx=%i \n", n);
 			TmpMap[idx] = Nx*Ny*Nz-1;
 		}
 	}
 	for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
-		int n = TmpMap[idx];
+		auto n = TmpMap[idx];
 		if ( n > Nx*Ny*Nz ){
-			printf("Bad value! idx=%i \n");
+			printf("Bad value! idx=%i \n",n);
 			TmpMap[idx] = Nx*Ny*Nz-1;
 		}
 	}
@@ -479,7 +485,8 @@ void ScaLBL_ColorModel::Initialize(){
 	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
 	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 
-	if (BoundaryCondition >0 ){
+	// establish reservoirs for external bC
+	if (BoundaryCondition == 1 || BoundaryCondition == 2 ||  BoundaryCondition == 3 || BoundaryCondition == 4 ){
 		if (Dm->kproc()==0){
 			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
@@ -502,6 +509,7 @@ void ScaLBL_ColorModel::Run(){
 	int IMAGE_COUNT = 0;
 	std::vector<std::string> ImageList;
 	bool SET_CAPILLARY_NUMBER = false;
+	bool RESCALE_FORCE = false;
 	bool MORPH_ADAPT = false;
 	bool USE_MORPH = false;
 	bool USE_SEED = false;
@@ -510,6 +518,7 @@ void ScaLBL_ColorModel::Run(){
 	int MAX_MORPH_TIMESTEPS = 50000; // maximum number of LBM timesteps to spend in morphological adaptation routine
 	int MIN_STEADY_TIMESTEPS = 100000;
 	int MAX_STEADY_TIMESTEPS = 200000;
+	int RESCALE_FORCE_AFTER_TIMESTEP = 0;
 	int RAMP_TIMESTEPS = 0;//50000;		 // number of timesteps to run initially (to get a reasonable velocity field before other pieces kick in)
 	int CURRENT_MORPH_TIMESTEPS=0;   // counter for number of timesteps spent in  morphological adaptation routine (reset each time)
 	int CURRENT_STEADY_TIMESTEPS=0;   // counter for number of timesteps spent in  morphological adaptation routine (reset each time)
@@ -524,13 +533,25 @@ void ScaLBL_ColorModel::Run(){
 	double initial_volume = 0.0;
 	double delta_volume = 0.0;
 	double delta_volume_target = 0.0;
-	double RESIDUAL_ENDPOINT_THRESHOLD = 0.04;
-	double NOISE_THRESHOLD = 0.0;
-	double BUMP_RATE = 2.0;
-	bool USE_BUMP_RATE = false;
-	int RESCALE_FORCE_COUNT = 0;
-	int RESCALE_FORCE_MAX = 0;
 	
+	/* history for morphological algoirthm */
+	double KRA_MORPH_FACTOR=0.5;
+	double volA_prev = 0.0; 
+	double log_krA_prev = 1.0;
+	double log_krA_target = 1.0;
+	double log_krA = 1.0;
+	double slope_krA_volume = 0.0;
+	if (color_db->keyExists( "vol_A_previous" )){
+		volA_prev  = color_db->getScalar<double>( "vol_A_previous" );
+	}
+	if (color_db->keyExists( "log_krA_previous" )){
+		log_krA_prev  = color_db->getScalar<double>( "log_krA_previous" );
+	}
+	if (color_db->keyExists( "krA_morph_factor" )){
+		KRA_MORPH_FACTOR  = color_db->getScalar<double>( "krA_morph_factor" );
+	}
+	
+	/* defaults for simulation protocols */
 	auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
 	if (protocol == "image sequence"){
 		// Get the list of images
@@ -542,45 +563,33 @@ void ScaLBL_ColorModel::Run(){
 		USE_MORPH = true;
 	}
 	else if (protocol == "seed water"){
-		morph_delta = 0.05;
+		morph_delta = -0.05;
 		seed_water = 0.01;
 		USE_SEED = true;
 		USE_MORPH = true;
 	}
 	else if (protocol == "open connected oil"){
-		morph_delta = 0.05;
+		morph_delta = -0.05;
 		USE_MORPH = true;
 		USE_MORPHOPEN_OIL = true;
 	}
 	else if (protocol == "shell aggregation"){
-		morph_delta = 0.05;
+		morph_delta = -0.05;
 		USE_MORPH = true;
 	}  
-	
-	if (color_db->keyExists( "residual_endpoint_threshold" )){
-		RESIDUAL_ENDPOINT_THRESHOLD  = color_db->getScalar<double>( "residual_endpoint_threshold" );
-	}
-	if (color_db->keyExists( "noise_threshold" )){
-		NOISE_THRESHOLD  = color_db->getScalar<double>( "noise_threshold" );
-		USE_BUMP_RATE = true;
-	}
-	if (color_db->keyExists( "bump_rate" )){
-		BUMP_RATE  = color_db->getScalar<double>( "bump_rate" );
-		USE_BUMP_RATE = true;
-	}
 	if (color_db->keyExists( "capillary_number" )){
 		capillary_number = color_db->getScalar<double>( "capillary_number" );
 		SET_CAPILLARY_NUMBER=true;
-		//RESCALE_FORCE_MAX = 1;
 	}
-	if (analysis_db->keyExists( "rescale_force_count" )){
-		RESCALE_FORCE_MAX = analysis_db->getScalar<int>( "rescale_force_count" );
+	if (color_db->keyExists( "rescale_force_after_timestep" )){
+		RESCALE_FORCE_AFTER_TIMESTEP = color_db->getScalar<int>( "rescale_force_after_timestep" );
+		RESCALE_FORCE = true;
 	}
 	if (color_db->keyExists( "timestep" )){
 		timestep = color_db->getScalar<int>( "timestep" );
 	}
-	if (BoundaryCondition != 0 && SET_CAPILLARY_NUMBER==true){
-		if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 \n");
+	if (BoundaryCondition != 0 && BoundaryCondition != 5 && SET_CAPILLARY_NUMBER==true){
+		if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 or 5 \n");
 		SET_CAPILLARY_NUMBER=false;
 	}
 	if (analysis_db->keyExists( "seed_water" )){
@@ -683,7 +692,7 @@ void ScaLBL_ColorModel::Run(){
 
 		// Perform the collision operation
 		ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
-		if (BoundaryCondition > 0){
+		if (BoundaryCondition > 0 && BoundaryCondition < 5){
 			ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
 			ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
 		}
@@ -704,9 +713,14 @@ void ScaLBL_ColorModel::Run(){
 			din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
 			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
 		}
+		else if (BoundaryCondition == 5){
+			ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
+			ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
+		}
 		ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
 				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
-		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+		ScaLBL_DeviceBarrier(); 
+		MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
 
 		// *************EVEN TIMESTEP*************
 		timestep++;
@@ -720,7 +734,7 @@ void ScaLBL_ColorModel::Run(){
 		// Perform the collision operation
 		ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
 		// Halo exchange for phase field
-		if (BoundaryCondition > 0){
+		if (BoundaryCondition > 0 && BoundaryCondition < 5){
 			ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
 			ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
 		}
@@ -739,21 +753,23 @@ void ScaLBL_ColorModel::Run(){
 			din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
 			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
 		}
+		else if (BoundaryCondition == 5){
+			ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
+			ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
+		}
 		ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
 				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
-		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+		ScaLBL_DeviceBarrier(); 
+		MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
 		//************************************************************************
-		
-		MPI_Barrier(comm);
 		PROFILE_STOP("Update");
 
-		if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){
+		if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition == 4){
 			printf("%i %f \n",timestep,din);
 		}
 		// Run the analysis
 		analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
 
-		
 		// allow initial ramp-up to get closer to steady state
 		if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
 			analysis.finish();
@@ -763,7 +779,7 @@ void ScaLBL_ColorModel::Run(){
 			double volA = Averages->gnb.V; 
 			volA /= Dm->Volume;
 			volB /= Dm->Volume;;
-			initial_volume = volA*Dm->Volume;
+			//initial_volume = volA*Dm->Volume;
 			double vA_x = Averages->gnb.Px/Averages->gnb.M; 
 			double vA_y = Averages->gnb.Py/Averages->gnb.M; 
 			double vA_z = Averages->gnb.Pz/Averages->gnb.M; 
@@ -795,28 +811,46 @@ void ScaLBL_ColorModel::Run(){
 					isSteady = true;
 				if (CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS)
 					isSteady = true;
-				
-				if (SET_CAPILLARY_NUMBER  && RESCALE_FORCE_COUNT < RESCALE_FORCE_MAX){
-					RESCALE_FORCE_COUNT++;
+				if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_STEADY_TIMESTEPS > RESCALE_FORCE_AFTER_TIMESTEP){
+					RESCALE_FORCE = false;
 					Fx *= capillary_number / Ca;
 					Fy *= capillary_number / Ca;
 					Fz *= capillary_number / Ca;
-
 					if (force_mag > 1e-3){
 						Fx *= 1e-3/force_mag;   // impose ceiling for stability
 						Fy *= 1e-3/force_mag;   
 						Fz *= 1e-3/force_mag;   
 					}
-					
 					if (rank == 0) printf("    -- adjust force by factor %f \n ",capillary_number / Ca);
 					Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
 					color_db->putVector<double>("F",{Fx,Fy,Fz});
 				}
-
 				if ( isSteady ){
 					MORPH_ADAPT = true;
 					CURRENT_MORPH_TIMESTEPS=0;
 					delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
+					//****** ENDPOINT ADAPTATION ********/
+					double krA_TMP= fabs(muA*flow_rate_A / force_mag);
+					double krB_TMP= fabs(muB*flow_rate_B / force_mag);
+					log_krA = log(krA_TMP);
+					if (krA_TMP < 0.0){
+						// cannot do endpoint adaptation if kr is negative
+						log_krA = log_krA_prev;
+					}
+					else if (krA_TMP < krB_TMP && morph_delta > 0.0){
+						/** morphological target based on relative permeability for A **/
+						log_krA_target = log(KRA_MORPH_FACTOR*(krA_TMP));
+						slope_krA_volume = (log_krA - log_krA_prev)/(Dm->Volume*(volA - volA_prev));
+						delta_volume_target=min(delta_volume_target,Dm->Volume*(volA+(log_krA_target - log_krA)/slope_krA_volume));
+						if (rank==0){
+							printf("    Enabling endpoint adaptation: krA = %f, krB = %f \n",krA_TMP,krB_TMP);	
+							printf("    log(kr)=%f, volume=%f, TARGET log(kr)=%f, volume change=%f \n",log_krA, volA, log_krA_target, delta_volume_target/(volA*Dm->Volume));							
+						}
+					}
+					log_krA_prev = log_krA;
+					volA_prev = volA;
+					//******************************** **/
+					/**  compute averages & write data **/
 					Averages->Full();
 					Averages->Write(timestep);
 					analysis.WriteVisData(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
@@ -879,7 +913,7 @@ void ScaLBL_ColorModel::Run(){
 							WriteHeader=true;
 						kr_log_file = fopen("relperm.csv","a");
 						if (WriteHeader)
-							fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility);
+							fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected cap.pressure cap.pressure.connected pressure.drop Ca M\n");
 
 						fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility);
 						fclose(kr_log_file);
@@ -890,26 +924,16 @@ void ScaLBL_ColorModel::Run(){
 						Fx *= capillary_number / Ca;
 						Fy *= capillary_number / Ca;
 						Fz *= capillary_number / Ca;
-						RESCALE_FORCE_COUNT = 1;
 						if (force_mag > 1e-3){
 							Fx *= 1e-3/force_mag;   // impose ceiling for stability
 							Fy *= 1e-3/force_mag;   
 							Fz *= 1e-3/force_mag;   
 						}
-						if (flow_rate_A < NOISE_THRESHOLD && USE_BUMP_RATE){
-							if (rank==0) printf("Hit noise threshold (%f): bumping capillary number by %f X \n",NOISE_THRESHOLD,BUMP_RATE);
-							Fx *= BUMP_RATE;   // impose bump condition
-							Fy *= BUMP_RATE;   
-							Fz *= BUMP_RATE;   
-							capillary_number *= BUMP_RATE;
-							color_db->putScalar<int>("capillary_number",capillary_number);
-							current_db->putDatabase("Color", color_db);
-							MORPH_ADAPT = false; // re-run current point if below noise threshold
-						}
 						if (rank == 0) printf("    -- adjust force by factor %f \n ",capillary_number / Ca);
 						Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
 						color_db->putVector<double>("F",{Fx,Fy,Fz});
 					}
+					
 					CURRENT_STEADY_TIMESTEPS = 0;
 				}
 				else{
@@ -943,7 +967,7 @@ void ScaLBL_ColorModel::Run(){
 					delta_volume = volA*Dm->Volume - initial_volume;
 					CURRENT_MORPH_TIMESTEPS += analysis_interval;
 					double massChange = SeedPhaseField(seed_water);
-					if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target);
+					if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", massChange, delta_volume, delta_volume_target);
 				}
 				else if (USE_MORPHOPEN_OIL){
 					delta_volume = volA*Dm->Volume - initial_volume;
@@ -961,41 +985,29 @@ void ScaLBL_ColorModel::Run(){
 					CURRENT_STEADY_TIMESTEPS=0;
 					initial_volume = volA*Dm->Volume;
 					delta_volume = 0.0;
-					if (USE_DIRECT){
-						//BoundaryCondition = 0;
-						//ScaLBL_Comm->BoundaryCondition = 0;
-						//ScaLBL_Comm_Regular->BoundaryCondition = 0;
-						//Fx = capillary_number*dir_x*force_mag / Ca;
-						//Fy = capillary_number*dir_y*force_mag / Ca;
-						//Fz = capillary_number*dir_z*force_mag / Ca;
-					}
+					if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
+						RESCALE_FORCE = true;
 				}
 				else if (!(USE_DIRECT) && CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
 					MORPH_ADAPT = false;
 					CURRENT_STEADY_TIMESTEPS=0;
 					initial_volume = volA*Dm->Volume;
 					delta_volume = 0.0;
+					RESCALE_FORCE = true;
+					if (RESCALE_FORCE_AFTER_TIMESTEP > 0)
+						RESCALE_FORCE = true;
 				}
-				if ( REVERSE_FLOW_DIRECTION ){
-					//if (rank==0) printf("*****REVERSE FLOW DIRECTION***** \n");
-					delta_volume = 0.0;
-					// flow direction will reverse after next steady point
-					MORPH_ADAPT = false;
-					CURRENT_STEADY_TIMESTEPS=0;
-					//morph_delta *= (-1.0);
-					REVERSE_FLOW_DIRECTION = false;
-				}
-				MPI_Barrier(comm);
 			}
 			morph_timesteps += analysis_interval;
 		}
+		MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
 	}
 	analysis.finish();
 	PROFILE_STOP("Loop");
 	PROFILE_SAVE("lbpm_color_simulator",1);
 	//************************************************************************
 	ScaLBL_DeviceBarrier();
-	MPI_Barrier(comm);
+	MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
 	stoptime = MPI_Wtime();
 	if (rank==0) printf("-------------------------------------------------------------------\n");
 	// Compute the walltime per timestep
@@ -1015,7 +1027,6 @@ void ScaLBL_ColorModel::Run(){
 
 double ScaLBL_ColorModel::ImageInit(std::string Filename){
 	
-	bool suppress = false;
 	if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str());
 	Mask->Decomp(Filename);
 	for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i];  // save what was read
@@ -1046,12 +1057,12 @@ double ScaLBL_ColorModel::ImageInit(std::string Filename){
 	
 	if (rank==0) printf("   new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
 	ScaLBL_CopyToDevice(Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double));
-	MPI_Barrier(comm);
+	MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
 	
 	ScaLBL_D3Q19_Init(fq, Np);
 	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
 	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
-	MPI_Barrier(comm);
+	MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
 	
 	ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
 
@@ -1085,10 +1096,9 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
 		ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm);
 		MPI_Barrier(Dm->Comm);
 
-		int count_oil=0;
-		int count_connected=0;
-		int count_porespace=0;
-		int count_water=0;
+		long long count_connected=0;
+		long long count_porespace=0;
+		long long count_water=0;
 		for (int k=1; k<nz-1; k++){
 			for (int j=1; j<ny-1; j++){
 				for (int i=1; i<nx-1; i++){
@@ -1188,7 +1198,7 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
 		ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
 		ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
 		ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
-		if (BoundaryCondition >0 ){
+		if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
 			if (Dm->kproc()==0){
 				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
@@ -1203,7 +1213,6 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
 	}
 	return(volume_change);
 }
-
 double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
   srand(time(NULL));
   double mass_loss =0.f;
@@ -1216,26 +1225,7 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
   ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double));
   ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double));
   
-/*	for (int k=1; k<Nz-1; k++){
-		for (int j=1; j<Ny-1; j++){
-			for (int i=1; i<Nx-1; i++){
-				double random_value = double(rand())/ RAND_MAX;
  
-				if (Averages->SDs(i,j,k) < 0.f){
-					// skip
-				}
-				else if (phase(i,j,k) > 0.f ){
-					phase(i,j,k) -= random_value*seed_water_in_oil;
-					mass_loss += random_value*seed_water_in_oil;
-					count++;
-				}
-				else {
-
-				}
-			}
-		}
-	}
-	*/
    for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){
     double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
     double dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
@@ -1304,6 +1294,8 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	double vF = 0.f;
 	double vS = 0.f;
 	double delta_volume;
+	double WallFactor = 0.0;
+	bool USE_CONNECTED_NWP = false;
 
 	DoubleArray phase(Nx,Ny,Nz);
 	IntArray phase_label(Nx,Ny,Nz);;
@@ -1316,8 +1308,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	// 1. Copy phase field to CPU
 	ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
 
-	double count,count_global,volume_initial,volume_final,volume_connected;
-	count = 0.f;
+	double count = 0.f;
 	for (int k=1; k<Nz-1; k++){
 		for (int j=1; j<Ny-1; j++){
 			for (int i=1; i<Nx-1; i++){
@@ -1325,7 +1316,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 			}
 		}
 	}
-	volume_initial = sumReduce( Dm->Comm, count);
+	double volume_initial = sumReduce( Dm->Comm, count);
 	/*
 	sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank);
 	FILE *INPUT = fopen(LocalRankFilename,"wb");
@@ -1333,13 +1324,16 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	fclose(INPUT);
 	*/
 	// 2. Identify connected components of phase field -> phase_label
+	
+	double volume_connected = 0.0;
+       	double second_biggest = 0.0;
+	if (USE_CONNECTED_NWP){
 		BlobIDstruct new_index;
 		ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
-	MPI_Barrier(comm);
+		MPI_Barrier(Dm->Comm);
 
 		// only operate on component "0"
 		count = 0.0;
-	double second_biggest = 0.0;
 
 		for (int k=0; k<Nz; k++){
 			for (int j=0; j<Ny; j++){
@@ -1359,14 +1353,35 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 		}	
 		volume_connected = sumReduce( Dm->Comm, count);
 		second_biggest = sumReduce( Dm->Comm, second_biggest);
+	}
+	else {
+		// use the whole NWP 
+		for (int k=0; k<Nz; k++){
+			for (int j=0; j<Ny; j++){
+				for (int i=0; i<Nx; i++){
+					if (Averages->SDs(i,j,k) > 0.f){
+						if (phase(i,j,k) > 0.f ){
+							phase_id(i,j,k) = 0;
+						}
+						else {
+							phase_id(i,j,k) = 1;
+						}
+					}
+					else {
+						phase_id(i,j,k) = 1;
+					}
+				}
+			}
+		}
+	}
 
-	int reach_x, reach_y, reach_z;
+	/*int reach_x, reach_y, reach_z;
 	for (int k=0; k<Nz; k++){
 		for (int j=0; j<Ny; j++){
 			for (int i=0; i<Nx; i++){
 			}
 		}
-	}
+	}*/
 
 	// 3. Generate a distance map to the largest object -> phase_distance
 	CalcDist(phase_distance,phase_id,*Dm);
@@ -1393,18 +1408,21 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 		}
 	}
 
+	if (USE_CONNECTED_NWP){
 	if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){
 		// if connected volume is less than 2% just delete the whole thing
 		if (rank==0) printf("Connected region has shrunk! \n");
 		REVERSE_FLOW_DIRECTION = true;
 	}
+	
 /*	else{*/
 		if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest );
+	}
 		if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
 		double target_delta_volume_incremental = target_delta_volume;
 		if (fabs(target_delta_volume) > 0.01*volume_initial)  
 			target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
-		delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental);
+		delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor);
 
 		for (int k=0; k<Nz; k++){
 			for (int j=0; j<Ny; j++){
@@ -1422,7 +1440,6 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 		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;
 					double d = phase_distance(i,j,k);
 					if (Averages->SDs(i,j,k) > 0.f){
 						if (d < 3.f){
@@ -1446,7 +1463,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 			}
 		}
 	}
-	volume_final= sumReduce( Dm->Comm, count);
+	double volume_final= sumReduce( Dm->Comm, count);
 
 	delta_volume = (volume_final-volume_initial);
 	if (rank == 0)  printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial);
@@ -1469,7 +1486,7 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
 	// 7. Re-initialize phase field and density
 	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
 	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
-	if (BoundaryCondition >0 ){
+	if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){
 		if (Dm->kproc()==0){
 			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
 			ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
@@ -1538,25 +1555,25 @@ void ScaLBL_ColorModel::WriteDebug(){
 	fwrite(PhaseField.data(),8,N,VELZ_FILE);
 	fclose(VELZ_FILE);
 
-//	ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
-//	FILE *CGX_FILE;
-//	sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank);
-//	CGX_FILE = fopen(LocalRankFilename,"wb");
-//	fwrite(PhaseField.data(),8,N,CGX_FILE);
-//	fclose(CGX_FILE);
-//
-//	ScaLBL_Comm->RegularLayout(Map,&ColorGrad[Np],PhaseField);
-//	FILE *CGY_FILE;
-//	sprintf(LocalRankFilename,"Gradient_Y.%05i.raw",rank);
-//	CGY_FILE = fopen(LocalRankFilename,"wb");
-//	fwrite(PhaseField.data(),8,N,CGY_FILE);
-//	fclose(CGY_FILE);
-//
-//	ScaLBL_Comm->RegularLayout(Map,&ColorGrad[2*Np],PhaseField);
-//	FILE *CGZ_FILE;
-//	sprintf(LocalRankFilename,"Gradient_Z.%05i.raw",rank);
-//	CGZ_FILE = fopen(LocalRankFilename,"wb");
-//	fwrite(PhaseField.data(),8,N,CGZ_FILE);
-//	fclose(CGZ_FILE);
+/*	ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField);
+	FILE *CGX_FILE;
+	sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank);
+	CGX_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,CGX_FILE);
+	fclose(CGX_FILE);
 
+	ScaLBL_Comm->RegularLayout(Map,&ColorGrad[Np],PhaseField);
+	FILE *CGY_FILE;
+	sprintf(LocalRankFilename,"Gradient_Y.%05i.raw",rank);
+	CGY_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,CGY_FILE);
+	fclose(CGY_FILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&ColorGrad[2*Np],PhaseField);
+	FILE *CGZ_FILE;
+	sprintf(LocalRankFilename,"Gradient_Z.%05i.raw",rank);
+	CGZ_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,CGZ_FILE);
+	fclose(CGZ_FILE);
+*/
 }

models/DFHModel.cpp

@@ -114,7 +114,6 @@ void ScaLBL_DFHModel::SetDomain(){
 }
 
 void ScaLBL_DFHModel::ReadInput(){
-	size_t readID;
 	//.......................................................................
 	if (rank == 0)    printf("Read input media... \n");
 	//.......................................................................

models/GreyscaleModel.cpp

@@ -0,0 +1,897 @@
+/*
+Greyscale lattice boltzmann model
+ */
+#include "models/GreyscaleModel.h"
+#include "analysis/distance.h"
+#include "analysis/morphology.h"
+#include <stdlib.h>
+#include <time.h>
+
+template<class TYPE>
+void DeleteArray( const TYPE *p )
+{
+    delete [] p;
+}
+
+ScaLBL_GreyscaleModel::ScaLBL_GreyscaleModel(int RANK, int NP, MPI_Comm COMM):
+rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),tau_eff(0),Den(0),Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),GreyPorosity(0),
+Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
+{
+	SignDist.resize(Nx,Ny,Nz);           
+    SignDist.fill(0);
+
+}
+ScaLBL_GreyscaleModel::~ScaLBL_GreyscaleModel(){
+
+}
+
+void ScaLBL_GreyscaleModel::ReadParams(string filename){
+	// read the input database 
+	db = std::make_shared<Database>( filename );
+	domain_db = db->getDatabase( "Domain" );
+	greyscale_db =  db->getDatabase( "Greyscale" );
+	analysis_db = db->getDatabase( "Analysis" );
+	vis_db = db->getDatabase( "Visualization" );
+
+	// set defaults
+	timestepMax = 100000;
+	tau = 1.0;
+    tau_eff = tau;
+    Den = 1.0;//constant density
+	tolerance = 0.01;
+	Fx = Fy = Fz = 0.0;
+	Restart=false;
+	din=dout=1.0;
+	flux=0.0;
+    dp = 10.0; //unit of 'dp': voxel
+    CollisionType = 1; //1: IMRT; 2: BGK
+	
+	// ---------------------- Greyscale Model parameters -----------------------//
+	if (greyscale_db->keyExists( "timestepMax" )){
+		timestepMax = greyscale_db->getScalar<int>( "timestepMax" );
+	}
+	if (greyscale_db->keyExists( "tau" )){
+		tau = greyscale_db->getScalar<double>( "tau" );
+	}
+	tau_eff = greyscale_db->getWithDefault<double>( "tau_eff", tau );
+	if (greyscale_db->keyExists( "Den" )){
+		Den = greyscale_db->getScalar<double>( "Den" );
+	}
+	if (greyscale_db->keyExists( "dp" )){
+		dp = greyscale_db->getScalar<double>( "dp" );
+	}
+	if (greyscale_db->keyExists( "F" )){
+		Fx = greyscale_db->getVector<double>( "F" )[0];
+		Fy = greyscale_db->getVector<double>( "F" )[1];
+		Fz = greyscale_db->getVector<double>( "F" )[2];
+	}
+	if (greyscale_db->keyExists( "Restart" )){
+		Restart = greyscale_db->getScalar<bool>( "Restart" );
+	}
+	if (greyscale_db->keyExists( "din" )){
+		din = greyscale_db->getScalar<double>( "din" );
+	}
+	if (greyscale_db->keyExists( "dout" )){
+		dout = greyscale_db->getScalar<double>( "dout" );
+	}
+	if (greyscale_db->keyExists( "flux" )){
+		flux = greyscale_db->getScalar<double>( "flux" );
+	}
+	if (greyscale_db->keyExists( "tolerance" )){
+		tolerance = greyscale_db->getScalar<double>( "tolerance" );
+	}
+	auto collision = greyscale_db->getWithDefault<std::string>( "collision", "IMRT" );
+	if (collision == "BGK"){
+        CollisionType=2;
+	}
+	// ------------------------------------------------------------------------//
+    
+    //------------------------ Other Domain parameters ------------------------//
+	BoundaryCondition = 0;
+	if (domain_db->keyExists( "BC" )){
+		BoundaryCondition = domain_db->getScalar<int>( "BC" );
+	}
+	// ------------------------------------------------------------------------//
+}
+
+void ScaLBL_GreyscaleModel::SetDomain(){
+	Dm  = std::shared_ptr<Domain>(new Domain(domain_db,comm));      // full domain for analysis
+	Mask  = std::shared_ptr<Domain>(new Domain(domain_db,comm));    // mask domain removes immobile phases
+	// domain parameters
+	Nx = Dm->Nx;
+	Ny = Dm->Ny;
+	Nz = Dm->Nz;
+	Lx = Dm->Lx;
+	Ly = Dm->Ly;
+	Lz = Dm->Lz;
+	N = Nx*Ny*Nz;
+
+	SignDist.resize(Nx,Ny,Nz);
+	Velocity_x.resize(Nx,Ny,Nz);
+	Velocity_y.resize(Nx,Ny,Nz);
+	Velocity_z.resize(Nx,Ny,Nz);
+	PorosityMap.resize(Nx,Ny,Nz);
+	Pressure.resize(Nx,Ny,Nz);
+
+	id = new signed char [N];
+	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1;               // initialize this way
+	MPI_Barrier(comm);
+	Dm->CommInit();
+	MPI_Barrier(comm);
+	// Read domain parameters
+	rank = Dm->rank();	
+	nprocx = Dm->nprocx();
+	nprocy = Dm->nprocy();
+	nprocz = Dm->nprocz();
+}
+
+void ScaLBL_GreyscaleModel::ReadInput(){
+	
+	sprintf(LocalRankString,"%05d",rank);
+	sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
+	sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
+	
+	if (domain_db->keyExists( "Filename" )){
+		auto Filename = domain_db->getScalar<std::string>( "Filename" );
+		Mask->Decomp(Filename);
+	}
+	else{
+        if (rank==0) printf("Filename of input image is not found, reading ID.0* instead.");
+		Mask->ReadIDs();
+	}
+	for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i];  // save what was read
+	
+	// Generate the signed distance map
+	// Initialize the domain and communication
+	Array<char> id_solid(Nx,Ny,Nz);
+	int count = 0;
+	// Solve for the position of the solid phase
+	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;
+				// Initialize the solid phase
+				signed char label = Mask->id[n];
+				if (label > 0)		id_solid(i,j,k) = 1;
+				else	     		id_solid(i,j,k) = 0;
+			}
+		}
+	}
+	// Initialize the signed distance function
+	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;
+				// Initialize distance to +/- 1
+				SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
+			}
+		}
+	}
+//	MeanFilter(SignDist);
+	if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
+	CalcDist(SignDist,id_solid,*Mask);
+	
+	if (rank == 0) cout << "Domain set." << endl;
+}
+
+/********************************************************
+ * AssignComponentLabels                                 *
+ ********************************************************/
+void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Permeability)
+{
+	size_t NLABELS=0;
+	signed char VALUE=0;
+	double POROSITY=0.f;
+	double PERMEABILITY=0.f;
+
+	auto LabelList = greyscale_db->getVector<int>( "ComponentLabels" );
+	auto PorosityList = greyscale_db->getVector<double>( "PorosityList" );
+	auto PermeabilityList = greyscale_db->getVector<double>( "PermeabilityList" );
+
+	NLABELS=LabelList.size();
+	if (NLABELS != PorosityList.size()){
+		ERROR("Error: ComponentLabels and PorosityList must be the same length! \n");
+	}
+
+	double label_count[NLABELS];
+	double label_count_global[NLABELS];
+	// Assign the labels
+
+	for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
+
+	for (int k=1;k<Nz-1;k++){
+		for (int j=1;j<Ny-1;j++){
+			for (int i=1;i<Nx-1;i++){
+				int n = k*Nx*Ny+j*Nx+i;
+				VALUE=id[n];
+				// Assign the affinity from the paired list
+				for (unsigned int idx=0; idx < NLABELS; idx++){
+				      //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
+					if (VALUE == LabelList[idx]){
+						POROSITY=PorosityList[idx];
+						label_count[idx] += 1.0;
+						idx = NLABELS;
+						//Mask->id[n] = 0; // set mask to zero since this is an immobile component
+					}
+				}
+				int idx = Map(i,j,k);
+				if (!(idx < 0)){
+                    if (POROSITY<=0.0){
+                        ERROR("Error: Porosity for grey voxels must be 0.0 < Porosity <= 1.0 !\n");
+                    }
+                    else{
+					    Porosity[idx] = POROSITY;
+                    }
+                }
+			}
+		}
+	}
+
+	if (NLABELS != PermeabilityList.size()){
+		ERROR("Error: ComponentLabels and PermeabilityList must be the same length! \n");
+	}
+	for (int k=1;k<Nz-1;k++){
+		for (int j=1;j<Ny-1;j++){
+			for (int i=1;i<Nx-1;i++){
+				int n = k*Nx*Ny+j*Nx+i;
+				VALUE=id[n];
+				// Assign the affinity from the paired list
+				for (unsigned int idx=0; idx < NLABELS; idx++){
+					//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
+					if (VALUE == LabelList[idx]){
+						PERMEABILITY=PermeabilityList[idx];
+						idx = NLABELS;
+						//Mask->id[n] = 0; // set mask to zero since this is an immobile component
+					}
+				}
+				int idx = Map(i,j,k);
+				if (!(idx < 0)){
+                    if (PERMEABILITY<=0.0){
+                        ERROR("Error: Permeability for grey voxel must be > 0.0 ! \n");
+                    }
+                    else{
+					    Permeability[idx] = PERMEABILITY/Dm->voxel_length/Dm->voxel_length;
+                    }
+                }
+			}
+		}
+	}
+
+
+	// Set Dm to match Mask
+	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i]; 
+	
+	for (int idx=0; idx<NLABELS; idx++)		label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
+    //Initialize a weighted porosity after considering grey voxels
+    GreyPorosity=0.0;
+	for (unsigned int idx=0; idx<NLABELS; idx++){
+		double volume_fraction  = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+        GreyPorosity+=volume_fraction*PorosityList[idx];
+    }
+
+	if (rank==0){
+        printf("Image resolution: %.5g [um/voxel]\n",Dm->voxel_length);
+		printf("Component labels: %lu \n",NLABELS);
+		for (unsigned int idx=0; idx<NLABELS; idx++){
+			VALUE=LabelList[idx];
+			POROSITY=PorosityList[idx];
+			PERMEABILITY=PermeabilityList[idx];
+			double volume_fraction  = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
+			printf("   label=%d: porosity=%.3g, permeability=%.3g [um^2] (=%.3g [voxel^2]), volume fraction=%.3g\n",
+                    VALUE,POROSITY,PERMEABILITY,PERMEABILITY/Dm->voxel_length/Dm->voxel_length,volume_fraction); 
+            printf("             effective porosity=%.3g\n",volume_fraction*POROSITY);
+		}
+        printf("The weighted porosity, considering both open and grey voxels, is %.3g\n",GreyPorosity);
+	}
+}
+
+
+void ScaLBL_GreyscaleModel::Create(){
+	/*
+	 *  This function creates the variables needed to run a LBM 
+	 */
+	//.........................................................
+	// don't perform computations at the eight corners
+	//id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
+	//id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
+
+	//.........................................................
+	// Initialize communication structures in averaging domain
+	for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
+	Mask->CommInit();
+	Np=Mask->PoreCount();
+	//...........................................................................
+	if (rank==0)    printf ("Create ScaLBL_Communicator \n");
+	// Create a communicator for the device (will use optimized layout)
+	// ScaLBL_Communicator ScaLBL_Comm(Mask); // original
+	ScaLBL_Comm  = std::shared_ptr<ScaLBL_Communicator>(new ScaLBL_Communicator(Mask));
+
+	int Npad=(Np/16 + 2)*16;
+	if (rank==0)    printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N);
+	Map.resize(Nx,Ny,Nz);       Map.fill(-2);
+	auto neighborList= new int[18*Npad];
+	Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
+	MPI_Barrier(comm);
+
+	//...........................................................................
+	//                MAIN  VARIABLES ALLOCATED HERE
+	//...........................................................................
+	// LBM variables
+	if (rank==0)    printf ("Allocating distributions \n");
+	//......................device distributions.................................
+	dist_mem_size = Np*sizeof(double);
+	neighborSize=18*(Np*sizeof(int));
+	//...........................................................................
+	ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
+	ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size);
+	ScaLBL_AllocateDeviceMemory((void **) &Permeability, sizeof(double)*Np);		
+	ScaLBL_AllocateDeviceMemory((void **) &Porosity, sizeof(double)*Np);		
+	ScaLBL_AllocateDeviceMemory((void **) &Pressure_dvc, sizeof(double)*Np);
+	ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
+	//...........................................................................
+	// Update GPU data structures
+	if (rank==0)	printf ("Setting up device map and neighbor list \n");
+	fflush(stdout);
+	int *TmpMap;
+	TmpMap=new int[Np];
+	for (int k=1; k<Nz-1; k++){
+		for (int j=1; j<Ny-1; j++){
+			for (int i=1; i<Nx-1; i++){
+				int idx=Map(i,j,k);
+				if (!(idx < 0))
+					TmpMap[idx] = k*Nx*Ny+j*Nx+i;
+			}
+		}
+	}
+	// check that TmpMap is valid
+	for (int idx=0; idx<ScaLBL_Comm->LastExterior(); idx++){
+		int n = TmpMap[idx];
+		if (n > Nx*Ny*Nz){
+			printf("Bad value! idx=%i \n");
+			TmpMap[idx] = Nx*Ny*Nz-1;
+		}
+	}
+	for (int idx=ScaLBL_Comm->FirstInterior(); idx<ScaLBL_Comm->LastInterior(); idx++){
+		int n = TmpMap[idx];
+		if (n > Nx*Ny*Nz){
+			printf("Bad value! idx=%i \n");
+			TmpMap[idx] = Nx*Ny*Nz-1;
+		}
+	}
+	ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
+	ScaLBL_DeviceBarrier();
+	delete [] TmpMap;
+	
+	// copy the neighbor list 
+	ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
+	// initialize phi based on PhaseLabel (include solid component labels)
+	double *Poros, *Perm;
+	Poros = new double[Np];
+	Perm = new double[Np];
+	AssignComponentLabels(Poros,Perm);
+	ScaLBL_CopyToDevice(Porosity, Poros, Np*sizeof(double));
+	ScaLBL_CopyToDevice(Permeability, Perm, Np*sizeof(double));
+}        
+
+
+void ScaLBL_GreyscaleModel::Initialize(){
+	if (rank==0)	printf ("Initializing distributions \n");
+    //TODO: for BGK, you need to consider voxel porosity
+    //      for IMRT, the whole set of feq is different
+    //      if in the future you have different collison mode, need to write two set of initialization functions
+    if (CollisionType==1){
+	    ScaLBL_D3Q19_GreyIMRT_Init(fq, Np, Den);
+        if (rank==0) printf("Collision model: Incompressible MRT.\n");
+    }
+    else if (CollisionType==2){
+	    ScaLBL_D3Q19_Init(fq, Np);
+        if (rank==0) printf("Collision model: BGK.\n");
+    }
+    else{
+        if (rank==0) printf("Unknown collison type! IMRT collision is used.\n"); 
+	    ScaLBL_D3Q19_GreyIMRT_Init(fq, Np, Den);
+        CollisionType=1;
+		greyscale_db->putScalar<std::string>( "collision", "IMRT" );
+    }
+
+	if (Restart == true){
+		if (rank==0){
+			printf("Initializing distributions from Restart! \n");
+		}
+		// Read in the restart file to CPU buffers
+        std::shared_ptr<double> cfq;
+        cfq = std::shared_ptr<double>(new double[19*Np],DeleteArray<double>);
+        FILE *File;
+        File=fopen(LocalRestartFile,"rb");
+        fread(cfq.get(),sizeof(double),19*Np,File);
+        fclose(File);
+
+		// Copy the restart data to the GPU
+		ScaLBL_CopyToDevice(fq,cfq.get(),19*Np*sizeof(double));
+		ScaLBL_DeviceBarrier();
+
+		MPI_Barrier(comm);
+	}
+}
+
+void ScaLBL_GreyscaleModel::Run(){
+	int nprocs=nprocx*nprocy*nprocz;
+	const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
+	
+	int analysis_interval = 1000; 	// number of timesteps in between in situ analysis 
+	int visualization_interval = 1000; 	 
+	int restart_interval = 10000; 	// number of timesteps in between in saving distributions for restart 
+	if (analysis_db->keyExists( "analysis_interval" )){
+		analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
+	}
+	if (analysis_db->keyExists( "visualization_interval" )){
+		visualization_interval = analysis_db->getScalar<int>( "visualization_interval" );
+	}
+	if (analysis_db->keyExists( "restart_interval" )){
+		restart_interval = analysis_db->getScalar<int>( "restart_interval" );
+	}
+	if (greyscale_db->keyExists( "timestep" )){
+		timestep = greyscale_db->getScalar<int>( "timestep" );
+	}
+
+	if (rank==0){
+		printf("********************************************************\n");
+		printf("No. of timesteps: %i \n", timestepMax);
+		fflush(stdout);
+	}
+
+	//.......create and start timer............
+	double starttime,stoptime,cputime;
+	ScaLBL_DeviceBarrier();
+	MPI_Barrier(comm);
+	starttime = MPI_Wtime();
+	//.........................................
+	
+	Minkowski Morphology(Mask);
+
+	//************ MAIN ITERATION LOOP ***************************************/
+	PROFILE_START("Loop");
+	auto current_db = db->cloneDatabase();
+	double rlx = 1.0/tau;
+    double rlx_eff = 1.0/tau_eff;
+	double error = 1.0;
+	double flow_rate_previous = 0.0;
+	while (timestep < timestepMax && error > tolerance) {
+		//************************************************************************/
+		// *************ODD TIMESTEP*************//
+		timestep++;
+		ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
+        switch (CollisionType){
+            case 1: 
+                    ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq,  ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+            case 2: 
+                    ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq,  ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
+                    break;
+            default: 
+                    ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq,  ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+        }
+		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		ScaLBL_DeviceBarrier();
+		// Set BCs
+		if (BoundaryCondition == 3){
+			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+        switch (CollisionType){
+            case 1: 
+		            ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+            case 2: 
+		            ScaLBL_D3Q19_AAodd_Greyscale(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
+                    break;
+            default: 
+		            ScaLBL_D3Q19_AAodd_Greyscale_IMRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+        }
+		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+
+		// *************EVEN TIMESTEP*************//
+		timestep++;
+		ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
+        switch (CollisionType){
+            case 1: 
+		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+            case 2: 
+		            ScaLBL_D3Q19_AAeven_Greyscale(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
+                    break;
+            default: 
+		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+        }
+        ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		ScaLBL_DeviceBarrier();
+		// Set BCs
+		if (BoundaryCondition == 3){
+			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+        switch (CollisionType){
+            case 1: 
+		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+            case 2: 
+		            ScaLBL_D3Q19_AAeven_Greyscale(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Pressure_dvc);
+                    break;
+            default: 
+		            ScaLBL_D3Q19_AAeven_Greyscale_IMRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx, rlx_eff, Fx, Fy, Fz,Porosity,Permeability,Velocity,Den,Pressure_dvc);
+                    break;
+        }
+        ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+		//************************************************************************/
+		
+		if (timestep%analysis_interval==0){
+			ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
+			ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
+			ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
+			//ScaLBL_Comm->RegularLayout(Map,Porosity,PorosityMap);
+			//ScaLBL_Comm->RegularLayout(Map,Pressure_dvc,Pressure);
+			
+			double count_loc=0;
+			double count;
+			double vax,vay,vaz;
+			double vax_loc,vay_loc,vaz_loc;
+            //double px_loc,py_loc,pz_loc;
+            //double px,py,pz;
+            //double mass_loc,mass_glb;
+            
+            //parameters for domain average
+	        int64_t i,j,k,n,imin,jmin,kmin,kmax;
+            // If external boundary conditions are set, do not average over the inlet and outlet
+            kmin=1; kmax=Nz-1;
+            //In case user forgets to specify the inlet/outlet buffer layers for BC>0
+            if (BoundaryCondition > 0 && Dm->kproc() == 0) kmin=4;
+            if (BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
+
+            imin=jmin=1;
+            // If inlet/outlet layers exist use these as default
+            //if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
+            //if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
+            if (BoundaryCondition > 0 && Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin = 1 + Dm->inlet_layers_z;//"1" indicates the halo layer
+            if (BoundaryCondition > 0 && Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax = Nz-1 - Dm->outlet_layers_z; 
+
+//			px_loc = py_loc = pz_loc = 0.f;
+//            mass_loc = 0.f;
+//			for (int k=kmin; k<kmax; k++){
+//				for (int j=jmin; j<Ny-1; j++){
+//					for (int i=imin; i<Nx-1; i++){
+//						if (SignDist(i,j,k) > 0){
+//							px_loc   += Velocity_x(i,j,k)*Den*PorosityMap(i,j,k);
+//							py_loc   += Velocity_y(i,j,k)*Den*PorosityMap(i,j,k);
+//							pz_loc   += Velocity_z(i,j,k)*Den*PorosityMap(i,j,k);
+//							mass_loc += Den*PorosityMap(i,j,k);
+//						}
+//					}
+//				}
+//			}
+//			MPI_Allreduce(&px_loc,  &px,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			MPI_Allreduce(&py_loc,  &py,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			MPI_Allreduce(&pz_loc,  &pz,      1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			MPI_Allreduce(&mass_loc,&mass_glb,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+//			
+//			vax = px/mass_glb;
+//			vay = py/mass_glb;
+//			vaz = pz/mass_glb;
+            
+			vax_loc = vay_loc = vaz_loc = 0.f;
+			for (int k=kmin; k<kmax; k++){
+				for (int j=jmin; j<Ny-1; j++){
+					for (int i=imin; i<Nx-1; i++){
+						if (SignDist(i,j,k) > 0){
+							vax_loc += Velocity_x(i,j,k);
+							vay_loc += Velocity_y(i,j,k);
+							vaz_loc += Velocity_z(i,j,k);
+							count_loc+=1.0;
+						}
+					}
+				}
+			}
+			//MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+			//MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+			//MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+			//MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+			
+            vax = sumReduce( Mask->Comm, vax_loc);
+            vay = sumReduce( Mask->Comm, vay_loc);
+            vaz = sumReduce( Mask->Comm, vaz_loc);
+            count = sumReduce( Mask->Comm, count_loc);
+
+			vax /= count;
+			vay /= count;
+			vaz /= count;
+
+			double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
+			double dir_x = Fx/force_mag;
+			double dir_y = Fy/force_mag;
+			double dir_z = Fz/force_mag;
+			if (force_mag == 0.0){
+				// default to z direction
+				dir_x = 0.0;
+				dir_y = 0.0;
+				dir_z = 1.0;
+				force_mag = 1.0;
+			}
+			//double flow_rate = (px*dir_x + py*dir_y + pz*dir_z)/mass_glb;
+			double flow_rate = (vax*dir_x + vay*dir_y + vaz*dir_z);
+			
+			error = fabs(flow_rate - flow_rate_previous) / fabs(flow_rate);
+			flow_rate_previous = flow_rate;
+			
+			//if (rank==0) printf("Computing Minkowski functionals \n");
+			Morphology.ComputeScalar(SignDist,0.f);
+			//Morphology.PrintAll();
+			double mu = (tau-0.5)/3.f;
+			double Vs = Morphology.V();
+			double As = Morphology.A();
+			double Hs = Morphology.H();
+			double Xs = Morphology.X();
+			Vs = sumReduce( Dm->Comm, Vs);
+			As = sumReduce( Dm->Comm, As);
+			Hs = sumReduce( Dm->Comm, Hs);
+			Xs = sumReduce( Dm->Comm, Xs);
+
+			double h = Dm->voxel_length;
+			//double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
+			double absperm = h*h*mu*GreyPorosity*flow_rate / force_mag;
+
+            if (rank==0){
+				printf("     AbsPerm = %.5g [micron^2]\n",absperm);
+                bool WriteHeader=false;
+                FILE * log_file = fopen("Permeability.csv","r");
+                if (log_file != NULL)
+                    fclose(log_file);
+                else
+                    WriteHeader=true;
+                log_file = fopen("Permeability.csv","a");
+                if (WriteHeader)
+                    fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n",
+                            timestep,Fx,Fy,Fz,mu,h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz,absperm);
+
+                fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, 
+                        h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm);
+                fclose(log_file);
+            }
+		}
+
+		if (timestep%visualization_interval==0){
+            VelocityField();
+        }
+
+		if (timestep%restart_interval==0){
+            //Use rank=0 write out Restart.db
+            if (rank==0) {
+                greyscale_db->putScalar<int>("timestep",timestep);    		
+                greyscale_db->putScalar<bool>( "Restart", true );
+                current_db->putDatabase("Greyscale", greyscale_db);
+                std::ofstream OutStream("Restart.db");
+                current_db->print(OutStream, "");
+                OutStream.close();
+      
+            }
+            //Write out Restart data.
+            std::shared_ptr<double> cfq;
+            cfq = std::shared_ptr<double>(new double[19*Np],DeleteArray<double>);
+            ScaLBL_CopyToHost(cfq.get(),fq,19*Np*sizeof(double));// Copy restart data to the CPU
+
+            FILE *RESTARTFILE;
+            RESTARTFILE=fopen(LocalRestartFile,"wb");
+            fwrite(cfq.get(),sizeof(double),19*Np,RESTARTFILE);
+            fclose(RESTARTFILE);
+		    MPI_Barrier(comm);
+        }
+	}
+
+	PROFILE_STOP("Loop");
+	PROFILE_SAVE("lbpm_greyscale_simulator",1);
+	//************************************************************************
+	ScaLBL_DeviceBarrier();
+	MPI_Barrier(comm);
+	stoptime = MPI_Wtime();
+	if (rank==0) printf("-------------------------------------------------------------------\n");
+	// Compute the walltime per timestep
+	cputime = (stoptime - starttime)/timestep;
+	// Performance obtained from each node
+	double MLUPS = double(Np)/cputime/1000000;
+
+	if (rank==0) printf("********************************************************\n");
+	if (rank==0) printf("CPU time = %f \n", cputime);
+	if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
+	MLUPS *= nprocs;
+	if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
+	if (rank==0) printf("********************************************************\n");
+
+	// ************************************************************************
+}
+
+void ScaLBL_GreyscaleModel::VelocityField(){
+
+/*	Minkowski Morphology(Mask);
+	int SIZE=Np*sizeof(double);
+	ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
+	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
+	ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
+
+	memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));
+	Morphology.Initialize();
+	Morphology.UpdateMeshValues();
+	Morphology.ComputeLocal();
+	Morphology.Reduce();
+	
+	double count_loc=0;
+	double count;
+	double vax,vay,vaz;
+	double vax_loc,vay_loc,vaz_loc;
+	vax_loc = vay_loc = vaz_loc = 0.f;
+	for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
+		vax_loc += VELOCITY[n];
+		vay_loc += VELOCITY[Np+n];
+		vaz_loc += VELOCITY[2*Np+n];
+		count_loc+=1.0;
+	}
+	
+	for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
+		vax_loc += VELOCITY[n];
+		vay_loc += VELOCITY[Np+n];
+		vaz_loc += VELOCITY[2*Np+n];
+		count_loc+=1.0;
+	}
+	MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	
+	vax /= count;
+	vay /= count;
+	vaz /= count;
+	
+	double mu = (tau-0.5)/3.f;
+	if (rank==0) printf("Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
+	if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Fx, Fy, Fz, mu, 
+						Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
+						*/
+	
+	std::vector<IO::MeshDataStruct> visData;
+	fillHalo<double> fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);
+
+	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 PressureVar = std::make_shared<IO::Variable>();
+
+	IO::initialize("","silo","false");
+	// Create the MeshDataStruct	
+	visData.resize(1);
+	visData[0].meshName = "domain";
+	visData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
+	SignDistVar->name = "SignDist";
+	SignDistVar->type = IO::VariableType::VolumeVariable;
+	SignDistVar->dim = 1;
+	SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	visData[0].vars.push_back(SignDistVar);
+	
+	VxVar->name = "Velocity_x";
+	VxVar->type = IO::VariableType::VolumeVariable;
+	VxVar->dim = 1;
+	VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	visData[0].vars.push_back(VxVar);
+	VyVar->name = "Velocity_y";
+	VyVar->type = IO::VariableType::VolumeVariable;
+	VyVar->dim = 1;
+	VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	visData[0].vars.push_back(VyVar);
+	VzVar->name = "Velocity_z";
+	VzVar->type = IO::VariableType::VolumeVariable;
+	VzVar->dim = 1;
+	VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	visData[0].vars.push_back(VzVar);
+	
+	PressureVar->name = "Pressure";
+	PressureVar->type = IO::VariableType::VolumeVariable;
+	PressureVar->dim = 1;
+	PressureVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	visData[0].vars.push_back(PressureVar);
+
+	Array<double>& SignData  = visData[0].vars[0]->data;
+	Array<double>& VelxData = visData[0].vars[1]->data;
+	Array<double>& VelyData = visData[0].vars[2]->data;
+	Array<double>& VelzData = visData[0].vars[3]->data;
+	Array<double>& PressureData = visData[0].vars[4]->data;
+	
+    ASSERT(visData[0].vars[0]->name=="SignDist");
+    ASSERT(visData[0].vars[1]->name=="Velocity_x");
+    ASSERT(visData[0].vars[2]->name=="Velocity_y");
+    ASSERT(visData[0].vars[3]->name=="Velocity_z");
+    ASSERT(visData[0].vars[4]->name=="Pressure");
+	
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
+	ScaLBL_Comm->RegularLayout(Map,Pressure_dvc,Pressure);
+
+    fillData.copy(SignDist,SignData);
+    fillData.copy(Velocity_x,VelxData);
+    fillData.copy(Velocity_y,VelyData);
+    fillData.copy(Velocity_z,VelzData);
+    fillData.copy(Pressure,PressureData);
+	
+    IO::writeData( timestep, visData, Dm->Comm );
+
+}
+
+void ScaLBL_GreyscaleModel::WriteDebug(){
+	// Copy back final phase indicator field and convert to regular layout
+	DoubleArray PhaseField(Nx,Ny,Nz);
+
+	//ScaLBL_CopyToHost(Porosity.data(), Poros, sizeof(double)*N);
+
+//	FILE *OUTFILE;
+//	sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
+//	OUTFILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,OUTFILE);
+//	fclose(OUTFILE);
+//
+//    ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
+//	FILE *AFILE;
+//	sprintf(LocalRankFilename,"A.%05i.raw",rank);
+//	AFILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,AFILE);
+//	fclose(AFILE);
+//
+//	ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
+//	FILE *BFILE;
+//	sprintf(LocalRankFilename,"B.%05i.raw",rank);
+//	BFILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,BFILE);
+//	fclose(BFILE);
+//
+//	ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField);
+//	FILE *PFILE;
+//	sprintf(LocalRankFilename,"Pressure.%05i.raw",rank);
+//	PFILE = fopen(LocalRankFilename,"wb");
+//	fwrite(PhaseField.data(),8,N,PFILE);
+//	fclose(PFILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
+	FILE *VELX_FILE;
+	sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank);
+	VELX_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,VELX_FILE);
+	fclose(VELX_FILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField);
+	FILE *VELY_FILE;
+	sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank);
+	VELY_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,VELY_FILE);
+	fclose(VELY_FILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField);
+	FILE *VELZ_FILE;
+	sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank);
+	VELZ_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,VELZ_FILE);
+	fclose(VELZ_FILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Porosity[0],PhaseField);
+	FILE *POROS_FILE;
+	sprintf(LocalRankFilename,"Porosity.%05i.raw",rank);
+	POROS_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,POROS_FILE);
+	fclose(POROS_FILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Permeability[0],PhaseField);
+	FILE *PERM_FILE;
+	sprintf(LocalRankFilename,"Permeability.%05i.raw",rank);
+	PERM_FILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,PERM_FILE);
+	fclose(PERM_FILE);
+}

models/GreyscaleModel.h

@@ -0,0 +1,92 @@
+/*
+Implementation of color lattice boltzmann model
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/stat.h>
+#include <iostream>
+#include <exception>
+#include <stdexcept>
+#include <fstream>
+
+#include "common/Communication.h"
+#include "common/MPI_Helpers.h"
+#include "common/Database.h"
+#include "common/ScaLBL.h"
+#include "ProfilerApp.h"
+#include "threadpool/thread_pool.h"
+
+class ScaLBL_GreyscaleModel{
+public:
+	ScaLBL_GreyscaleModel(int RANK, int NP, MPI_Comm COMM);
+	~ScaLBL_GreyscaleModel();	
+	
+	// functions in they should be run
+	void ReadParams(string filename);
+	void ReadParams(std::shared_ptr<Database> db0);
+	void SetDomain();
+	void ReadInput();
+	void Create();
+	void Initialize();
+	void Run();
+	void WriteDebug();
+	void VelocityField();
+	
+	bool Restart,pBC;
+	int timestep,timestepMax;
+	int BoundaryCondition;
+    int CollisionType;
+	double tau;
+    double tau_eff;
+    double Den;//constant density
+	double tolerance;
+	double Fx,Fy,Fz,flux;
+	double din,dout;
+    double dp;//solid particle diameter, unit in voxel
+    double GreyPorosity;
+	
+	int Nx,Ny,Nz,N,Np;
+	int rank,nprocx,nprocy,nprocz,nprocs;
+	double Lx,Ly,Lz;
+
+	std::shared_ptr<Domain> Dm;   // this domain is for analysis
+	std::shared_ptr<Domain> Mask; // this domain is for lbm
+	std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
+    
+    // input database
+    std::shared_ptr<Database> db;
+    std::shared_ptr<Database> domain_db;
+    std::shared_ptr<Database> greyscale_db;
+    std::shared_ptr<Database> analysis_db;
+    std::shared_ptr<Database> vis_db;
+
+    signed char *id;    
+	int *NeighborList;
+	int *dvcMap;
+	double *fq;
+	double *Permeability;//grey voxel permeability
+	double *Porosity;
+	double *Velocity;
+	double *Pressure_dvc;
+    IntArray Map;
+    DoubleArray SignDist;
+    DoubleArray Velocity_x;
+    DoubleArray Velocity_y;
+    DoubleArray Velocity_z;
+    DoubleArray PorosityMap;
+    DoubleArray Pressure;
+		
+private:
+	MPI_Comm comm;
+    
+	int dist_mem_size;
+	int neighborSize;
+	// filenames
+    char LocalRankString[8];
+    char LocalRankFilename[40];
+    char LocalRestartFile[40];
+   
+    void AssignComponentLabels(double *Porosity, double *Permeablity);
+    
+};
+

models/MRTModel.cpp

@@ -18,6 +18,7 @@
  */
 #include "models/MRTModel.h"
 #include "analysis/distance.h"
+#include "common/ReadMicroCT.h"
 
 ScaLBL_MRTModel::ScaLBL_MRTModel(int RANK, int NP, MPI_Comm COMM):
 rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tau(0),
@@ -108,20 +109,36 @@ void ScaLBL_MRTModel::SetDomain(){
 }
 
 void ScaLBL_MRTModel::ReadInput(){
-    int rank=Dm->rank();
-    size_t readID;
-    //.......................................................................
-    //.......................................................................
-    Mask->ReadIDs();
     
     sprintf(LocalRankString,"%05d",Dm->rank());
     sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
     sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
 
+    
+    if (domain_db->keyExists( "Filename" )){
+    	auto Filename = domain_db->getScalar<std::string>( "Filename" );
+    	Mask->Decomp(Filename);
+    }
+    else if (domain_db->keyExists( "GridFile" )){
+    	// Read the local domain data
+    	auto input_id = readMicroCT( *domain_db, comm );
+    	// Fill the halo (assuming GCW of 1)
+    	array<int,3> size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) };
+    	ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz };
+    	ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
+    	fillHalo<signed char> fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
+    	Array<signed char> id_view;
+    	id_view.viewRaw( size1, Mask->id );
+    	fill.copy( input_id, id_view );
+    	fill.fill( id_view );
+    }
+    else{
+    	Mask->ReadIDs();
+    }
+
     // Generate the signed distance map
 	// Initialize the domain and communication
 	Array<char> id_solid(Nx,Ny,Nz);
-	int count = 0;
 	// Solve for the position of the solid phase
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
@@ -137,7 +154,6 @@ void ScaLBL_MRTModel::ReadInput(){
 	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;
 				// Initialize distance to +/- 1
 				Distance(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
 			}
@@ -206,7 +222,6 @@ void ScaLBL_MRTModel::Run(){
 	double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
 	
 	Minkowski Morphology(Mask);
-	int SIZE=Np*sizeof(double);
 
 	if (rank==0){
 		bool WriteHeader=false;
@@ -238,12 +253,38 @@ void ScaLBL_MRTModel::Run(){
 		ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
 		ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq,  ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
 		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		// Set boundary conditions
+		if (BoundaryCondition == 3){
+			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		else if (BoundaryCondition == 4){
+			din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		else if (BoundaryCondition == 5){
+			ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
+			ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
+		}
 		ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
 		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 		timestep++;
 		ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
 		ScaLBL_D3Q19_AAeven_MRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
 		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
+		// Set boundary conditions
+		if (BoundaryCondition == 3){
+			ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		else if (BoundaryCondition == 4){
+			din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
+			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
+		}
+		else if (BoundaryCondition == 5){
+			ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
+			ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
+		}
 		ScaLBL_D3Q19_AAeven_MRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
 		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 		//************************************************************************/

tests/CMakeLists.txt

@@ -3,6 +3,7 @@
 #ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
+ADD_LBPM_EXECUTABLE( lbpm_greyscale_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator )

tests/GenerateSphereTest.cpp

@@ -145,8 +145,8 @@ inline void MorphOpen(DoubleArray SignDist, char *id, Domain &Dm, int nx, int ny
 
 	// Increase the critical radius until the target saturation is met
 	double deltaR=0.05; // amount to change the radius in voxel units
-	double Rcrit_old;
-	double Rcrit_new;
+	double Rcrit_old=0.0;
+	double Rcrit_new=0.0;
 
 	double GlobalNumber = 1.f;
 	int imin,jmin,kmin,imax,jmax,kmax;
@@ -363,7 +363,7 @@ int main(int argc, char **argv)
 	        nspheres = domain_db->getScalar<int>( "nspheres");
 
 		//printf("Set domain \n");
-		int BoundaryCondition=1;
+		//int BoundaryCondition=1;
 		//Nz += 2;
 		//Nx = Ny = Nz;	// Cubic domain
 		int N = Nx*Ny*Nz;
@@ -396,7 +396,7 @@ int main(int argc, char **argv)
 		int sum = 0;
 		double sum_local;
 		double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
-		double porosity, pore_vol;
+		double porosity;
 		//...........................................................................
 		DoubleArray SignDist(Nx,Ny,Nz);
 		//.......................................................................
@@ -450,7 +450,6 @@ int main(int argc, char **argv)
 			}
 		}
 		sum=0;
-		pore_vol = 0.0;
 		for ( k=1;k<Nz-1;k++){
 			for ( j=1;j<Ny-1;j++){
 				for ( i=1;i<Nx-1;i++){

tests/TestBubbleDFH.cpp

@@ -74,11 +74,6 @@ int main(int argc, char **argv)
         ScaLBL_DeviceBarrier();
         MPI_Barrier(comm);
 
-        PROFILE_ENABLE(1);
-        //PROFILE_ENABLE_TRACE();
-        //PROFILE_ENABLE_MEMORY();
-        PROFILE_SYNCHRONIZE();
-        PROFILE_START("Main");
         Utilities::setErrorHandlers();
 
         // Variables that specify the computational domain  
@@ -164,14 +159,14 @@ int main(int argc, char **argv)
             pBC=false;
 
         // Full domain used for averaging (do not use mask for analysis)
-        std::shared_ptr<Domain> Dm(new Domain(domain_db,comm));
+        auto Dm = std::make_shared<Domain>(domain_db,comm);
         for (int i=0; i<Dm->Nx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = 1;
-        std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
+        auto Averages = std::make_shared<TwoPhase>(Dm);
         //   TwoPhase Averages(Dm);
         Dm->CommInit();
 
         // Mask that excludes the solid phase
-        std::shared_ptr<Domain> Mask(new Domain(domain_db,comm));
+        auto Mask = std::make_shared<Domain>(domain_db,comm);
         MPI_Barrier(comm);
 
         Nx+=2; Ny+=2; Nz += 2;
@@ -191,9 +186,8 @@ int main(int argc, char **argv)
 		//	printf("Local File Name =  %s \n",LocalRankFilename);
 		// .......... READ THE INPUT FILE .......................................
 		//	char value;
-		char *id;
-		id = new char[N];
-		double sum, sum_local;
+		auto id = new char[N];
+		double sum;
 		//...........................................................................
 		if (rank == 0) cout << "Setting up bubble..." << endl;
 	    double BubbleRadius = 15.5; // Radius of the capillary tube
@@ -244,19 +238,17 @@ int main(int argc, char **argv)
 		// Initialize communication structures in averaging domain
 		for (i=0; i<Mask->Nx*Mask->Ny*Mask->Nz; i++) Mask->id[i] = id[i];
 		Mask->CommInit();
-		double *PhaseLabel;
-		PhaseLabel = new double[N];
+		auto PhaseLabel = new double[N];
 		
 		//...........................................................................
 		if (rank==0)	printf ("Create ScaLBL_Communicator \n");
 		// Create a communicator for the device (will use optimized layout)
-		std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm(new ScaLBL_Communicator(Mask));
+		auto ScaLBL_Comm = std::make_shared<ScaLBL_Communicator>(Mask);
 		
 		int Npad=(Np/16 + 2)*16;
 		if (rank==0)	printf ("Set up memory efficient layout Npad=%i \n",Npad);
-		int *neighborList;
 		IntArray Map(Nx,Ny,Nz);
-		neighborList= new int[18*Npad];
+		auto neighborList= new int[18*Npad];
 		Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
 		MPI_Barrier(comm);
 
@@ -515,9 +507,8 @@ int main(int argc, char **argv)
 		// Copy back final phase indicator field and convert to regular layout
 		DoubleArray PhaseField(Nx,Ny,Nz);
         ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField);
-    	FILE *OUTFILE;
-		sprintf(LocalRankFilename,"Phase.raw",rank);
-		OUTFILE = fopen(LocalRankFilename,"wb");
+		sprintf(LocalRankFilename,"Phase.raw");
+		auto OUTFILE = fopen(LocalRankFilename,"wb");
     	fwrite(PhaseField.data(),8,N,OUTFILE);
     	fclose(OUTFILE);
     	
@@ -535,9 +526,8 @@ int main(int argc, char **argv)
 				}
 			}
 		}
-    	FILE *GFILE;
 		sprintf(LocalRankFilename,"Gradient.raw");
-		GFILE = fopen(LocalRankFilename,"wb");
+		auto GFILE = fopen(LocalRankFilename,"wb");
     	fwrite(GradNorm.data(),8,N,GFILE);
     	fclose(GFILE);
     	
@@ -545,14 +535,12 @@ int main(int argc, char **argv)
 		DoubleArray Rho2(Nx,Ny,Nz);
         ScaLBL_Comm->RegularLayout(Map,&Den[0],Rho1);
         ScaLBL_Comm->RegularLayout(Map,&Den[Np],Rho2);
-    	FILE *RFILE1;
 		sprintf(LocalRankFilename,"Rho1.raw");
-		RFILE1 = fopen(LocalRankFilename,"wb");
+		auto RFILE1 = fopen(LocalRankFilename,"wb");
     	fwrite(Rho1.data(),8,N,RFILE1);
     	fclose(RFILE1);
-    	FILE *RFILE2;
 		sprintf(LocalRankFilename,"Rho2.raw");
-		RFILE2 = fopen(LocalRankFilename,"wb");
+		auto RFILE2 = fopen(LocalRankFilename,"wb");
     	fwrite(Rho2.data(),8,N,RFILE2);
     	fclose(RFILE2);
     	

tests/TestColorGradDFH.cpp

@@ -53,9 +53,6 @@ int main(int argc, char **argv)
 	    int Nx = db->getVector<int>( "n" )[0];
 	    int Ny = db->getVector<int>( "n" )[1];
 	    int Nz = db->getVector<int>( "n" )[2];
-	    int nprocx = db->getVector<int>( "nproc" )[0];
-	    int nprocy = db->getVector<int>( "nproc" )[1];
-	    int nprocz = db->getVector<int>( "nproc" )[2];
 
 	    if (rank==0){
 	    	printf("********************************************************\n");
@@ -64,7 +61,7 @@ int main(int argc, char **argv)
 	    }
 
 	    // Get the rank info
-	    std::shared_ptr<Domain> Dm(new Domain(db,comm));
+		auto Dm = std::make_shared<Domain>(db,comm);
 		Nx += 2;
 		Ny += 2;
 		Nz += 2;
@@ -111,7 +108,6 @@ int main(int argc, char **argv)
 		MPI_Barrier(comm);
 
 		//......................device distributions.................................
-		int dist_mem_size = Np*sizeof(double);
 		int neighborSize=18*Np*sizeof(int);
 		if (rank==0)	printf ("Allocating distributions \n");
 		int *NeighborList;

tests/TestCommD3Q19.cpp

@@ -49,7 +49,7 @@ extern void GlobalFlipScaLBL_D3Q19_Init(double *dist, IntArray Map, int Np, int
 	{1,1,0},{-1,-1,0},{1,-1,0},{-1,1,0},{1,0,1},{-1,0,-1},{1,0,-1},{-1,0,1},
 	{0,1,1},{0,-1,-1},{0,1,-1},{0,-1,1}};
 	
-	int q,i,j,k,n,N;
+	int q,i,j,k;
 	int Cqx,Cqy,Cqz; // Discrete velocity
 	int x,y,z;		// Global indices
 	int xn,yn,zn; 	// Global indices of neighbor 
@@ -59,8 +59,6 @@ extern void GlobalFlipScaLBL_D3Q19_Init(double *dist, IntArray Map, int Np, int
 	Y = Ny*nprocy;
 	Z = Nz*nprocz;
     NULL_USE(Z);
-	N = (Nx+2)*(Ny+2)*(Nz+2);	// size of the array including halo
-
 
 	for (k=0; k<Nz; k++){ 
 		for (j=0; j<Ny; j++){
@@ -104,16 +102,13 @@ extern int GlobalCheckDebugDist(double *dist, IntArray Map, int Np, int Nx, int
 {
 
 	int returnValue = 0;
-	int q,i,j,k,n,N,idx;
-	int Cqx,Cqy,Cqz; // Discrete velocity
+	int q,i,j,k,idx;
 	int x,y,z;		// Global indices
-	int xn,yn,zn; 	// Global indices of neighbor 
 	int X,Y,Z;		// Global size
 	X = Nx*nprocx;
 	Y = Ny*nprocy;
 	Z = Nz*nprocz;
 	NULL_USE(Z);
-	N = (Nx+2)*(Ny+2)*(Nz+2);	// size of the array including halo
 	for (k=0; k<Nz; k++){ 
 		for (j=0; j<Ny; j++){
 			for (i=0; i<Nx; i++){
@@ -168,9 +163,6 @@ inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
 //***************************************************************************************
 int main(int argc, char **argv)
 {
-	//*****************************************
-	// ***** MPI STUFF ****************
-	//*****************************************
 	// Initialize MPI
 	int rank,nprocs;
 	MPI_Init(&argc,&argv);
@@ -178,10 +170,7 @@ int main(int argc, char **argv)
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
 	int check;
-
 	{
-		MPI_Request req1[18],req2[18];
-		MPI_Status stat1[18],stat2[18];
 
 		if (rank == 0){
 			printf("********************************************************\n");
@@ -191,11 +180,8 @@ int main(int argc, char **argv)
 
 		// BGK Model parameters
 		string FILENAME;
-		unsigned int nBlocks, nthreads;
-		int timestepMax, interval;
-		double tau,Fx,Fy,Fz,tol;
 		// Domain variables
-		int i,j,k,n;
+		int i,j,k;
 
         // Load inputs
         auto db = loadInputs( nprocs );
@@ -223,8 +209,7 @@ int main(int argc, char **argv)
 		char LocalRankFilename[40];
 		sprintf(LocalRankFilename,"ID.%05i",rank);
 
-		char *id;
-		id = new char[Nx*Ny*Nz];
+		auto id = new char[Nx*Ny*Nz];
 
 		/*		if (rank==0) printf("Assigning phase ID from file \n");
 		if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
@@ -237,7 +222,7 @@ int main(int argc, char **argv)
 		for (k=0;k<Nz;k++){
 			for (j=0;j<Ny;j++){
 				for (i=0;i<Nx;i++){
-					n = k*Nx*Ny+j*Nx+i;
+					int n = k*Nx*Ny+j*Nx+i;
 					id[n] = 1;
 					Dm->id[n] = id[n];
 				}
@@ -270,7 +255,7 @@ int main(int argc, char **argv)
 		for (k=1;k<Nz-1;k++){
 			for (j=1;j<Ny-1;j++){
 				for (i=1;i<Nx-1;i++){
-					n = k*Nx*Ny+j*Nx+i;
+					int n = k*Nx*Ny+j*Nx+i;
 					if (id[n] == component){
 						sum_local+=1.0;
 					}

tests/TestFluxBC.cpp

@@ -32,21 +32,14 @@ int main (int argc, char **argv)
 	}
 	{
 	  int i,j,k,n,Np;
-		bool pBC=true;
-		double Lx,Ly,Lz;
-		Lx = Ly = Lz = 1.f;
 		double din,dout;
-		int BC=1;
 
 	    // Load inputs
 	    auto db = loadInputs( nprocs );
 	    int Nx = db->getVector<int>( "n" )[0];
 	    int Ny = db->getVector<int>( "n" )[1];
 	    int Nz = db->getVector<int>( "n" )[2];
-	    int nprocx = db->getVector<int>( "nproc" )[0];
-	    int nprocy = db->getVector<int>( "nproc" )[1];
-	    int nprocz = db->getVector<int>( "nproc" )[2];
-		std::shared_ptr<Domain> Dm(new Domain(db,comm));
+		auto Dm = std::make_shared<Domain>(db,comm);
 		
 		Nx += 2;   Ny+=2;	Nz += 2;
 		Nx = Ny = Nz;	// Cubic domain
@@ -55,8 +48,7 @@ int main (int argc, char **argv)
 		//.......................................................................
 		// Assign the phase ID
 		//.......................................................................
-		char *id;
-		id = new char[N];
+		auto id = new char[N];
 		for (k=0;k<Nz;k++){
 			for (j=0;j<Ny;j++){
 				for (i=0;i<Nx;i++){
@@ -160,9 +152,7 @@ int main (int argc, char **argv)
     	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
     	ScaLBL_CopyToHost(&VEL[0],&dvc_vel[0],SIZE);
 
-		double err,value,Q;
-
-    	Q = 0.f;    	
+    	double Q = 0.f;    	
     	k=1;
     	for (j=1;j<Ny-1;j++){
     		for (i=1;i<Nx-1;i++){
@@ -176,7 +166,7 @@ int main (int argc, char **argv)
 
     	// respect backwards read / write!!!
 		printf("Inlet Flux: input=%f, output=%f \n",flux,Q);
-		err = fabs(flux + Q);
+		double err = fabs(flux + Q);
 		if (err > 1e-12){
 			error = 1;
 			printf("  Inlet error %f \n",err);
@@ -185,7 +175,7 @@ int main (int argc, char **argv)
 		// Consider a larger number of timesteps and simulate flow
 		double Fx, Fy, Fz;
 		double tau = 1.0;
-		double mu=(tau-0.5)/3.0;
+		//double mu=(tau-0.5)/3.0;
 		double rlx_setA=1.0/tau;
 		double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
 		dout=1.f;

tests/TestForceD3Q19.cpp

@@ -458,23 +458,15 @@ int main (int argc, char **argv)
 		ASSERT(x!=NULL);
 	}
     
-    // set the error code
-    // Note: the error code should be consistent across all processors
-    int error = 0;
-    
     int Np = 1;
-    int Q = 9;
+    //int Q = 9;
 
     double Fx = 1.0;
     double Fy = 1.0;
     double Fz = 1.0;
     
-    double *dist;    
-    double * Velocity;
-    
-    dist = new double [19*Np];
-    Velocity = new double [3*Np];
-
+    auto dist = new double [19*Np];
+    //auto Velocity  = new double [3*Np
 
     for (int n=0; n<Np; n++){
     	dist[n] = 0.3333333333333333;

tests/TestInterfaceSpeed.cpp

@@ -24,7 +24,7 @@ int main (int argc, char *argv[])
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
 
-	int i,j,k,n;
+	int i,j,k;
 
     // Load inputs
 	string FILENAME = argv[1];
@@ -36,7 +36,7 @@ int main (int argc, char *argv[])
     int Ny = domain_db->getVector<int>( "n" )[1];
     int Nz = domain_db->getVector<int>( "n" )[2];
 
-    std::shared_ptr<Domain> Dm(new Domain(domain_db,comm));
+    auto Dm = std::make_shared<Domain>(domain_db,comm);
 
     Nx+=2; Ny+=2; Nz+=2;
 
@@ -44,7 +44,7 @@ int main (int argc, char *argv[])
 
 	Dm->CommInit();
 
-	std::shared_ptr<TwoPhase> Averages(new TwoPhase(Dm));
+	auto Averages = std::make_shared<TwoPhase>(Dm);
 	int timestep=0;
 
 	double Cx,Cy,Cz;

tests/TestMap.cpp

@@ -56,11 +56,7 @@ int main(int argc, char **argv)
 	    int Nx = db->getVector<int>( "n" )[0];
 	    int Ny = db->getVector<int>( "n" )[1];
 	    int Nz = db->getVector<int>( "n" )[2];
-	    int nprocx = db->getVector<int>( "nproc" )[0];
-	    int nprocy = db->getVector<int>( "nproc" )[1];
-	    int nprocz = db->getVector<int>( "nproc" )[2];
-
-		std::shared_ptr<Domain> Dm(new Domain(db,comm));
+		auto Dm = std::make_shared<Domain>(db,comm);
 
 		Nx += 2;
 		Ny += 2;

tests/TestMassConservationD3Q7.cpp

@@ -41,7 +41,7 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 				int jglobal= j+(Ny-2)*ColorModel.Mask->jproc();
 				int kglobal= k+(Nz-2)*ColorModel.Mask->kproc();
 				// Initialize phase position field for parallel bubble test
-				if (jglobal < 40){
+				if (kglobal < 40){
 					ColorModel.Mask->id[n] = 0;
 				}
 				else if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
@@ -66,9 +66,6 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 
 int main(int argc, char **argv)
 {
-	//*****************************************
-	// ***** MPI STUFF ****************
-	//*****************************************
 	// Initialize MPI
 	int rank,nprocs;
 	MPI_Init(&argc,&argv);
@@ -76,19 +73,7 @@ int main(int argc, char **argv)
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
 	// parallel domain size (# of sub-domains)
-	int nprocx,nprocy,nprocz;
-	int iproc,jproc,kproc;
-	int sendtag,recvtag;
-	//*****************************************
-	// MPI ranks for all 18 neighbors
-	//**********************************
-	int rank_x,rank_y,rank_z,rank_X,rank_Y,rank_Z;
-	int rank_xy,rank_XY,rank_xY,rank_Xy;
-	int rank_xz,rank_XZ,rank_xZ,rank_Xz;
-	int rank_yz,rank_YZ,rank_yZ,rank_Yz;
-	//**********************************
-	MPI_Request req1[18],req2[18];
-	MPI_Status stat1[18],stat2[18];
+	int CleanCheck = 0;
 
 	if (rank == 0){
 		printf("********************************************************\n");
@@ -110,7 +95,6 @@ int main(int argc, char **argv)
 		Ny = CM.Ny;
 		Nz = CM.Nz;
 		N = Nx*Ny*Nz;
-	int dist_mem_size = N*sizeof(double);
 
 		//CM.ReadInput(); 
 		double radius=0.4*double(Nx);
@@ -139,15 +123,13 @@ int main(int argc, char **argv)
 		 */
 
 		CM.timestepMax = 2;
+		CM.timestep = 0;
 		CM.Run();
 		int D3Q7[7][3]={{0,0,0},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1}};
 		// Compare and make sure mass is conserved at every lattice site
-	double *Error;
-	Error = new double [N];
-	double *A_q, *B_q;
-	A_q = new double [7*Np];
-	B_q = new double [7*Np];
-	bool CleanCheck = true;
+		auto Error = new double[N];
+		auto A_q = new double[7*Np];
+		//auto B_q = new double[7*Np];
 		double original,final, sum_q;
 		double total_mass_A_0 = 0.0;
 		double total_mass_B_0= 0.0;
@@ -188,14 +170,14 @@ int main(int argc, char **argv)
 						}
 						Error[n] = sum_q - original;
 
-					       /*if (fabs(DenFinal[idx] - DenOriginal[idx]) > 1e-15){		      
+						if (fabs(DenFinal[idx] - DenOriginal[idx]) > 1e-15){		      
 							//if (CM.Dm->id[n] == 0) printf("Solid phase! \n");
 							//if (CM.Dm->id[n] == 1) printf("Wetting phase! \n");
 							//if (CM.Dm->id[n] == 2) printf("Non-wetting phase! \n");							
-						printf("Mass not conserved: WP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
+							//printf("Mass not conserved: WP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
 							CleanCheck=false;
 							Error[n] += final-original;
-						}*/
+						}
 						final = DenFinal[Np+idx];
 						if (final < 0.0) count_negative_B++;
 						original = DenOriginal[Np+idx];
@@ -205,10 +187,11 @@ int main(int argc, char **argv)
 							//if (CM.Dm->id[n] == 0) printf("Solid phase! \n");
 							//if (CM.Dm->id[n] == 1) printf("Wetting phase! \n");
 							//if (CM.Dm->id[n] == 2) printf("Non-wetting phase! \n");
-						printf("Mass not conserved: NWP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
+							//printf("Mass not conserved: NWP density, site=%i,%i,%i, original = %f, final = %f \n",i,j,k,original,final);
 							CleanCheck=false;
 							Error[n] += final-original;
-						}*/
+						} 
+						*/
 					}
 				}
 			}
@@ -219,13 +202,13 @@ int main(int argc, char **argv)
 		printf("Global mass difference B = %.5g\n",total_mass_B_1-total_mass_B_0);
 
 		if (count_negative_A > 0 ||count_negative_B > 0) CleanCheck=1;
-	if (fabs(total_mass_A_1-total_mass_A_0) > 1.0e-15||fabs(total_mass_B_1-total_mass_B_0) > 1.0e-15 ) CleanCheck=2;
+		if (fabs(total_mass_A_1-total_mass_A_0) > 1.0e-8 || fabs(total_mass_B_1-total_mass_B_0) > 1.0e-8) CleanCheck=2;
 
-	/*	
 		FILE *OUTFILE;
 		OUTFILE = fopen("error.raw","wb");
 		fwrite(Error,8,N,OUTFILE);
 		fclose(OUTFILE);	
+		/*	
 
 	if (rank==0) printf("Checking that the correct velocity is retained \n");
 	// Swap convention is observed -- velocity is negative
@@ -276,7 +259,7 @@ int main(int argc, char **argv)
 		}
 	}
 		 */
-	if (CleanCheck){
+		if (CleanCheck == 0){
 			if (rank==0) printf("Test passed: mass conservation for D3Q7 \n");
 		}
 		else {

tests/TestNetcdf.cpp

@@ -14,7 +14,6 @@ void load( const std::string& );
 void test_NETCDF( UnitTest& ut )
 {
     const int rank = comm_rank( MPI_COMM_WORLD );
-    const int size = comm_size( MPI_COMM_WORLD );
     int nprocx = 2;
     int nprocy = 2;
     int nprocz = 2;

tests/TestSubphase.cpp

@@ -60,13 +60,11 @@ int main(int argc, char **argv)
 		}
 
 		// Get the rank info
-		std::shared_ptr<Domain> Dm(new Domain(db,comm));
-		//   const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
-		std::shared_ptr<SubPhase> Averages(new SubPhase(Dm));
+		auto Dm = std::make_shared<Domain>(db,comm);
+		auto Averages = std::make_shared<SubPhase>(Dm);
 		Nx += 2;
 		Ny += 2;
 		Nz += 2;
-		int N = Nx*Ny*Nz;
 		//.......................................................................
 		for ( k=1;k<Nz-1;k++){
 			for ( j=1;j<Ny-1;j++){

tests/TestTopo3D.cpp

@@ -60,12 +60,11 @@ int main(int argc, char **argv)
 		}
 
 		// Get the rank info
-		std::shared_ptr<Domain> Dm(new Domain(db,comm));
+		auto Dm = std::make_shared<Domain>(db,comm);
 
 		Nx += 2;
 		Ny += 2;
 		Nz += 2;
-		int N = Nx*Ny*Nz;
 		//.......................................................................
 		for ( k=1;k<Nz-1;k++){
 			for ( j=1;j<Ny-1;j++){

tests/TestTorus.cpp

@@ -60,13 +60,12 @@ int main(int argc, char **argv)
     }
 
     // Get the rank info
-		std::shared_ptr<Domain> Dm(new Domain(db,comm));
+    auto Dm = std::make_shared<Domain>(db,comm);
  //   const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
-		std::shared_ptr<TwoPhase> Averages(new TwoPhase(Dm));
+    auto Averages = std::make_shared<TwoPhase>(Dm);
 	Nx += 2;
 	Ny += 2;
 	Nz += 2;
-	int N = Nx*Ny*Nz;
 	//.......................................................................
 	for ( k=1;k<Nz-1;k++){
 		for ( j=1;j<Ny-1;j++){
@@ -142,7 +141,7 @@ int main(int argc, char **argv)
 		}
 	}
 
-    double beta = 0.95;
+    //double beta = 0.95;
 	if (rank==0) printf("initializing the system \n");
 
 	Averages->UpdateSolid();

tests/TestTorusEvolve.cpp

@@ -60,12 +60,11 @@ int main(int argc, char **argv)
     }
 
     // Get the rank info
-    std::shared_ptr<Domain> Dm(new Domain(db,comm));
+    auto Dm = std::make_shared<Domain>(db,comm);
 
     Nx += 2;
     Ny += 2;
     Nz += 2;
-    int N = Nx*Ny*Nz;
     //.......................................................................
     for ( k=1;k<Nz-1;k++){
     	for ( j=1;j<Ny-1;j++){
@@ -98,14 +97,13 @@ int main(int argc, char **argv)
 	//.......................................................................
 	// Assign the phase ID field based and the signed distance
 	//.......................................................................
-    double R1,R2,R;
     double CX,CY,CZ; //CY1,CY2;
     CX=Nx*nprocx*0.5;
     CY=Ny*nprocy*0.5;
     CZ=Nz*nprocz*0.5;
-    R1 = (Nx-2)*nprocx*0.3; // middle radius
-    R2 = (Nx-2)*nprocx*0.1; // donut thickness
-    R = 0.4*nprocx*(Nx-2);
+    auto R1 = (Nx-2)*nprocx*0.3; // middle radius
+    auto R2 = (Nx-2)*nprocx*0.1; // donut thickness
+    //auto R = 0.4*nprocx*(Nx-2);
     
     Minkowski Object(Dm);
 

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,16 +23,15 @@ using namespace std;
 
 int main(int argc, char **argv)
 {
-  // Initialize MPI
-  int provided_thread_support = -1;
-  MPI_Init_thread(&argc,&argv,MPI_THREAD_MULTIPLE,&provided_thread_support);
+  // Initialize MPI and error handlers
+  Utilities::startup( argc, argv );
+
+  { // Limit scope so variables that contain communicators will free before MPI_Finialize
+
     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){
 	    printf("********************************************************\n");
@@ -40,7 +39,7 @@ int main(int argc, char **argv)
 	    printf("********************************************************\n");
     }
     // Initialize compute device
-	int device=ScaLBL_SetDevice(rank);
+    ScaLBL_SetDevice(rank);
     ScaLBL_DeviceBarrier();
     MPI_Barrier(comm);
 
@@ -64,10 +63,13 @@ int main(int argc, char **argv)
     PROFILE_STOP("Main");
     PROFILE_SAVE("lbpm_color_simulator",1);
     // ****************************************************
+
     MPI_Barrier(comm);
-  } // Limit scope so variables that contain communicators will free before MPI_Finialize
     MPI_Comm_free(&comm);
-  MPI_Finalize();
+
+  } // Limit scope so variables that contain communicators will free before MPI_Finialize
+
+  Utilities::shutdown();
 }
 
 

tests/lbpm_greyscale_simulator.cpp

@@ -0,0 +1,63 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/stat.h>
+#include <iostream>
+#include <exception>
+#include <stdexcept>
+#include <fstream>
+
+#include "common/ScaLBL.h"
+#include "common/Communication.h"
+#include "common/MPI_Helpers.h"
+#include "models/GreyscaleModel.h"
+//#define WRITE_SURFACES
+
+/*
+ * Simulator for two-phase flow in porous media
+ * James E. McClure 2013-2014
+ */
+
+using namespace std;
+
+
+int main(int argc, char **argv)
+{
+	//*****************************************
+	// ***** MPI STUFF ****************
+	//*****************************************
+	// 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);
+	{
+		// parallel domain size (# of sub-domains)
+
+		if (rank == 0){
+			printf("********************************************************\n");
+			printf("Running Greyscale Single Phase Permeability Calculation \n");
+			printf("********************************************************\n");
+		}
+		// Initialize compute device
+		int device=ScaLBL_SetDevice(rank);
+		NULL_USE(device);
+		ScaLBL_DeviceBarrier();
+		MPI_Barrier(comm);
+		
+		ScaLBL_GreyscaleModel Greyscale(rank,nprocs,comm);
+		auto filename = argv[1];
+		Greyscale.ReadParams(filename);
+		Greyscale.SetDomain();    // this reads in the domain 
+		Greyscale.ReadInput();
+		Greyscale.Create();       // creating the model will create data structure to match the pore structure and allocate variables
+		Greyscale.Initialize();   // initializing the model will set initial conditions for variables
+		Greyscale.Run();	 
+		//Greyscale.VelocityField();
+		//Greyscale.WriteDebug();
+	}
+	// ****************************************************
+	MPI_Barrier(comm);
+	MPI_Finalize();
+	// ****************************************************
+}

tests/lbpm_morph_pp.cpp

@@ -34,7 +34,6 @@ int main(int argc, char **argv)
 		//.......................................................................
 	        int n, nx, ny, nz;
 		char LocalRankFilename[40];
-		char FILENAME[128];
 
 		string filename;
 		double SW;
@@ -239,10 +238,9 @@ int main(int argc, char **argv)
 		}
 		MPI_Barrier(comm);
 
-		sprintf(FILENAME,READFILE.c_str());
-		sprintf(FILENAME+strlen(FILENAME),".morph.raw");
-		if (rank==0) printf("Writing file to: %s \n", FILENAME);
-		Mask->AggregateLabels(FILENAME);
+        auto filename2 = READFILE + ".morph.raw";
+		if (rank==0) printf("Writing file to: %s \n", filename2.c_str());
+		Mask->AggregateLabels(filename2);
 	}
 
 	MPI_Barrier(comm);

tests/lbpm_morphdrain_pp.cpp

@@ -32,10 +32,7 @@ int main(int argc, char **argv)
 		//.......................................................................
 		// Reading the domain information file
 		//.......................................................................
-		int n, nprocx, nprocy, nprocz, nx, ny, nz;
-		char LocalRankString[8];
 		char LocalRankFilename[40];
-		char FILENAME[128];
 
 		string filename;
 		double SW,Rcrit_new;
@@ -43,8 +40,10 @@ int main(int argc, char **argv)
 			filename=argv[1];
 			Rcrit_new=0.f; 
 			//SW=strtod(argv[2],NULL);
+		} else {
+            ERROR("No input database provided\n");
         }
-		else ERROR("No input database provided\n");
+        NULL_USE( Rcrit_new );
 		// read the input database 
 		auto db = std::make_shared<Database>( filename );
 		auto domain_db = db->getDatabase( "Domain" );
@@ -62,19 +61,16 @@ int main(int argc, char **argv)
 		if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
 		//	GenerateResidual(id,nx,ny,nz,Saturation);
 		
-		nx = size[0];
-		ny = size[1];
-		nz = size[2];
-		nprocx = nproc[0];
-		nprocy = nproc[1];
-		nprocz = nproc[2];
+		int nx = size[0];
+		int ny = size[1];
+		int nz = size[2];
 
-		int N = (nx+2)*(ny+2)*(nz+2);
+		size_t N = (nx+2)*(ny+2)*(nz+2);
 
 		std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
 		std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
 		//		std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
-		for (n=0; n<N; n++) Dm->id[n]=1;
+		for (size_t n=0; n<N; n++) Dm->id[n]=1;
 		Dm->CommInit();
 
 		signed char *id;
@@ -116,7 +112,6 @@ int main(int argc, char **argv)
 		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;
 					// Initialize distance to +/- 1
 					SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
 				}
@@ -158,7 +153,6 @@ int main(int argc, char **argv)
 			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;
 						// Initialize distance to +/- 1
 						SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
 					}
@@ -204,10 +198,9 @@ int main(int argc, char **argv)
 		}
 		MPI_Barrier(comm);
 
-		sprintf(FILENAME,READFILE.c_str());
-		sprintf(FILENAME+strlen(FILENAME),".morphdrain.raw");
-		if (rank==0) printf("Writing file to: %s \n", FILENAME);
-		Mask->AggregateLabels(FILENAME);
+        auto filename2 = READFILE + ".morphdrain.raw";
+		if (rank==0) printf("Writing file to: %s \n", filename2.data() );
+		Mask->AggregateLabels( filename2 );
 	}
 
 	MPI_Barrier(comm);

tests/lbpm_morphopen_pp.cpp

@@ -32,10 +32,7 @@ int main(int argc, char **argv)
 		//.......................................................................
 		// Reading the domain information file
 		//.......................................................................
-		int n, nprocx, nprocy, nprocz, nx, ny, nz;
-		char LocalRankString[8];
 		char LocalRankFilename[40];
-		char FILENAME[128];
 
 		string filename;
 		double SW,Rcrit_new;
@@ -45,6 +42,7 @@ int main(int argc, char **argv)
 			//SW=strtod(argv[2],NULL);
 		}
 		else ERROR("No input database provided\n");
+        NULL_USE( Rcrit_new );
 		// read the input database 
 		auto db = std::make_shared<Database>( filename );
 		auto domain_db = db->getDatabase( "Domain" );
@@ -69,19 +67,16 @@ int main(int argc, char **argv)
 		if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
 		//	GenerateResidual(id,nx,ny,nz,Saturation);
 		
-		nx = size[0];
-		ny = size[1];
-		nz = size[2];
-		nprocx = nproc[0];
-		nprocy = nproc[1];
-		nprocz = nproc[2];
+		int nx = size[0];
+		int ny = size[1];
+		int nz = size[2];
 
-		int N = (nx+2)*(ny+2)*(nz+2);
+		size_t N = (nx+2)*(ny+2)*(nz+2);
 
 		std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
 		std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
 		//		std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
-		for (n=0; n<N; n++) Dm->id[n]=1;
+		for (size_t n=0; n<N; n++) Dm->id[n]=1;
 		Dm->CommInit();
 
 		signed char *id;
@@ -119,7 +114,6 @@ int main(int argc, char **argv)
 		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;
 					// Initialize distance to +/- 1
 					SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
 				}
@@ -161,7 +155,6 @@ int main(int argc, char **argv)
 			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;
 						// Initialize distance to +/- 1
 						SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
 					}
@@ -207,10 +200,9 @@ int main(int argc, char **argv)
 		}
 		MPI_Barrier(comm);
 
-		sprintf(FILENAME,READFILE.c_str());
-		sprintf(FILENAME+strlen(FILENAME),".morphopen.raw");
-		if (rank==0) printf("Writing file to: %s \n", FILENAME);
-		Mask->AggregateLabels(FILENAME);
+        auto filename2 = READFILE + ".morphopen.raw";
+		if (rank==0) printf("Writing file to: %s \n", filename2.data());
+		Mask->AggregateLabels(filename2);
 	}
 
 	MPI_Barrier(comm);

tests/lbpm_permeability_simulator.cpp

@@ -23,9 +23,6 @@ using namespace std;
 
 int main(int argc, char **argv)
 {
-	//*****************************************
-	// ***** MPI STUFF ****************
-	//*****************************************
 	// Initialize MPI
 	int rank,nprocs;
 	MPI_Init(&argc,&argv);
@@ -33,10 +30,6 @@ int main(int argc, char **argv)
 	MPI_Comm_rank(comm,&rank);
 	MPI_Comm_size(comm,&nprocs);
 	{
-		// parallel domain size (# of sub-domains)
-		int nprocx,nprocy,nprocz;
-		int iproc,jproc,kproc;
-
 		if (rank == 0){
 			printf("********************************************************\n");
 			printf("Running Single Phase Permeability Calculation \n");
@@ -44,10 +37,10 @@ int main(int argc, char **argv)
 		}
 		// Initialize compute device
 		int device=ScaLBL_SetDevice(rank);
+        NULL_USE( device );
 		ScaLBL_DeviceBarrier();
 		MPI_Barrier(comm);
 		
-		
 		ScaLBL_MRTModel MRT(rank,nprocs,comm);
 		auto filename = argv[1];
 		MRT.ReadParams(filename);

tests/lbpm_refine_pp.cpp

@@ -12,6 +12,7 @@
 #include "common/Communication.h"
 #include "common/Domain.h"
 #include "analysis/pmmc.h"
+#include "analysis/distance.h"
 
 int main(int argc, char **argv)
 {
@@ -26,10 +27,9 @@ int main(int argc, char **argv)
 		//.......................................................................
 		// Reading the domain information file
 		//.......................................................................
-		int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
 		double Lx, Ly, Lz;
+        Lx = Ly = Lz = 1.0;
 		int i,j,k,n;
-		int BC=0;
 
 		string filename;
 		if (argc > 1){
@@ -41,18 +41,18 @@ int main(int argc, char **argv)
 		auto domain_db = db->getDatabase( "Domain" );
 
 		// Read domain parameters
-		auto L = domain_db->getVector<double>( "L" );
 		auto size = domain_db->getVector<int>( "n" );
 		auto nproc = domain_db->getVector<int>( "nproc" );
 		auto ReadValues = domain_db->getVector<char>( "ReadValues" );
 		auto WriteValues = domain_db->getVector<char>( "WriteValues" );
 		
-		nx = size[0];
-		ny = size[1];
-		nz = size[2];
-		nprocx = nproc[0];
-		nprocy = nproc[1];
-		nprocz = nproc[2];
+		int nx = size[0];
+		int ny = size[1];
+		int nz = size[2];
+		int nprocx = nproc[0];
+		int nprocy = nproc[1];
+		int nprocz = nproc[2];
+		int BoundaryCondition=0;
 
 		// Check that the number of processors >= the number of ranks
 		if ( rank==0 ) {
@@ -64,16 +64,26 @@ int main(int argc, char **argv)
 			ERROR("Insufficient number of processors");
 		}
 
-		char LocalRankFilename[40];
+		//Domain Mask(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
+		Domain Mask(domain_db,MPI_COMM_WORLD);
+		if (domain_db->keyExists( "Filename" )){
+			auto Filename = domain_db->getScalar<std::string>( "Filename" );
+		        if (rank==0) printf("Reading domain from %s \n",Filename.c_str());
+			Mask.Decomp(Filename);
+			if (rank==0) printf("Complete. \n");
+		}
+		else{
+			Mask.ReadIDs();
+		}
+		Mask.CommInit();
 
-		int rnx,rny,rnz;
-		rnx=2*nx;
-		rny=2*ny;
-		rnz=2*nz;
+		char LocalRankFilename[40];
+		int rnx=2*nx;
+		int rny=2*ny;
+		int rnz=2*nz;
 
 		if (rank==0) printf("Refining mesh to %i x %i x %i \n",rnx,rny,rnz);
 
-		int BoundaryCondition=0;
 		Domain Dm(rnx,rny,rnz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
 
 		// Communication the halos
@@ -85,6 +95,7 @@ int main(int argc, char **argv)
 		int N = nx*ny*nz;
 
 		// Define communication sub-domain -- everywhere
+		if (rank==0) printf("Initialize refined domain \n");
 		for (int k=0; k<rnz; k++){
 			for (int j=0; j<rny; j++){
 				for (int i=0; i<rnx; i++){
@@ -95,8 +106,37 @@ int main(int argc, char **argv)
 		}
 		Dm.CommInit();
 		
+		// Generate the signed distance map
+		// Initialize the domain and communication
+		Array<char> Labels(nx,ny,nz);
 		DoubleArray SignDist(nx,ny,nz);
-		// Read the signed distance from file
+
+		// Solve for the position of the solid phase
+		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;
+					// Initialize the solid phase
+					signed char label = Mask.id[n];
+					if (label > 0)		Labels(i,j,k) = 1;
+					else	     		Labels(i,j,k) = 0;
+				}
+			}
+		}
+		// Initialize the signed distance function
+		for (int k=0;k<nz;k++){
+			for (int j=0;j<ny;j++){
+				for (int i=0;i<nx;i++){
+					// Initialize distance to +/- 1
+					SignDist(i,j,k) = 2.0*double(Labels(i,j,k))-1.0;
+				}
+			}
+		}
+	//	MeanFilter(Averages->SDs);
+		if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
+		CalcDist(SignDist,Labels,Mask);
+		
+		/*		// Read the signed distance from file
 		sprintf(LocalRankFilename,"SignDist.%05i",rank);
 		FILE *DIST = fopen(LocalRankFilename,"rb");
 		size_t ReadSignDist;
@@ -112,7 +152,7 @@ int main(int argc, char **argv)
 		ReadLabels=fread(Labels,1,N,LABELS);
 		if (ReadLabels != size_t(N)) printf("lbpm_refine_pp: Error reading ID  (rank=%i)\n",rank);
 		fclose(LABELS);
-
+*/
 		if ( rank==0 )   printf("Set up Domain, read input distance \n");
 
 		DoubleArray RefinedSignDist(rnx,rny,rnz);
@@ -128,13 +168,12 @@ int main(int argc, char **argv)
 			}
 		}
 
-		int ri,rj,rk,rn; //refined mesh indices
 		//char *RefineLabel;
 		//RefineLabel = new char [rnx*rny*rnz];
 		Array <char> RefineLabel(rnx,rny,rnz);
-		for (rk=1; rk<rnz-1; rk++){
-			for (rj=1; rj<rny-1; rj++){
-				for (ri=1; ri<rnx-1; ri++){
+		for (int rk=1; rk<rnz-1; rk++){
+			for (int rj=1; rj<rny-1; rj++){
+				for (int ri=1; ri<rnx-1; ri++){
 					n = rk*rnx*rny+rj*rnx+ri;
 					// starting node for each processor matches exactly
 					i = (ri-1)/2+1;
@@ -148,23 +187,24 @@ int main(int argc, char **argv)
 					pt.y=0.5*(rj-1)+1.f;
 					pt.z=0.5*(rk-1)+1.f;
 					RefinedSignDist(ri,rj,rk) = LocalApprox.eval(pt);
-					RefineLabel(ri,rj,rk) = Labels[k*nx*ny+j*nx+i]; 
+					RefineLabel(ri,rj,rk) = Labels(i,j,k);
+					Dm.id[n] = Labels(i,j,k);
 				}
 			}
 		}
 		fillData.fill(RefinedSignDist);
-		//	sprintf(LocalRankFilename,"ID.%05i",rank);
-		//FILE *ID = fopen(LocalRankFilename,"wb");
-		//fwrite(id,1,N,ID);
-		//fclose(ID);
-/*
-		sprintf(LocalRankFilename,"RefineDist.%05i",rank);
-		FILE *REFINEDIST = fopen(LocalRankFilename,"wb");
-		fwrite(RefinedSignDist.data(),8,rnx*rny*rnz,REFINEDIST);
-		fclose(REFINEDIST);
-*/
-		if ( rank==0 )   printf("Write output \n");
 
+
+		if (domain_db->keyExists( "Filename" )){
+			auto Filename = domain_db->getScalar<std::string>( "Filename" );
+			if ( rank==0 )   printf("Write output \n");
+			Dm.AggregateLabels("id_2x.raw");
+			Mask.AggregateLabels("id.raw");
+			//FILE *WRITEID = fopen("refine.raw","wb");
+			//fwrite(RefineLabel.data(),1,rnx*rny*rnz,WRITEID);
+			//fclose(WRITEID);
+		}
+		else{
 			DoubleArray BlockDist(nx,ny,nz);
 			FILE *WRITEID, *REFINEDIST;
 			char * id;
@@ -188,17 +228,7 @@ int main(int argc, char **argv)
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
 
-/*		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]= 0;
-				}
-			}
-		}
-		*/
+
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
@@ -219,17 +249,7 @@ int main(int argc, char **argv)
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
 
-/*		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]=0;
-				}
-			}
-		}
-		*/
+
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
@@ -251,17 +271,7 @@ int main(int argc, char **argv)
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
 
-/*		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]=0;
-				}
-			}
-		}
-		*/
+
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
@@ -281,18 +291,7 @@ int main(int argc, char **argv)
 			REFINEDIST = fopen(LocalRankFilename,"wb");
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
-/*
-		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]=0;
-				}
-			}
-		}
-		*/
+
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
@@ -312,18 +311,7 @@ int main(int argc, char **argv)
 			REFINEDIST = fopen(LocalRankFilename,"wb");
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
-/*
-		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]=0;
-				}
-			}
-		}
-		*/
+
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
@@ -344,17 +332,6 @@ int main(int argc, char **argv)
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
 
-/*		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]=0;
-				}
-			}
-		}
-		*/
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
@@ -374,18 +351,7 @@ int main(int argc, char **argv)
 			REFINEDIST = fopen(LocalRankFilename,"wb");
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
-/*
-		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]=0;
-				}
-			}
-		}
-		*/
+
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
@@ -407,23 +373,11 @@ int main(int argc, char **argv)
 			fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 			fclose(REFINEDIST);
 
-/*		for (int k=0; k<nz; k++){
-			for (int j=0; j<ny; j++){
-				for (int i=0; i<nx; i++){
-					if (BlockDist(i,j,k) > 0.f)
-						id[k*nx*ny + j*nx + i]=2;
-					else
-						id[k*nx*ny + j*nx + i]=0;
-				}
-			}
-		}
-		*/
 			sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 			WRITEID = fopen(LocalRankFilename,"wb");
 			fwrite(id,1,nx*ny*nz,WRITEID);
 			fclose(WRITEID);
-
-
+		}
 	}
 	MPI_Barrier(comm);
 	MPI_Finalize();

tests/lbpm_serial_decomp.cpp

@@ -47,11 +47,7 @@ int main(int argc, char **argv)
 	//.......................................................................
 	// Reading the domain information file
 	//.......................................................................
-	int nprocs, nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
-	double Lx, Ly, Lz;
-	int64_t Nx,Ny,Nz;
 	int64_t i,j,k,n;
-	int BC=0;
 	int64_t xStart,yStart,zStart;
 	int checkerSize;
 	int inlet_count_x, inlet_count_y, inlet_count_z;
@@ -112,25 +108,25 @@ int main(int argc, char **argv)
 		ReadType = "8bit";
 	}
 
-	nx = size[0];
-	ny = size[1];
-	nz = size[2];
-	nprocx = nproc[0];
-	nprocy = nproc[1];
-	nprocz = nproc[2];
-	Nx = SIZE[0];
-	Ny = SIZE[1];
-	Nz = SIZE[2];
+	int nx = size[0];
+	int ny = size[1];
+	int nz = size[2];
+	int nprocx = nproc[0];
+	int nprocy = nproc[1];
+	int nprocz = nproc[2];
+	long int Nx = SIZE[0];
+	long int Ny = SIZE[1];
+	long int Nz = SIZE[2];
 
 	printf("Input media: %s\n",Filename.c_str());
 	printf("Relabeling %lu values\n",ReadValues.size());
-	for (int idx=0; idx<ReadValues.size(); idx++){
+	for (size_t idx=0; idx<ReadValues.size(); idx++){
 		int oldvalue=ReadValues[idx];
 		int newvalue=WriteValues[idx];
 		printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
 	}
 
-	nprocs=nprocx*nprocy*nprocz;
+	int nprocs=nprocx*nprocy*nprocz;
 
 	char *SegData = NULL;
 	// Rank=0 reads the entire segmented data and distributes to worker processes
@@ -172,7 +168,7 @@ int main(int argc, char **argv)
 				n = k*Nx*Ny+j*Nx+i;
 				//char locval = loc_id[n];
 				char locval = SegData[n];
-				for (int idx=0; idx<ReadValues.size(); idx++){
+				for (size_t idx=0; idx<ReadValues.size(); idx++){
 					signed char oldvalue=ReadValues[idx];
 					signed char newvalue=WriteValues[idx];
 					if (locval == oldvalue){
@@ -185,10 +181,10 @@ int main(int argc, char **argv)
 		}
 	}
 	if (rank==0){
-		for (int idx=0; idx<ReadValues.size(); idx++){
+		for (size_t idx=0; idx<ReadValues.size(); idx++){
 			long int label=ReadValues[idx];
 			long int count=LabelCount[idx];
-			printf("Label=%d, Count=%d \n",label,count);
+			printf("Label=%ld, Count=%ld \n",label,count);
 		}
 	}
 	
@@ -215,7 +211,7 @@ int main(int argc, char **argv)
 		printf("Checkerboard pattern at y inlet for %i layers \n",inlet_count_y);
 		// use checkerboard pattern
 		for (int k = 0; k<Nz; k++){
-			for (int j = yStart; i < yStart+inlet_count_y; j++){
+			for (int j = yStart; j < yStart+inlet_count_y; j++){
 				for (int i = 0; i<Nx; i++){
 					if ( (i/checkerSize + k/checkerSize)%2 == 0){
 						// void checkers
@@ -272,7 +268,7 @@ int main(int argc, char **argv)
 		printf("Checkerboard pattern at y outlet for %i layers \n",outlet_count_y);
 		// use checkerboard pattern
 		for (int k = 0; k<Nz; k++){
-			for (int j = yStart + ny*nprocy - outlet_count_y; i < yStart + ny*nprocy; j++){
+			for (int j = yStart + ny*nprocy - outlet_count_y; j < yStart + ny*nprocy; j++){
 				for (int i = 0; i<Nx; i++){
 					if ( (i/checkerSize + k/checkerSize)%2 == 0){
 						// void checkers
@@ -306,9 +302,6 @@ int main(int argc, char **argv)
 		}
 	}
 
-	// Get the rank info
-	int64_t N = (nx+2)*(ny+2)*(nz+2);
-
 	// number of sites to use for periodic boundary condition transition zone
 	int64_t z_transition_size = (nprocz*nz - (Nz - zStart))/2;
 	if (z_transition_size < 0) z_transition_size=0;

tests/lbpm_uCT_pp.cpp

@@ -61,7 +61,7 @@ int main(int argc, char **argv)
         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 BoundaryCondition = domain_db->getScalar<int>( "BC" );
         int nx = size[0];
         int ny = size[1];
         int nz = size[2];
@@ -200,10 +200,10 @@ int main(int argc, char **argv)
         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++) {
-		  float x= float(Dm[0]->iproc()*Nx[0]+i-1);
+                  //float x= float(Dm[0]->iproc()*Nx[0]+i-1);
                   float y= float (Dm[0]->jproc()*Ny[0]+j-1);
                   float z= float(Dm[0]->kproc()*Nz[0]+k-1);
-		  float cx = float(center[0] - offset[0]);
+                  //float cx = float(center[0] - offset[0]);
                   float cy = float(center[1] - offset[1]);
                   float cz = float(center[2] - offset[2]);
                     // distance from the center line 

tests/pmmc_cylinder.cpp

@@ -11,17 +11,11 @@
 
 int main (int argc, char **argv)
 {
-	
 	//  printf("Radius = %s \n,"RADIUS);  
-	int SIZE = N*N*N;
 	int Nx,Ny,Nz;
 	Nx = Ny = Nz = N;
 	int i,j,k,p,q,r;
 	
-//	double *Solid; // cylinder
-//	double *Phase; // region of the cylinder	
-//	Solid = new double [SIZE];
-//	Phase = new double [SIZE];
 	DoubleArray SignDist(Nx,Ny,Nz);
 	DoubleArray Phase(Nx,Ny,Nz);
 	double fluid_isovalue = 0.0;
@@ -36,9 +30,6 @@ int main (int argc, char **argv)
 	//...........................................................................
 	double awn,ans,aws,lwns,nwp_volume;
 	double As;
-	double dEs,dAwn,dAns;			 // Global surface energy (calculated by rank=0)
-	double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;	
-	double As_global;
 	//	bool add=1;			// Set to false if any corners contain nw-phase ( F > fluid_isovalue)
 	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};  // cube corners
 	//	int count_in=0,count_out=0;
@@ -75,7 +66,6 @@ int main (int argc, char **argv)
 	int n_local_nws_pts;
 	
 	int c;
-	int newton_steps = 0;
 	//...........................................................................
 	int ncubes = (Nx-2)*(Ny-2)*(Nz-2);	// Exclude the "upper" halo
 	IntArray cubeList(3,ncubes);