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

analysis/SubPhase.cpp

@@ -28,8 +28,15 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
 
 	//.........................................
 	if (Dm->rank()==0){
+		bool WriteHeader=false;
+		SUBPHASE = fopen("subphase.csv","r");
+		if (SUBPHASE != NULL)
+			fclose(SUBPHASE);
+		else
+			WriteHeader=true;
+
 		SUBPHASE = fopen("subphase.csv","a+");
-		if (ftell(SUBPHASE) == 0)
+		if (WriteHeader)
 		{
 			// If timelog is empty, write a short header to list the averages
 			//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
@@ -141,10 +148,11 @@ void SubPhase::Basic(){
 	if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
 
 	imin=jmin=1;
-	// If inlet 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) kmin = Dm->inlet_layers_z;
+	if (Dm->outlet_layers_z > 0) kmax = Dm->outlet_layers_z;
 	
 	nb.reset(); wb.reset();
 /*

analysis/runAnalysis.cpp

@@ -389,9 +389,9 @@ public:
             // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
         }
         if ( matches(type,AnalysisType::ComputeAverages) ) {
-            PROFILE_START("Compute subphase",1);
+            PROFILE_START("Compute basic averages",1);
             Averages.Basic();
-            PROFILE_STOP("Compute subphase",1);
+            PROFILE_STOP("Compute basic averages",1);
         }
     }
 private:
@@ -402,6 +402,29 @@ private:
     double beta;
 };
 
+class SubphaseWorkItem: public ThreadPool::WorkItemRet<void>
+{
+public:
+	SubphaseWorkItem( AnalysisType type_, int timestep_, SubPhase& Averages_ ):
+                type(type_), timestep(timestep_), Averages(Averages_){ }
+    ~SubphaseWorkItem() { }
+    virtual void run() {
+
+        if ( matches(type,AnalysisType::ComputeAverages) ) {
+            PROFILE_START("Compute subphase",1);
+            Averages.Full();
+            PROFILE_STOP("Compute subphase",1);
+        }
+    }
+private:
+    SubphaseWorkItem();
+    AnalysisType type;
+    int timestep;
+    SubPhase& Averages;
+    double beta;
+};
+
+
 
 /******************************************************************
  *  MPI comm wrapper for use with analysis                         *
@@ -468,6 +491,7 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
     ThreadPool::thread_id_t d_wait_analysis;
     ThreadPool::thread_id_t d_wait_vis;
     ThreadPool::thread_id_t d_wait_restart;
+    ThreadPool::thread_id_t d_wait_subphase;
 
     char rankString[20];
     sprintf(rankString,"%05d",Dm->rank());
@@ -479,6 +503,12 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
     d_blobid_interval = db->getScalar<int>( "blobid_interval" );
     d_visualization_interval = db->getScalar<int>( "visualization_interval" );
     auto restart_file = db->getScalar<std::string>( "restart_file" );
+	if (db->keyExists( "subphase_analysis_interval" )){
+		d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
+	}
+	else{
+		d_subphase_analysis_interval = INT_MAX;
+	}
     d_restartFile = restart_file + "." + rankString;
     d_rank = MPI_WORLD_RANK();
     writeIDMap(ID_map_struct(),0,id_map_filename);
@@ -892,10 +922,17 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
     // if ( matches(type,AnalysisType::ComputeAverages) ) {
     if ( timestep%d_analysis_interval == 0 ) {
         auto work = new BasicWorkItem(type,timestep,Averages);
-        work->add_dependency(d_wait_analysis);    // Make sure we are done using analysis before modifying
+        work->add_dependency(d_wait_subphase);    // Make sure we are done using analysis before modifying
+        work->add_dependency(d_wait_analysis);    
         d_wait_analysis = d_tpool.add_work(work);
     }
     
+    if ( timestep%d_subphase_analysis_interval == 0 ) {
+        auto work = new BasicWorkItem(type,timestep,Averages);
+        work->add_dependency(d_wait_subphase);    // Make sure we are done using analysis before modifying
+        d_wait_subphase = d_tpool.add_work(work);
+    }
+
     if (timestep%d_restart_interval==0){
     	std::shared_ptr<double> cfq,cDen;
     	// Copy restart data to the CPU

analysis/runAnalysis.h

@@ -7,7 +7,7 @@
 #include "common/Communication.h"
 #include "common/ScaLBL.h"
 #include "threadpool/thread_pool.h"
-
+#include <limits.h>
 
 typedef std::shared_ptr<std::pair<int,IntArray>> BlobIDstruct;
 typedef std::shared_ptr<std::vector<BlobIDType>> BlobIDList;
@@ -88,6 +88,7 @@ private:
     int d_Np;
     int d_rank;
     int d_restart_interval, d_analysis_interval, d_blobid_interval, d_visualization_interval;
+    int d_subphase_analysis_interval;
     double d_beta;
     bool d_regular;
     ThreadPool d_tpool;
@@ -107,6 +108,7 @@ private:
     // Ids of work items to use for dependencies
     ThreadPool::thread_id_t d_wait_blobID;
     ThreadPool::thread_id_t d_wait_analysis;
+    ThreadPool::thread_id_t d_wait_subphase;
     ThreadPool::thread_id_t d_wait_vis;
     ThreadPool::thread_id_t d_wait_restart;
 

common/Domain.cpp

@@ -77,6 +77,7 @@ Domain::Domain( std::shared_ptr<Database> db, MPI_Comm Communicator):
 	Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0),
 	Comm(MPI_COMM_NULL),
 	inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
+	outlet_layers_x(0), outlet_layers_y(0), outlet_layers_z(0),
 	sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0),
 	sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0),
 	sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0),
@@ -126,6 +127,12 @@ void Domain::initialize( std::shared_ptr<Database> db )
 		inlet_layers_y = InletCount[1];
 		inlet_layers_z = InletCount[2];
 	}
+	if (d_db->keyExists( "OutletLayers" )){
+		auto OutletCount = d_db->getVector<int>( "OutletLayers" );
+		outlet_layers_x = OutletCount[0];
+		outlet_layers_y = OutletCount[1];
+		outlet_layers_z = OutletCount[2];
+	}
 	
     ASSERT( n.size() == 3u );
     ASSERT( L.size() == 3u );
@@ -144,9 +151,13 @@ void Domain::initialize( std::shared_ptr<Database> db )
     MPI_Comm_rank( Comm, &myrank );
 	rank_info = RankInfoStruct(myrank,nproc[0],nproc[1],nproc[2]);
 	// inlet layers only apply to lower part of domain
-	if (nproc[0] > 0) inlet_layers_x = 0;
+	if (rank_info.ix > 0) inlet_layers_x = 0;
-	if (nproc[1] > 0) inlet_layers_y = 0;
+	if (rank_info.jy > 0) inlet_layers_y = 0;
-	if (nproc[2] > 0) inlet_layers_z = 0;
+	if (rank_info.kz > 0) inlet_layers_z = 0;
+	// outlet layers only apply to top part of domain
+	if (rank_info.ix < nproc[0]-1 ) outlet_layers_x = 0;
+	if (rank_info.jy < nproc[1]-1) outlet_layers_y = 0;
+	if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0;
     // Fill remaining variables
 	N = Nx*Ny*Nz;
 	Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
@@ -531,6 +542,7 @@ void Domain::ReadIDs(){
 	readID=fread(id,1,N,IDFILE);
 	if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank());
 	fclose(IDFILE);
+
 	// Compute the porosity
 	double sum;
 	double porosity;
@@ -540,7 +552,8 @@ void Domain::ReadIDs(){
 	//.........................................................
 	// If external boundary conditions are applied remove solid
 	if (BoundaryCondition >  0  && kproc() == 0){
- 		for (int k=0; k<3; k++){
+    	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++){
 				for (int i=0;i<Nx;i++){
 					int n = k*Nx*Ny+j*Nx+i;
@@ -550,7 +563,8 @@ void Domain::ReadIDs(){
  		}
  	}
     if (BoundaryCondition >  0  && kproc() == nprocz()-1){
- 		for (int k=Nz-3; k<Nz; k++){
+    	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++){
  				for (int i=0;i<Nx;i++){
  					int n = k*Nx*Ny+j*Nx+i;
@@ -559,7 +573,7 @@ void Domain::ReadIDs(){
  			}
  		}
  	}
-    for (int k=1;k<Nz-1;k++){
+    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;

common/Domain.h

@@ -110,6 +110,7 @@ public: // Public variables (need to create accessors instead)
     double Lx,Ly,Lz,Volume;
     int Nx,Ny,Nz,N;
     int inlet_layers_x, inlet_layers_y, inlet_layers_z;
+    int outlet_layers_x, outlet_layers_y, outlet_layers_z;
     RankInfoStruct rank_info;
 
     MPI_Comm Comm;        // MPI Communicator for this domain