From 0b2f116fb034c1ab1e4c530e7e54b102e3c8d605 Mon Sep 17 00:00:00 2001
From: Rex Zhe Li <zhe.rex.li@gmail.com>
Date: Wed, 20 Dec 2017 18:17:00 +1100
Subject: [PATCH] fix the inlet cross-sectional area calculations

---
 tests/lbpm_color_macro_simulator.cpp | 56 ++++++++++------------------
 1 file changed, 19 insertions(+), 37 deletions(-)

diff --git a/tests/lbpm_color_macro_simulator.cpp b/tests/lbpm_color_macro_simulator.cpp
index fbdc93f8..defac936 100644
--- a/tests/lbpm_color_macro_simulator.cpp
+++ b/tests/lbpm_color_macro_simulator.cpp
@@ -370,7 +370,6 @@ int main(int argc, char **argv)
 //	char value;
 	char *id;
 	id = new char[N];
-	int sum = 0;
 	double sum_local;
 	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));
@@ -387,41 +386,6 @@ int main(int argc, char **argv)
 	ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N);
 	MPI_Barrier(comm);
 	if (rank == 0) cout << "Domain set." << endl;
-	
-	//.......................................................................
-	// Assign the phase ID field based on the signed distance
-	//.......................................................................
-	double inlet_area_local=0.f;
-	double InletArea;
-	for (k=0;k<Nz;k++){
-		for (j=0;j<Ny;j++){
-			for (i=0;i<Nx;i++){
-				int n = k*Nx*Ny+j*Nx+i;
-				id[n] = 0;
-			}
-		}
-	}
-	sum=0;
-	pore_vol = 0.0;
-	for ( k=0;k<Nz;k++){
-		for ( j=0;j<Ny;j++){
-			for ( i=0;i<Nx;i++){
-				int n = k*Nx*Ny+j*Nx+i;
-				if (Averages->SDs(n) > 0.0){
-					id[n] = 2;	
-				}
-				// compute the porosity (actual interface location used)
-				if (Averages->SDs(n) > 0.0){
-					sum++;	
-					if (k==1){
-					  inlet_area_local+=1.f;
-					}
-
-				}
-			}
-		}
-	}
-	MPI_Allreduce(&inlet_area_local,&InletArea,1,MPI_DOUBLE,MPI_SUM,comm);
 
 	if (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
 	sprintf(LocalRankFilename,"ID.%05i",rank);
@@ -430,9 +394,27 @@ int main(int argc, char **argv)
 	if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
 	readID=fread(id,1,N,IDFILE);
 	if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank);
-
 	fclose(IDFILE);
 
+    // Calculate the actual inlet cross-sectional area
+	double inlet_area_local=0.f;
+	double InletArea;
+    if (Dm.kproc==0)
+    {
+        for ( k=1;k<2;k++){
+            for ( j=1;j<Ny-1;j++){
+                for ( i=1;i<Nx-1;i++){
+                    int n = k*Nx*Ny+j*Nx+i;
+                    if (id[n] > 0){
+                        inlet_area_local+=1.f;
+                    }
+                }
+            }
+        }
+    
+    }
+	MPI_Allreduce(&inlet_area_local,&InletArea,1,MPI_DOUBLE,MPI_SUM,comm);
+
 	// Set up kstart, kfinish so that the reservoirs are excluded from averaging
 	int kstart,kfinish;
 	kstart = 1;