From c76e549cabd361900b8b41f83e515c6cd968194a Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Fri, 30 Nov 2018 15:53:28 -0500
Subject: [PATCH] compute volume fractions for components

---
 models/ColorModel.cpp | 33 ++++++++++++++++++++++-----------
 1 file changed, 22 insertions(+), 11 deletions(-)

diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 85d60891..a72cd266 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -169,15 +169,12 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 		ERROR("Error: ComponentLabels and ComponentAffinity must be the same length! \n");
 	}
 
-	if (rank==0){
-		printf("Components labels: %lu \n",NLABELS);
-		for (unsigned int idx=0; idx<NLABELS; idx++){
-			VALUE=LabelList[idx];
-			AFFINITY=AffinityList[idx];
-			printf("   label=%i, affinity=%f\n",int(VALUE),AFFINITY); 
-		}
-	}
+	int label_count[NLABELS];
+	int label_count_global[NLABELS];
 	// Assign the labels
+
+	for (int idx=0; idx<NLABELS; idx++) label_count[idx]=0;
+
 	for (int k=0;k<Nz;k++){
 		for (int j=0;j<Ny;j++){
 			for (int i=0;i<Nx;i++){
@@ -185,9 +182,10 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
 				VALUE=id[n];
 				// Assign the affinity from the paired list
 				for (unsigned int idx=0; idx < NLABELS; idx++){
-					//printf("rank=%i, idx=%i, value=%i, %i, \n",rank(),idx, VALUE,LabelList[idx]);
+				      //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
 					if (VALUE == LabelList[idx]){
 						AFFINITY=AffinityList[idx];
+						label_count[idx]++;
 						idx = NLABELS;
 						Mask->id[n] = 0; // set mask to zero since this is an immobile component
 					}
@@ -201,6 +199,19 @@ 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]; 
+
+	MPI_Allreduce(&label_count[0],&label_count_global[0],NLABELS,MPI_INT,MPI_SUM,Dm->Comm);
+
+	if (rank==0){
+		printf("Components labels: %lu \n",NLABELS);
+		for (unsigned int idx=0; idx<NLABELS; idx++){
+			VALUE=LabelList[idx];
+			AFFINITY=AffinityList[idx];
+			double volume_fraction  = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocz);
+			printf("   label=%i, affinity=%f, volume fraction==%f\n",int(VALUE),AFFINITY,volume_fraction); 
+		}
+	}
+
 }
 
 
@@ -577,8 +588,8 @@ void ScaLBL_ColorModel::Run(){
 					if (rank==0){
 						printf("** WRITE STEADY POINT *** ");
 						printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
-						volA /= double(Nx*Ny*Nz*nprocs);
-						volB /= double(Nx*Ny*Nz*nprocs);
+						volA /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocz);
+						volB /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocz);
 						FILE * kr_log_file = fopen("relperm.csv","a");
 						fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g ",timestep-analysis_interval+20,muA,muB,5.796*alpha,Fx,Fy,Fz);
 						fprintf(kr_log_file,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z);