From ebabdcaef0c47b1a2fc01bfe8e58cd499f663e78 Mon Sep 17 00:00:00 2001
From: James E McClure <jemcclur@gmail.com>
Date: Tue, 26 Mar 2019 10:26:02 -0400
Subject: [PATCH] fix interface averaging for mass /momentum

---
 analysis/SubPhase.cpp     | 15 +++++++++++++--
 tests/TestColorBubble.cpp | 34 ++++++++++++++++++++--------------
 2 files changed, 33 insertions(+), 16 deletions(-)

diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index e59d91ef..e24ffe4f 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -290,7 +290,7 @@ void SubPhase::Full(){
 	if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
 	if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z;
 		
-	nd.reset();	nc.reset(); wd.reset();	wc.reset();	iwn.reset();
+	nd.reset();	nc.reset(); wd.reset();	wc.reset();	iwn.reset();	iwnc.reset();
 
  	Dm->CommunicateMeshHalo(Phi);
 	for (int k=1; k<Nz-1; k++){
@@ -543,7 +543,18 @@ void SubPhase::Full(){
 	gwc.Pz=sumReduce( Dm->Comm, wc.Pz);
 	gwc.K=sumReduce( Dm->Comm, wc.K);
 	gwc.p=sumReduce( Dm->Comm, wc.p);
-
+	
+	giwn.Mn=sumReduce( Dm->Comm, iwn.Mn);
+	giwn.Pnx=sumReduce( Dm->Comm, iwn.Pnx);
+	giwn.Pny=sumReduce( Dm->Comm, iwn.Pny);
+	giwn.Pnz=sumReduce( Dm->Comm, iwn.Pnz);
+	giwn.Kn=sumReduce( Dm->Comm, iwn.Kn);
+	giwn.Mw=sumReduce( Dm->Comm, iwn.Mw);
+	giwn.Pwx=sumReduce( Dm->Comm, iwn.Pwx);
+	giwn.Pwy=sumReduce( Dm->Comm, iwn.Pwy);
+	giwn.Pwz=sumReduce( Dm->Comm, iwn.Pwz);
+	giwn.Kw=sumReduce( Dm->Comm, iwn.Kw);
+	
 	// pressure averaging
 	if (vol_wc_bulk > 0.0)
 		wc.p = wc.p /vol_wc_bulk;
diff --git a/tests/TestColorBubble.cpp b/tests/TestColorBubble.cpp
index 10f343b1..0e6ea25a 100644
--- a/tests/TestColorBubble.cpp
+++ b/tests/TestColorBubble.cpp
@@ -14,13 +14,13 @@ using namespace std;
 inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
 	// initialize a bubble
 	int i,j,k,n;
-	int rank = ColorModel.Mask->rank();
-	int nprocx = ColorModel.Mask->nprocx();
-	int nprocy = ColorModel.Mask->nprocy();
-	int nprocz = ColorModel.Mask->nprocz();
-	int Nx = ColorModel.Mask->Nx;
-	int Ny = ColorModel.Mask->Ny;
-	int Nz = ColorModel.Mask->Nz;
+	int rank = ColorModel.Dm->rank();
+	int nprocx = ColorModel.Dm->nprocx();
+	int nprocy = ColorModel.Dm->nprocy();
+	int nprocz = ColorModel.Dm->nprocz();
+	int Nx = ColorModel.Dm->Nx;
+	int Ny = ColorModel.Dm->Ny;
+	int Nz = ColorModel.Dm->Nz;
 	if (rank == 0) cout << "Setting up bubble..." << endl;
 	for (k=0;k<Nz;k++){
 		for (j=0;j<Ny;j++){
@@ -34,14 +34,12 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 	for (k=0;k<Nz;k++){
 		for (j=0;j<Ny;j++){
 			for (i=0;i<Nx;i++){
-				n = k*Nx*Ny + j*Nz + i;
-				int iglobal= i+(Nx-2)*ColorModel.Mask->iproc();
-				int jglobal= j+(Ny-2)*ColorModel.Mask->jproc();
-				int kglobal= k+(Nz-2)*ColorModel.Mask->kproc();
+				n = k*Nx*Ny + j*Nx + i;
+				double iglobal= double(i+(Nx-2)*ColorModel.Dm->iproc())-double((Nx-2)*nprocx)*0.5;
+				double jglobal= double(j+(Ny-2)*ColorModel.Dm->jproc())-double((Ny-2)*nprocy)*0.5;
+				double kglobal= double(k+(Nz-2)*ColorModel.Dm->kproc())-double((Nz-2)*nprocz)*0.5;
 				// Initialize phase position field for parallel bubble test
-				if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
-						+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
-						+(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
+				if ((iglobal*iglobal)+(jglobal*jglobal)+(kglobal*kglobal) < BubbleRadius*BubbleRadius){
 					ColorModel.Mask->id[n] = 2;
 					ColorModel.Mask->id[n] = 2;
 				}
@@ -50,9 +48,17 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
 					ColorModel.Mask->id[n]=1;
 				}
 				ColorModel.id[n] = ColorModel.Mask->id[n];
+				ColorModel.Dm->id[n] = ColorModel.Mask->id[n];
 			}
 		}
 	}
+	
+	FILE *OUTFILE;
+	char LocalRankFilename[40];
+	sprintf(LocalRankFilename,"Bubble.%05i.raw",rank);
+	OUTFILE = fopen(LocalRankFilename,"wb");
+	fwrite(ColorModel.id,1,Nx*Ny*Nz,OUTFILE);
+	fclose(OUTFILE);
 	// initialize the phase indicator field
 }
 //***************************************************************************************