From 43f1d53fa7aefdb46799b79c732b0638c4e1b2d5 Mon Sep 17 00:00:00 2001
From: James E McClure <mcclurej@vt.edu>
Date: Fri, 18 Mar 2016 10:25:06 -0400
Subject: [PATCH] Adding domain decomposition to lbpm_uCT_pp (in progress)

---
 tests/lbpm_segmented_decomp.cpp |   2 +-
 tests/lbpm_uCT_pp.cpp           | 260 ++++++++++++++++++++++++--------
 2 files changed, 198 insertions(+), 64 deletions(-)

diff --git a/tests/lbpm_segmented_decomp.cpp b/tests/lbpm_segmented_decomp.cpp
index a54e8a1c..24ec27d4 100644
--- a/tests/lbpm_segmented_decomp.cpp
+++ b/tests/lbpm_segmented_decomp.cpp
@@ -23,7 +23,7 @@ int main(int argc, char **argv)
 	MPI_Comm_size(comm,&nprocs);
 
 
-        int SOLID=atoi(argv[1]);
+    int SOLID=atoi(argv[1]);
 	int NWP=atoi(argv[2]);
 	if (rank==0){
 	  printf("Solid Label %i \n",SOLID);
diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp
index bfe929d8..08128ae9 100644
--- a/tests/lbpm_uCT_pp.cpp
+++ b/tests/lbpm_uCT_pp.cpp
@@ -3,81 +3,215 @@
 // and averaging on a blob-by-blob basis
 // James E. McClure 2014
 
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "common/Array.h"
+#include "common/Domain.h"
 #include "IO/netcdf.h"
 
 #include "ProfilerApp.h"
 
-
-void load( const std::string filename )
-{
-    int fid = netcdf::open( filename );
-
-    std::vector<std::string> vars = netcdf::getVarNames( fid );
-    for (size_t i=0; i<vars.size(); i++) {
-        printf("Reading variable %s\n",vars[i].c_str());
-        netcdf::VariableType type = netcdf::getVarType( fid, vars[i] );
-        if ( type == netcdf::STRING )
-            Array<std::string> tmp = netcdf::getVar<std::string>( fid, vars[i] );
-        else if ( type == netcdf::SHORT )
-            Array<short> tmp = netcdf::getVar<short>( fid, vars[i] );
-        else
-            Array<double> tmp = netcdf::getVar<double>( fid, vars[i] );
-    }
-
-    std::vector<std::string> attr = netcdf::getAttNames( fid );
-    for (size_t i=0; i<attr.size(); i++) {
-        printf("Reading attribute %s\n",attr[i].c_str());
-        netcdf::VariableType type = netcdf::getAttType( fid, attr[i] );
-        if ( type == netcdf::STRING )
-            Array<std::string> tmp = netcdf::getAtt<std::string>( fid, attr[i] );
-        else
-            Array<double> tmp = netcdf::getAtt<double>( fid, attr[i] );
-    }
-    netcdf::close( fid );
-}
-
-
 int main(int argc, char **argv)
 {
-    PROFILE_START("Main");
 
-    //std::vector<std::string> filenames;
-    std::string filename;
-    if ( argc==0 ) {
-        printf("At least one filename must be specified\n");
-        return 1;
-    }
-    else {
-      filename=std::string(argv[1]);
-      printf("Reading file: %s\n",filename.c_str());
-    }
-						
+	// 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);
 
-   int fid = netcdf::open(filename);
-
-    // Read all of the attributes
-    std::vector<std::string> attr = netcdf::getAttNames( fid );
-    for (size_t i=0; i<attr.size(); i++) {
-        printf("Reading attribute %s\n",attr[i].c_str());
-        netcdf::VariableType type = netcdf::getAttType( fid, attr[i] );
-        if ( type == netcdf::STRING ){
-	  Array<std::string> tmp = netcdf::getAtt<std::string>( fid, attr[i] );
+	//std::vector<std::string> filenames;
+	std::string filename;
+	if ( argc==0 ) {
+		printf("At least one filename must be specified\n");
+		return 1;
 	}
-        else{
-	  //Array<double> tmp = netcdf::getAtt<double>( fid, attr[i] );
+	else {
+		filename=std::string(argv[1]);
+		printf("Input data file: %s\n",filename.c_str());
 	}
+
+    //.......................................................................
+    // Reading the domain information file
+    //.......................................................................
+    int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
+    double Lx, Ly, Lz;
+    int Nx,Ny,Nz;
+    int i,j,k,n;
+	int BC=0;
+
+    if (rank==0){
+    	ifstream domain("Domain.in");
+    	domain >> nprocx;
+    	domain >> nprocy;
+    	domain >> nprocz;
+    	domain >> nx;
+    	domain >> ny;
+    	domain >> nz;
+    	domain >> nspheres;
+    	domain >> Lx;
+    	domain >> Ly;
+    	domain >> Lz;
+
+    }
+	MPI_Barrier(comm);
+	// Computational domain
+	//.................................................
+	MPI_Bcast(&nx,1,MPI_INT,0,comm);
+	MPI_Bcast(&ny,1,MPI_INT,0,comm);
+	MPI_Bcast(&nz,1,MPI_INT,0,comm);
+	MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
+	MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
+	MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
+	MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
+	MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
+	MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
+/*	//.................................................
+	MPI_Bcast(&Ny,1,MPI_INT,0,comm);
+	MPI_Bcast(&Ny,1,MPI_INT,0,comm);
+	MPI_Bcast(&Nz,1,MPI_INT,0,comm);
+	MPI_Bcast(&xStart,1,MPI_INT,0,comm);
+	MPI_Bcast(&yStart,1,MPI_INT,0,comm);
+	MPI_Bcast(&zStart,1,MPI_INT,0,comm);
+*/	//.................................................
+	MPI_Barrier(comm);
+
+    // Check that the number of processors >= the number of ranks
+    if ( rank==0 ) {
+        printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
+        printf("Number of MPI ranks used: %i \n", nprocs);
+        printf("Full domain size: %i x %i x %i  \n",nx*nprocx,ny*nprocy,nz*nprocz);
+    }
+    if ( nprocs < nprocx*nprocy*nprocz ){
+        ERROR("Insufficient number of processors");
     }
 
-    // Read the data array
-    std::string varname("VOLUME");
-    printf("Reading %s\n",varname.c_str());
-    Array<float> VOLUME = netcdf::getVar<float>( fid, varname);
-    //printf("VOLUME size = %zu \n",VOLUME.length());	
-    printf("VOLUME dims =  %zu x %zu x %zu \n",VOLUME.size(0),VOLUME.size(1),VOLUME.size(2));	
-    printf("Sucess!! \n");
-    netcdf::close( fid );
 
-    PROFILE_SAVE("TestNetcdf");
-    return 0;  
+	PROFILE_START("ReadVolume");
+
+	// Read the input volume to rank 0 only, then distribute pieces to workers
+	if (rank==0){
+		// Open the netcdf file
+		int fid = netcdf::open(filename);
+
+		// Read all of the attributes
+		std::vector<std::string> attr = netcdf::getAttNames( fid );
+		for (size_t i=0; i<attr.size(); i++) {
+			printf("Reading attribute %s\n",attr[i].c_str());
+			netcdf::VariableType type = netcdf::getAttType( fid, attr[i] );
+			if ( type == netcdf::STRING ){
+				Array<std::string> tmp = netcdf::getAtt<std::string>( fid, attr[i] );
+			}
+			else{
+				//Array<double> tmp = netcdf::getAtt<double>( fid, attr[i] );
+			}
+		}
+
+		// Read the VOLUME data array
+		std::string varname("VOLUME");
+		printf("Reading %s\n",varname.c_str());
+		Array<float> VOLUME = netcdf::getVar<float>( fid, varname);
+		printf("VOLUME dims =  %zu x %zu x %zu \n",VOLUME.size(0),VOLUME.size(1),VOLUME.size(2));
+		printf("Sucess!! \n");
+		netcdf::close( fid );
+	}
+	PROFILE_SAVE("ReadVolume");
+
+    MPI_Barrier(comm);
+
+    // Get the rank info
+    int N = (nx+2)*(ny+2)*(nz+2);
+	Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
+	for (k=0;k<nz+2;k++){
+		for (j=0;j<ny+2;j++){
+			for (i=0;i<nx+2;i++){
+				n = k*(nx+2)*(ny+2)+j*(nx+2)+i;
+				Dm.id[n] = 1;
+			}
+		}
+	}
+	Dm.CommInit(comm);
+/*
+	// Set up the sub-domains
+	if (rank==0){
+		printf("Distributing subdomains across %i processors \n",nprocs);
+		printf("Process grid: %i x %i x %i \n",Dm.nprocx,Dm.nprocy,Dm.nprocz);
+		printf("Subdomain size: %i \n",N);
+		printf("Size of transition region: %i \n", z_transition_size);
+		char *tmp;
+		tmp = new char[N];
+		for (int kp=0; kp<nprocz; kp++){
+			for (int jp=0; jp<nprocy; jp++){
+				for (int ip=0; ip<nprocx; ip++){
+					// rank of the process that gets this subdomain
+					int rnk = kp*Dm.nprocx*Dm.nprocy + jp*Dm.nprocx + ip;
+					// Pack and send the subdomain for rnk
+					for (k=0;k<nz+2;k++){
+						for (j=0;j<ny+2;j++){
+							for (i=0;i<nx+2;i++){
+								int x = xStart + ip*nx + i-1;
+								int y = yStart + jp*ny + j-1;
+						//		int z = zStart + kp*nz + k-1;
+								int z = zStart + kp*nz + k-1 - z_transition_size;
+								if (z<zStart) 	z=zStart;
+								if (!(z<Nz))	z=Nz-1;
+								int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
+								int nglobal = z*Nx*Ny+y*Nx+x;
+								tmp[nlocal] = SegData[nglobal];
+							}
+						}
+					}
+					if (rnk==0){
+						for (k=0;k<nz+2;k++){
+							for (j=0;j<ny+2;j++){
+								for (i=0;i<nx+2;i++){
+									int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
+									Dm.id[nlocal] = tmp[nlocal];
+								}
+							}
+						}
+					}
+					else{
+						printf("Sending data to process %i \n", rnk);
+						MPI_Send(tmp,N,MPI_CHAR,rnk,15,comm);
+					}
+				}
+			}
+		}
+	}
+	else{
+		// Recieve the subdomain from rank = 0
+		printf("Ready to recieve data %i at process %i \n", N,rank);
+		MPI_Recv(Dm.id,N,MPI_CHAR,0,15,comm,MPI_STATUS_IGNORE);
+	}
+	MPI_Barrier(comm);
+
+	nx+=2; ny+=2; nz+=2;
+	N=nx*ny*nz;
+
+	if (rank==0) printf("All sub-domains recieved \n");
+    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;
+			        if (Dm.id[n]==char(SOLID))     Dm.id[n] = 0;
+			       	else if (Dm.id[n]==char(NWP))  Dm.id[n] = 1;
+			       	else                           Dm.id[n] = 2;
+
+			}
+		}
+	}
+	if (rank==0) printf("Domain set \n");
+*/
+
+    MPI_Barrier(comm);
+    MPI_Finalize();
+	return 0;
 }