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

analysis/Minkowski.cpp

@@ -82,6 +82,22 @@ void Minkowski::Initialize()
 	vol_n = euler = Jn = An = Kn = 0.0;
 }
 
+double Minkowski::V(){
+	return vol_n_global;
+}
+
+double Minkowski::A(){
+	return An_global;
+}
+
+double Minkowski::J(){
+	return Jn_global;
+}
+
+double Minkowski::X(){
+	return euler_global;
+}
+
 
 void Minkowski::UpdateMeshValues()
 {
@@ -148,11 +164,9 @@ void Minkowski::ComputeLocal()
 				// Compute volume averages
 				for (int p=0;p<8;p++){
 					n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
-					if ( Dm->id[n] != 0 ){
-						// 1-D index for this cube corner
-						if ( SDn(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < 0 ){
-							vol_n += 0.125;
-						}
+					// 1-D index for this cube corner
+					if ( SDn(i+cube[p][0],j+cube[p][1],k+cube[p][2]) < 0 ){
+						vol_n += 0.125;
 					}
 				}
 

analysis/Minkowski.h

@@ -84,6 +84,11 @@ public:
 	DoubleArray SDn_y;
 	DoubleArray SDn_z;
 
+	double V();
+	double A();
+	double J();
+	double X();
+	
 	//...........................................................................
 	Minkowski(std::shared_ptr <Domain> Dm);
 	~Minkowski();
@@ -95,6 +100,7 @@ public:
 	void PrintAll();
 	int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel);
 	void SortBlobs();
+	
 };
 
 

analysis/TwoPhase.cpp

@@ -424,7 +424,7 @@ void TwoPhase::UpdateMeshValues()
 		for (j=0; j<Ny; j++){
 			for (i=0; i<Nx; i++){
 				n = k*Nx*Ny+j*Nx+i;
-				if (Dm->id[n] == 0){
+				if (!(Dm->id[n] > 0)){
 					// Solid phase
 					PhaseID(i,j,k) = 0;
 				}
@@ -461,7 +461,7 @@ void TwoPhase::ComputeLocal()
 				// Compute volume averages
 				for (int p=0;p<8;p++){
 					n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
-					if ( Dm->id[n] != 0 ){
+					if ( Dm->id[n] > 0 ){
 						// 1-D index for this cube corner
 						// compute the norm of the gradient of the phase indicator field
 						// Compute the non-wetting phase volume contribution
@@ -654,7 +654,7 @@ void TwoPhase::ComponentAverages()
 				// Compute volume averages
 				for (int p=0;p<8;p++){
 					n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
-					if ( Dm->id[n] != 0 ){
+					if ( Dm->id[n] > 0 ){
 						// 1-D index for this cube corner
 						// compute the norm of the gradient of the phase indicator field
 						// Compute the non-wetting phase volume contribution

analysis/distance.cpp

@@ -168,7 +168,7 @@ void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm,
     for (int i=0; i<N[0]; i++) {
         for (int j=0; j<N[1]; j++) {
             ID(i,j,0) = ID(i,j,1);
-            ID(i,j,N[1]-1) = ID(i,j,N[2]-2);
+            ID(i,j,N[2]-1) = ID(i,j,N[2]-2);
         }
     }
     // Communicate ghosts

analysis/runAnalysis.cpp

@@ -28,7 +28,6 @@
 #include "ProfilerApp.h"
 
 
-
 AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs)  
 {
     lhs = static_cast<AnalysisType>(
@@ -178,15 +177,24 @@ public:
         PROFILE_START("Save Vis",1);
         ASSERT(visData[0].vars[0]->name=="phase");
         ASSERT(visData[0].vars[1]->name=="Pressure");
-        ASSERT(visData[0].vars[2]->name=="SignDist");
-        ASSERT(visData[0].vars[3]->name=="BlobID");
+        ASSERT(visData[0].vars[2]->name=="Velocity_x");
+        ASSERT(visData[0].vars[3]->name=="Velocity_y");
+        ASSERT(visData[0].vars[4]->name=="Velocity_z");
+        ASSERT(visData[0].vars[5]->name=="SignDist");
+        ASSERT(visData[0].vars[6]->name=="BlobID");
         Array<double>& PhaseData = visData[0].vars[0]->data;
         Array<double>& PressData = visData[0].vars[1]->data;
-        Array<double>& SignData  = visData[0].vars[2]->data;
-        Array<double>& BlobData  = visData[0].vars[3]->data;
+        Array<double>& VelxData = visData[0].vars[2]->data;
+        Array<double>& VelyData = visData[0].vars[3]->data;
+        Array<double>& VelzData = visData[0].vars[4]->data;
+        Array<double>& SignData  = visData[0].vars[5]->data;
+        Array<double>& BlobData  = visData[0].vars[6]->data;
         fillData.copy(Averages.SDn,PhaseData);
         fillData.copy(Averages.Press,PressData);
         fillData.copy(Averages.SDs,SignData);
+        fillData.copy(Averages.Vel_x,VelxData);
+        fillData.copy(Averages.Vel_y,VelyData);
+        fillData.copy(Averages.Vel_z,VelzData);
         fillData.copy(Averages.Label_NWP,BlobData);
         IO::writeData( timestep, visData, comm.comm );
         PROFILE_STOP("Save Vis",1);
@@ -330,6 +338,11 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
 	d_meshData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
 	auto PhaseVar = std::make_shared<IO::Variable>();
 	auto PressVar = std::make_shared<IO::Variable>();
+
+	auto VxVar = std::make_shared<IO::Variable>();
+	auto VyVar = std::make_shared<IO::Variable>();
+	auto VzVar = std::make_shared<IO::Variable>();
+
 	auto SignDistVar = std::make_shared<IO::Variable>();
 	auto BlobIDVar = std::make_shared<IO::Variable>();
 	PhaseVar->name = "phase";
@@ -342,6 +355,23 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
 	PressVar->dim = 1;
 	PressVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
 	d_meshData[0].vars.push_back(PressVar);
+	
+	VxVar->name = "Velocity_x";
+	VxVar->type = IO::VariableType::VolumeVariable;
+	VxVar->dim = 1;
+	VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	d_meshData[0].vars.push_back(VxVar);
+	VyVar->name = "Velocity_y";
+	VyVar->type = IO::VariableType::VolumeVariable;
+	VyVar->dim = 1;
+	VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	d_meshData[0].vars.push_back(VyVar);
+	VzVar->name = "Velocity_z";
+	VzVar->type = IO::VariableType::VolumeVariable;
+	VzVar->dim = 1;
+	VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
+	d_meshData[0].vars.push_back(VzVar);
+	
 	SignDistVar->name = "SignDist";
 	SignDistVar->type = IO::VariableType::VolumeVariable;
 	SignDistVar->dim = 1;
@@ -564,35 +594,35 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
     }
     //if ( matches(type,AnalysisType::CopySimState) ) {
     if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
-       // Copy the members of Averages to the cpu (phase was copied above)
-        PROFILE_START("Copy-Pressure",1);
-        ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
-        ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
-        ScaLBL_DeviceBarrier();
-        PROFILE_STOP("Copy-Pressure",1);
-        PROFILE_START("Copy-Wait",1);
-        PROFILE_STOP("Copy-Wait",1);
-        PROFILE_START("Copy-State",1);
-        //memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
-	if (d_regular)
-	  d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
-	else
-	  ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
-	// copy other variables
-        d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
-        d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
-        d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
-        d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
-        PROFILE_STOP("Copy-State",1);
+    	// Copy the members of Averages to the cpu (phase was copied above)
+    	PROFILE_START("Copy-Pressure",1);
+    	ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
+    	//ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
+    	ScaLBL_DeviceBarrier();
+    	PROFILE_STOP("Copy-Pressure",1);
+    	PROFILE_START("Copy-Wait",1);
+    	PROFILE_STOP("Copy-Wait",1);
+    	PROFILE_START("Copy-State",1);
+    	//memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
+    	if (d_regular)
+    		d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
+    	else
+    		ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
+    	// copy other variables
+    	d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
+    	d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
+    	d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
+    	d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
+    	PROFILE_STOP("Copy-State",1);
     }
     std::shared_ptr<double> cfq,cPhi;
     //if ( matches(type,AnalysisType::CreateRestart) ) {
     if (timestep%d_restart_interval==0){
-        // Copy restart data to the CPU
-        cPhi = std::shared_ptr<double>(new double[d_Np],DeleteArray<double>);
-        cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
-        ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
-        ScaLBL_CopyToHost(cPhi.get(),Phi,d_Np*sizeof(double));
+    	// Copy restart data to the CPU
+    	cPhi = std::shared_ptr<double>(new double[d_Np],DeleteArray<double>);
+    	cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
+    	ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
+    	ScaLBL_CopyToHost(cPhi.get(),Phi,d_Np*sizeof(double));
     }
     PROFILE_STOP("Copy data to host",1);
 

common/ScaLBL.cpp

@@ -1136,14 +1136,12 @@ void ScaLBL_Communicator::RecvGrad(double *phi, double *grad){
 	ScaLBL_DeviceBarrier();
 
 	//...................................................................................
-	// NOTE: AA Routine writes to opposite 
 	// Unpack the gradributions on the device
 	//...................................................................................
 	//...Unpacking for x face(2,8,10,12,14)................................
 	ScaLBL_Gradient_Unpack(1.0,-1,0,0,dvcRecvDist_x,0,recvCount_x,recvbuf_x,phi,grad,N);
 	ScaLBL_Gradient_Unpack(0.5,-1,-1,0,dvcRecvDist_x,recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
 	ScaLBL_Gradient_Unpack(0.5,-1,1,0,dvcRecvDist_x,2*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
-	ScaLBL_Gradient_Unpack(0.5,-1,0,-1,dvcRecvDist_x,3*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
 	ScaLBL_Gradient_Unpack(0.5,-1,0,1,dvcRecvDist_x,4*recvCount_x,recvCount_x,recvbuf_x,phi,grad,N);
 	//...................................................................................
 	//...Packing for X face(1,7,9,11,13)................................
@@ -1546,7 +1544,7 @@ void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, Double
 	int Nz = map.size(2);
 
 	// initialize the array
-	regdata.fill(-1.f);
+	regdata.fill(0.f);
 	
 	double *TmpDat;
 	double value;

cpu/Color.cpp

@@ -1267,7 +1267,7 @@ extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, int
 //		double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 //		double Fx, double Fy, double Fz, int start, int finish, int Np){
 extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
-		double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
+		double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
 		double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
 
 	int ijk,nn,n;
@@ -1785,9 +1785,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do
 		ux = jx / rho0;
 		uy = jy / rho0;
 		uz = jz / rho0;
-		//Velocity[n] = ux;
-		//Velocity[Np+n] = uy;
-		//Velocity[2*Np+n] = uz;
+		Vel[n] = ux;
+		Vel[Np+n] = uy;
+		Vel[2*Np+n] = uz;
 
 		// Instantiate mass transport distributions
 		// Stationary value - distribution 0
@@ -2433,9 +2433,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *di
 		ux = jx / rho0;
 		uy = jy / rho0;
 		uz = jz / rho0;
-		//Velocity[n] = ux;
-		//Velocity[Np+n] = uy;
-		//Velocity[2*Np+n] = uz;
+		Vel[n] = ux;
+		Vel[Np+n] = uy;
+		Vel[2*Np+n] = uz;
 
 		// Instantiate mass transport distributions
 		// Stationary value - distribution 0

cpu/dfh.cpp

@@ -549,9 +549,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, double
 		//........................................................................
 
 		// write the velocity 
-		ux = jx / rho0;
-		uy = jy / rho0;
-		uz = jz / rho0;
+		ux = (jx + force_x) / rho0;
+		uy = (jy + force_y) / rho0;
+		uz = (jz + force_z) / rho0;
 		//Velocity[n] = ux;
 		//Velocity[Np+n] = uy;
 		//Velocity[2*Np+n] = uz;
@@ -1139,9 +1139,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, double *
 		dist[nread] = fq;
 
 		// write the velocity 
-		ux = jx / rho0;
-		uy = jy / rho0;
-		uz = jz / rho0;
+		ux = (jx + force_x) / rho0;
+		uy = (jy + force_y) / rho0;
+		uz = (jz + force_z) / rho0;
 		//Velocity[n] = ux;
 		//Velocity[Np+n] = uy;
 		//Velocity[2*Np+n] = uz;

example/Sph1896/input.db

@@ -1,3 +1,13 @@
+MRT {
+    tau = 1.0
+    F = 0, 0, 1e-4
+    timestepMax = 1000
+    din = 1.0
+    dout = 1.0
+    Restart = false
+    flux = 0.0
+}
+
 Color {
     tau   = 1.0;  
     alpha = 1e-2;
@@ -18,8 +28,8 @@ Color {
 
 Domain {
     nproc = 1, 1, 1     // Number of processors (Npx,Npy,Npz)
-    n = 80, 80, 80      // Size of local domain (Nx,Ny,Nz)
-    n_spheres = 1       // Number of spheres
+    n = 120, 120, 120      // Size of local domain (Nx,Ny,Nz)
+    nspheres = 1896       // Number of spheres
     L = 1, 1, 1         // Length of domain (x,y,z)
     BC = 0              // Boundary condition type
 }

gpu/Color.cu

@@ -1777,9 +1777,9 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A
 			ux = jx / rho0;
 			uy = jy / rho0;
 			uz = jz / rho0;
-			//Velocity[n] = ux;
-			//Velocity[Np+n] = uy;
-			//Velocity[2*Np+n] = uz;
+			Velocity[n] = ux;
+			Velocity[Np+n] = uy;
+			Velocity[2*Np+n] = uz;
 
 			// Instantiate mass transport distributions
 			// Stationary value - distribution 0
@@ -2426,9 +2426,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double
 			ux = jx / rho0;
 			uy = jy / rho0;
 			uz = jz / rho0;
-			//Velocity[n] = ux;
-			//Velocity[Np+n] = uy;
-			//Velocity[2*Np+n] = uz;
+			Velocity[n] = ux;
+			Velocity[Np+n] = uy;
+			Velocity[2*Np+n] = uz;
 
 			// Instantiate mass transport distributions
 			// Stationary value - distribution 0

gpu/dfh.cu

@@ -566,9 +566,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_DFH(int *neighborList, double *dist, dou
 			//........................................................................
 
 			// write the velocity 
-			ux = jx / rho0;
-			uy = jy / rho0;
-			uz = jz / rho0;
+			ux = (jx + force_x) / rho0;
+			uy = (jy + force_y) / rho0;
+			uz = (jz + force_z) / rho0;
 			//Velocity[n] = ux;
 			//Velocity[Np+n] = uy;
 			//Velocity[2*Np+n] = uz;
@@ -1159,9 +1159,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_DFH(int *neighborList, double *dist, doub
 			dist[nread] = fq;
 
 			// write the velocity 
-			ux = jx / rho0;
-			uy = jy / rho0;
-			uz = jz / rho0;
+			ux = (jx + force_x) / rho0;
+			uy = (jy + force_y) / rho0;
+			uz = (jz + force_z) / rho0;
 			//Velocity[n] = ux;
 			//Velocity[Np+n] = uy;
 			//Velocity[2*Np+n] = uz;

models/ColorModel.cpp

@@ -261,8 +261,8 @@ void ScaLBL_ColorModel::Create(){
 		}
 	}
 	// check that TmpMap is valid
-	for (int idx=0; idx<ScaLBL_Comm->last_interior; idx++){
-		if (idx == ScaLBL_Comm->next) idx = ScaLBL_Comm->first_interior;
+	for (int idx=0; idx<ScaLBL_Comm->LastInterior(); idx++){
+		if (idx == ScaLBL_Comm->LastExterior()) idx = ScaLBL_Comm->FirstInterior();
 		int n = TmpMap[idx];
 		if (n > Nx*Ny*Nz){
 			printf("Bad value! idx=%i \n");
@@ -331,8 +331,8 @@ void ScaLBL_ColorModel::Initialize(){
 	if (rank==0)	printf ("Initializing distributions \n");
 	ScaLBL_D3Q19_Init(fq, Np);
 	if (rank==0)	printf ("Initializing phase field \n");
-	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->next, Np);
-	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
+	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
+	ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 
 	if (BoundaryCondition >0 ){
 		if (Dm->kproc()==0){
@@ -380,9 +380,9 @@ void ScaLBL_ColorModel::Run(){
 		// Compute the Phase indicator field
 		// Read for Aq, Bq happens in this routine (requires communication)
 		ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
-		ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
+		ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
-		ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
+		ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
 		
 		// Perform the collision operation
 		ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
@@ -390,7 +390,7 @@ void ScaLBL_ColorModel::Run(){
 		ScaLBL_Comm_Regular->SendHalo(Phi);
 
 		ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
-				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
+				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm_Regular->RecvHalo(Phi);
 		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
 		// Set BCs
@@ -407,23 +407,23 @@ void ScaLBL_ColorModel::Run(){
 			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
 		}
 		ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
-				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->next, Np);
+				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
 		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 
 		// *************EVEN TIMESTEP*************
 		timestep++;
 		// Compute the Phase indicator field
 		ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
-		ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
+		ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
-		ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->next, Np);
+		ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
 
 		// Perform the collision operation
 		ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
 		// Halo exchange for phase field
 		ScaLBL_Comm_Regular->SendHalo(Phi);
 		ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
-				alpha, beta, Fx, Fy, Fz,  Nx, Nx*Ny, ScaLBL_Comm->first_interior, ScaLBL_Comm->last_interior, Np);
+				alpha, beta, Fx, Fy, Fz,  Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
 		ScaLBL_Comm_Regular->RecvHalo(Phi);
 		ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
 		// Set boundary conditions
@@ -440,7 +440,7 @@ void ScaLBL_ColorModel::Run(){
 			ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
 		}
 		ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
-				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->next, Np);
+				alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
 		ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
 		//************************************************************************
 		

models/DFHModel.cpp

@@ -606,4 +606,19 @@ void ScaLBL_DFHModel::WriteDebug(){
 	OUTFILE = fopen(LocalRankFilename,"wb");
 	fwrite(PhaseField.data(),8,N,OUTFILE);
 	fclose(OUTFILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
+	FILE *AFILE;
+	sprintf(LocalRankFilename,"A.%05i.raw",rank);
+	AFILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,AFILE);
+	fclose(AFILE);
+
+	ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
+	FILE *BFILE;
+	sprintf(LocalRankFilename,"B.%05i.raw",rank);
+	BFILE = fopen(LocalRankFilename,"wb");
+	fwrite(PhaseField.data(),8,N,BFILE);
+	fclose(BFILE);
+
 }

models/MRTModel.cpp

@@ -139,6 +139,7 @@ void ScaLBL_MRTModel::Create(){
 	// copy the neighbor list 
 	ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
 	MPI_Barrier(comm);
+	
 }        
 
 void ScaLBL_MRTModel::Initialize(){
@@ -177,21 +178,66 @@ void ScaLBL_MRTModel::Run(){
 	}
 	//************************************************************************/
 	stoptime = MPI_Wtime();
-	//	cout << "CPU time: " << (stoptime - starttime) << " seconds" << endl;
-	cputime = stoptime - starttime;
-	//	cout << "Lattice update rate: "<< double(Nx*Ny*Nz*timestep)/cputime/1000000 <<  " MLUPS" << endl;
-	double MLUPS = double(Np*timestep)/cputime/1000000;
-//	if (rank==0) printf("********************************************************\n");
-//	if (rank==0) printf("CPU time = %f \n", cputime);
-//	if (rank==0) printf("Lattice update rate (per process)= %f MLUPS \n", MLUPS);
+	if (rank==0) printf("-------------------------------------------------------------------\n");
+	// Compute the walltime per timestep
+	cputime = (stoptime - starttime)/timestep;
+	// Performance obtained from each node
+	double MLUPS = double(Np)/cputime/1000000;
+
+	if (rank==0) printf("********************************************************\n");
+	if (rank==0) printf("CPU time = %f \n", cputime);
+	if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
 	MLUPS *= nprocs;
+	if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
+	if (rank==0) printf("********************************************************\n");
+
 }
 
-void ScaLBL_MRTModel::VelocityField(double *Vz){
+void ScaLBL_MRTModel::VelocityField(double *VELOCITY){
 
+	Minkowski Morphology(Mask);
 	int SIZE=Np*sizeof(double);
 	ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
 	ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
-	ScaLBL_CopyToHost(&Vz[0],&Velocity[2*Np],SIZE);
+	ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
+
+	memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));
+	Morphology.Initialize();
+	Morphology.UpdateMeshValues();
+	Morphology.ComputeLocal();
+	Morphology.Reduce();
+	
+	double count_loc=0;
+	double count;
+	double vax,vay,vaz;
+	double vax_loc,vay_loc,vaz_loc;
+	vax_loc = vay_loc = vaz_loc = 0.f;
+	for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
+		vax_loc += VELOCITY[n];
+		vay_loc += VELOCITY[Np+n];
+		vaz_loc += VELOCITY[2*Np+n];
+		count_loc+=1.0;
+	}
+	
+	for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
+		vax_loc += VELOCITY[n];
+		vay_loc += VELOCITY[Np+n];
+		vaz_loc += VELOCITY[2*Np+n];
+		count_loc+=1.0;
+	}
+	MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
+	
+	vax /= count;
+	vay /= count;
+	vaz /= count;
+	
+	double mu = (tau-0.5)/3.f;
+	if (rank==0) printf("Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
+	if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Fx, Fy, Fz, mu, 
+						Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
+
 	
 }

models/MRTModel.h

@@ -27,6 +27,7 @@
 #include "common/ScaLBL.h"
 #include "common/Communication.h"
 #include "common/MPI_Helpers.h"
+#include "analysis/Minkowski.h"
 #include "ProfilerApp.h"
 
 class ScaLBL_MRTModel{
@@ -69,6 +70,8 @@ public:
     double *fq;
     double *Velocity;
     double *Pressure;
+    
+    //Minkowski Morphology;
 		
 private:
 	MPI_Comm comm;

tests/CMakeLists.txt

@@ -2,7 +2,7 @@
 #ADD_LBPM_EXECUTABLE( lbpm_nonnewtonian_simulator ) 
 #ADD_LBPM_EXECUTABLE( lbpm_nondarcy_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_color_simulator )
-#ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
+ADD_LBPM_EXECUTABLE( lbpm_permeability_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_BGK_simulator )
 #ADD_LBPM_EXECUTABLE( lbpm_color_macro_simulator )
 ADD_LBPM_EXECUTABLE( lbpm_dfh_simulator )

tests/GenerateSphereTest.cpp

@@ -494,7 +494,7 @@ int main(int argc, char **argv)
 		sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
 		FILE *DIST = fopen(LocalRankFilename,"wb");
 		if (DIST==NULL) ERROR("Error opening file: ID.xxxxx");
-		fwrite(SignDist.data(),1,N,DIST);
+		fwrite(SignDist.data(),1,N*sizeof(double),DIST);
 		fclose(DIST);
 		//......................................................................
 

tests/TestPoiseuille.cpp

@@ -57,7 +57,7 @@ int main(int argc, char **argv)
 		MRT.Create();       // creating the model will create data structure to match the pore structure and allocate variables
 		MRT.Initialize();   // initializing the model will set initial conditions for variables
 		MRT.Run();	 
-		double *Vz;  	Vz= new double [MRT.Np];
+		double *Vz;  	Vz= new double [3*MRT.Np];
 		MRT.VelocityField(Vz);
 
 		if (rank == 0) printf("Force: %f,%f,%f \n",MRT.Fx,MRT.Fy,MRT.Fz);
@@ -80,7 +80,7 @@ int main(int argc, char **argv)
 				n = MRT.Map(i,j,k);
 				//printf("%i,%i,%i; %i :",i,j,k,n);
 				if (n<0) {vz =0.f; printf(" b    "); }
-				else { vz=Vz[n]; printf(" a    "); }
+				else { vz=Vz[n+2*MRT.Np]; printf(" a    "); }
 				printf("%f ",vz);
 				//Analytical solution
 				double x=1.f*i-1.5;
@@ -97,7 +97,7 @@ int main(int argc, char **argv)
 				n = MRT.Map(i,j,k);
 				//printf("%i,%i,%i; %i :",i,j,k,n);
 				if (n<0) {vz =0.f; printf(" b    "); }
-				else { vz=Vz[n]; printf(" a    "); }
+				else { vz=Vz[n+2*MRT.Np]; printf(" a    "); }
 				printf("%f ",vz);
 				//Analytical solution
 				double x=1.f*i-1.5;

tests/lbpm_permeability_simulator.cpp

@@ -10,7 +10,7 @@
 #include "common/Communication.h"
 #include "analysis/TwoPhase.h"
 #include "common/MPI_Helpers.h"
-
+#include "models/MRTModel.h"
 //#define WRITE_SURFACES
 
 /*
@@ -42,396 +42,17 @@ int main(int argc, char **argv)
 			printf("Running Single Phase Permeability Calculation \n");
 			printf("********************************************************\n");
 		}
-
-		// Variables that specify the computational domain  
-		string FILENAME;
-		int Nx,Ny,Nz;		// local sub-domain size
-		int nspheres;		// number of spheres in the packing
-		double Lx,Ly,Lz;	// Domain length
-		double D = 1.0;		// reference length for non-dimensionalization
-		// Color Model parameters
-		int timestepMax, interval;
-		double tau,Fx,Fy,Fz,tol,err;
-		double din,dout;
-		bool pBC,Restart;
-		int i,j,k,n;
-
-		int RESTART_INTERVAL=20000;
-
-		if (rank==0){
-			//.............................................................
-			//		READ SIMULATION PARMAETERS FROM INPUT FILE
-			//.............................................................
-			ifstream input("Permeability.in");
-			// Line 1: model parameters (tau, alpha, beta, das, dbs)
-			input >> tau;			// Viscosity parameter
-			// Line 2: External force components (Fx,Fy, Fz)
-			input >> Fx;
-			input >> Fy;
-			input >> Fz;
-			// Line 3: Pressure Boundary conditions
-			input >> Restart;
-			input >> pBC;
-			input >> din;
-			input >> dout;
-			// Line 4: time-stepping criteria
-			input >> timestepMax;		// max no. of timesteps
-			input >> interval;			// restart interval
-			input >> tol;				// error tolerance
-			//.............................................................
-
-			//.......................................................................
-			// Reading the domain information file
-			//.......................................................................
-			ifstream domain("Domain.in");
-			domain >> nprocx;
-			domain >> nprocy;
-			domain >> nprocz;
-			domain >> Nx;
-			domain >> Ny;
-			domain >> Nz;
-			//domain >> nspheres;
-			domain >> Lx;
-			domain >> Ly;
-			domain >> Lz;
-			//.......................................................................
-
-		}
-		// **************************************************************
-		// Broadcast simulation parameters from rank 0 to all other procs
-		MPI_Barrier(comm);
-		//.................................................
-		MPI_Bcast(&tau,1,MPI_DOUBLE,0,comm);
-		//MPI_Bcast(&pBC,1,MPI_LOGICAL,0,comm);
-		//	MPI_Bcast(&Restart,1,MPI_LOGICAL,0,comm);
-		MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
-		MPI_Bcast(&dout,1,MPI_DOUBLE,0,comm);
-		MPI_Bcast(&Fx,1,MPI_DOUBLE,0,comm);
-		MPI_Bcast(&Fy,1,MPI_DOUBLE,0,comm);
-		MPI_Bcast(&Fz,1,MPI_DOUBLE,0,comm);
-		MPI_Bcast(&timestepMax,1,MPI_INT,0,comm);
-		MPI_Bcast(&interval,1,MPI_INT,0,comm);
-		MPI_Bcast(&tol,1,MPI_DOUBLE,0,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_Barrier(comm);
-
-		RESTART_INTERVAL=interval;
-		// **************************************************************
-		// **************************************************************
-		double rlx_setA = 1.f/tau;
-		double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
-
-		if (nprocs != nprocx*nprocy*nprocz){
-			printf("nprocx =  %i \n",nprocx);
-			printf("nprocy =  %i \n",nprocy);
-			printf("nprocz =  %i \n",nprocz);
-			INSIST(nprocs == nprocx*nprocy*nprocz,"Fatal error in processor count!");
-		}
-
-		if (rank==0){
-			printf("********************************************************\n");
-			printf("tau = %f \n", tau);
-			printf("Force(x) = %.5g \n", Fx);
-			printf("Force(y) = %.5g \n", Fy);
-			printf("Force(z) = %.5g \n", Fz);
-			printf("Sub-domain size = %i x %i x %i\n",Nx,Ny,Nz);
-			printf("Process grid = %i x %i x %i\n",nprocx,nprocy,nprocz);
-			printf("********************************************************\n");
-		}
-
-		double viscosity=(tau-0.5)/3.0;
-		// Initialized domain and averaging framework for Two-Phase Flow
-		int BC=pBC;
-		Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
-		for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
-		Dm.CommInit();
-		TwoPhase Averages(Dm);
-
-		// Mask that excludes the solid phase
-		Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
-		MPI_Barrier(comm);
-
-		Nx += 2;	Ny += 2;	Nz += 2;
-		int N = Nx*Ny*Nz;
-
-		//.......................................................................
-		if (rank == 0)	printf("Read input media... \n");
-		//.......................................................................
-
-		//.......................................................................
-		// Filenames used
-		char LocalRankString[8];
-		char LocalRankFilename[40];
-		char LocalRestartFile[40];
-		char tmpstr[10];
-		sprintf(LocalRankString,"%05d",rank);
-		sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
-		sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
-
-		//	printf("Local File Name =  %s \n",LocalRankFilename);
-		// .......... READ THE INPUT FILE .......................................
-		//	char value;
-		char *id;
-		id = new char[N];
-		double sum, 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));
-		double porosity, pore_vol;
-		//...........................................................................
-		if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
-		//.......................................................................
-		// Read the signed distance
-		sprintf(LocalRankString,"%05d",rank);
-		sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
-		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
-		//.......................................................................
-
-		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.f;
-		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 (rank==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
-		sprintf(LocalRankFilename,"ID.%05i",rank);
-		size_t readID;
-		FILE *IDFILE = fopen(LocalRankFilename,"rb");
-		if (IDFILE==NULL) ERROR("lbpm_permeability_simulator: Error opening file: ID.xxxxx");
-		readID=fread(id,1,N,IDFILE);
-		if (readID != size_t(N)) printf("lbpm_permeability_simulator: Error reading ID (rank=%i) \n",rank);
-		fclose(IDFILE);
 		
-		//.......................................................................
-		// Compute the media porosity, assign phase labels and solid composition
-		//.......................................................................
-		sum_local=0.0;
-		int Np=0;  // number of local pore nodes
-		//.......................................................................
-		for (k=1;k<Nz-1;k++){
-			for (j=1;j<Ny-1;j++){
-				for (i=1;i<Nx-1;i++){
-					n = k*Nx*Ny+j*Nx+i;
-					if (id[n] > 0){
-						sum_local+=1.0;
-						Np++;
-					}
-				}
-			}
-		}
-		MPI_Allreduce(&sum_local,&sum,1,MPI_DOUBLE,MPI_SUM,comm);
-		porosity = sum*iVol_global;
-		if (rank==0) printf("Media porosity = %f \n",porosity);
-
-		//.........................................................
-		// don't perform computations at the eight corners
-		id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
-		id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
-		//.........................................................
-		MPI_Barrier(comm);
-
-		// Initialize communication structures in averaging domain
-		for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
-		Mask.CommInit(comm);
-
-		//...........................................................................
-		if (rank==0)	printf ("Create ScaLBL_Communicator \n");
-		// Create a communicator for the device
-
-		ScaLBL_Communicator ScaLBL_Comm(Mask);
-
-		// LBM variables
-		if (rank==0)	printf ("Allocating distributions \n");
-		int Npad=(Np/16 + 2)*16;
-
-		int *neighborList;
-		IntArray Map(Nx,Ny,Nz);
-		neighborList= new int[18*Npad];
-		Np = ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Mask.id,Np);
-		MPI_Barrier(comm);
-
-		//......................device distributions.................................
-		int dist_mem_size = Np*sizeof(double);
-		int neighborSize=18*(Np*sizeof(int));
-
-		int *NeighborList;
-		//		double *f_even,*f_odd;
-		double * dist;
-		double * Velocity;
-		double * Pressure;
-		//...........................................................................
-		ScaLBL_AllocateDeviceMemory((void **) &dist, 19*dist_mem_size);
-		ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
-		ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np);
-		ScaLBL_AllocateDeviceMemory((void **) &Pressure, 3*sizeof(double)*Np);
-		ScaLBL_CopyToDevice(NeighborList,     neighborList, neighborSize);
-		//...........................................................................
-
-		//...........................................................................
-		if (rank==0)	printf("Setting the distributions, size = %i\n", N);
-		//...........................................................................
-
-		// Finalize setup for averaging domain
-		//Averages.SetupCubes(Dm);
-		Averages.UpdateSolid();
-		// Initialize two phase flow variables (all wetting phase)
-		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;
-					Averages.Phase(i,j,k) = -1.0;
-					Averages.SDn(i,j,k) = Averages.Phase(i,j,k);
-					Averages.Phase_tplus(i,j,k) = Averages.SDn(i,j,k);
-					Averages.Phase_tminus(i,j,k) = Averages.SDn(i,j,k);
-					Averages.DelPhi(i,j,k) = 0.0;
-					Averages.Press(i,j,k) = 0.0;
-					Averages.Vel_x(i,j,k) = 0.0;
-					Averages.Vel_y(i,j,k) = 0.0;
-					Averages.Vel_z(i,j,k) = 0.0;
-				}
-			}
-		}
-
-		//.......................................................................
-
-		ScaLBL_D3Q19_Init(dist, Np);
-
-		int timestep = 0;
-		if (rank==0) printf("********************************************************\n");
-		if (rank==0)	printf("No. of timesteps: %i \n", timestepMax);
-
-		//.......create and start timer............
-		double starttime,stoptime,cputime;
-		MPI_Barrier(comm);
-		starttime = MPI_Wtime();
-		//.........................................
-
-		double D32,Fo,Re,velocity,err1D,mag_force,vel_prev;
-		err = vel_prev = 1.0;
-		if (rank==0) printf("Begin timesteps: error tolerance is %f \n", tol);
-		//************ MAIN ITERATION LOOP ***************************************/
-		while (timestep < timestepMax && err > tol ){
-
-			timestep++;
-			ScaLBL_Comm.SendD3Q19AA(dist); //READ FROM NORMAL
-			ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
-			ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
-			ScaLBL_D3Q19_AAodd_MRT(NeighborList, dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
-			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
-
-			timestep++;
-			ScaLBL_Comm.SendD3Q19AA(dist); //READ FORM NORMAL
-			ScaLBL_D3Q19_AAeven_MRT(dist, ScaLBL_Comm.first_interior, ScaLBL_Comm.last_interior, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
-			ScaLBL_Comm.RecvD3Q19AA(dist); //WRITE INTO OPPOSITE
-			ScaLBL_D3Q19_AAeven_MRT(dist, 0, ScaLBL_Comm.next, Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
-			ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
-			//************************************************************************/
-
-			if (timestep%500 == 0){
-				//...........................................................................
-				// Copy the data for for the analysis timestep
-				//...........................................................................
-				// Copy the phase from the GPU -> CPU
-				//...........................................................................
-				ScaLBL_DeviceBarrier();
-		        ScaLBL_D3Q19_Pressure(dist,Pressure,Np);
-		        ScaLBL_D3Q19_Momentum(dist,Velocity,Np);
-		        
-		        ScaLBL_Comm.RegularLayout(Map,Pressure,Averages.Press);
-		        ScaLBL_Comm.RegularLayout(Map,&Velocity[0],Averages.Vel_x);
-		        ScaLBL_Comm.RegularLayout(Map,&Velocity[Np],Averages.Vel_y);
-		        ScaLBL_Comm.RegularLayout(Map,&Velocity[2*Np],Averages.Vel_z);
-
-				// Way more work than necessary -- this is just to get the solid interfacial area!!
-				Averages.Initialize();
-				Averages.UpdateMeshValues();
-				Averages.ComputeLocal();
-				Averages.Reduce();
-
-				double vawx = Averages.vaw_global(0);
-				double vawy = Averages.vaw_global(1);
-				double vawz = Averages.vaw_global(2);
-				if (rank==0){
-					// ************* DIMENSIONLESS FORCHEIMER EQUATION *************************
-					//  Dye, A.L., McClure, J.E., Gray, W.G. and C.T. Miller
-					//  Description of Non-Darcy Flows in Porous Medium Systems
-					//  Physical Review E 87 (3), 033012
-					//  Fo := density*D32^3*(density*force) / (viscosity^2)
-					//	Re := density*D32*velocity / viscosity
-					//  Fo = a*Re + b*Re^2
-					// *************************************************************************
-					//viscosity = (tau-0.5)*0.333333333333333333;
-					D32 = 6.0*(Dm.Volume-Averages.vol_w_global)/Averages.As_global;
-					printf("Sauter Mean Diameter = %f \n",D32);
-					mag_force = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
-					Fo = D32*D32*D32*mag_force/viscosity/viscosity;
-					// .... 1-D flow should be aligned with force ...
-					velocity = vawx*Fx/mag_force + vawy*Fy/mag_force + vawz*Fz/mag_force;
-					err1D = fabs(velocity-sqrt(vawx*vawx+vawy*vawy+vawz*vawz))/velocity;
-					//.......... Computation of the Reynolds number Re ..............
-					Re = D32*velocity/viscosity;
-					printf("Force: %.5g,%.5g,%.5g \n",Fx,Fy,Fz);
-					printf("Velocity: %.5g,%.5g,%.5g \n",vawx,vawy,vawz);
-					printf("Relative error for 1D representation: %.5g \n",err1D);
-					printf("Dimensionless force: %5g \n", Fo);
-					printf("Reynolds number: %.5g \n", Re);
-					printf("Dimensionless Permeability (k/D^2): %.5g \n", Re/Fo);
-				}
-
-			}
-		}
-		//************************************************************************/
-		ScaLBL_DeviceBarrier();
-		MPI_Barrier(comm);
-		stoptime = MPI_Wtime();
-		if (rank==0) printf("-------------------------------------------------------------------\n");
-		// Compute the walltime per timestep
-		cputime = (stoptime - starttime)/timestep;
-		// Performance obtained from each node
-		double MLUPS = double(Np)/cputime/1000000;
-
-		if (rank==0) printf("********************************************************\n");
-		if (rank==0) printf("CPU time = %f \n", cputime);
-		if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS);
-		MLUPS *= nprocs;
-		if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS);
-		if (rank==0) printf("********************************************************\n");
-
-		NULL_USE(RESTART_INTERVAL);
+		ScaLBL_MRTModel MRT(rank,nprocs,comm);
+		auto filename = argv[1];
+		MRT.ReadParams(filename);
+		MRT.SetDomain();    // this reads in the domain 
+		MRT.ReadInput();
+		MRT.Create();       // creating the model will create data structure to match the pore structure and allocate variables
+		MRT.Initialize();   // initializing the model will set initial conditions for variables
+		MRT.Run();	 
+		double *Velocity;  	Velocity= new double [3*MRT.Np];
+		MRT.VelocityField(Velocity);
 	}
 	// ****************************************************
 	MPI_Barrier(comm);

tests/lbpm_refine_pp.cpp

@@ -31,37 +31,28 @@ int main(int argc, char **argv)
 		int i,j,k,n;
 		int BC=0;
 
-		if ( rank==0 ) {
-			printf("lbpm_refine_pp: Refining signed distance mesh (x2) \n");
+		string filename;
+		if (argc > 1){
+			filename=argv[1];
 		}
-		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;
+		else ERROR("No input database provided\n");
+		// read the input database 
+		auto db = std::make_shared<Database>( filename );
+		auto domain_db = db->getDatabase( "Domain" );
 
-		}
-		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_Barrier(comm);
+		// Read domain parameters
+		auto L = domain_db->getVector<double>( "L" );
+		auto size = domain_db->getVector<int>( "n" );
+		auto nproc = domain_db->getVector<int>( "nproc" );
+		auto ReadValues = domain_db->getVector<char>( "ReadValues" );
+		auto WriteValues = domain_db->getVector<char>( "WriteValues" );
+		
+		nx = size[0];
+		ny = size[1];
+		nz = size[2];
+		nprocx = nproc[0];
+		nprocy = nproc[1];
+		nprocz = nproc[2];
 
 		// Check that the number of processors >= the number of ranks
 		if ( rank==0 ) {
@@ -112,7 +103,15 @@ int main(int argc, char **argv)
 		ReadSignDist=fread(SignDist.data(),8,N,DIST);
 		if (ReadSignDist != size_t(N)) printf("lbpm_refine_pp: Error reading signed distance function (rank=%i)\n",rank);
 		fclose(DIST);
-
+		
+		char *Labels;
+		Labels = new char[N];
+		sprintf(LocalRankFilename,"ID.%05i",rank);
+		FILE *LABELS = fopen(LocalRankFilename,"rb");
+		size_t ReadLabels;
+		ReadLabels=fread(Labels,1,N,LABELS);
+		if (ReadLabels != size_t(N)) printf("lbpm_refine_pp: Error reading ID  (rank=%i)\n",rank);
+		fclose(LABELS);
 
 		if ( rank==0 )   printf("Set up Domain, read input distance \n");
 
@@ -130,6 +129,9 @@ int main(int argc, char **argv)
 		}
 
 		int ri,rj,rk,rn; //refined mesh indices
+		//char *RefineLabel;
+		//RefineLabel = new char [rnx*rny*rnz];
+		Array <char> RefineLabel(rnx,rny,rnz);
 		for (rk=1; rk<rnz-1; rk++){
 			for (rj=1; rj<rny-1; rj++){
 				for (ri=1; ri<rnx-1; ri++){
@@ -146,6 +148,7 @@ int main(int argc, char **argv)
 					pt.y=0.5*(rj-1)+1.f;
 					pt.z=0.5*(rk-1)+1.f;
 					RefinedSignDist(ri,rj,rk) = LocalApprox.eval(pt);
+					RefineLabel(ri,rj,rk) = Labels[k*nx*ny+j*nx+i]; 
 				}
 			}
 		}
@@ -175,6 +178,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i,j,k);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i,j,k);
 				}
 			}
 		}
@@ -183,16 +188,17 @@ int main(int argc, char **argv)
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
 
-		for (int k=0; k<nz; k++){
+/*		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					if (BlockDist(i,j,k) > 0.f)
 						id[k*nx*ny + j*nx + i]=2;
 					else
-						id[k*nx*ny + j*nx + i]=0;
+						id[k*nx*ny + j*nx + i]= 0;
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);
@@ -203,6 +209,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j,k);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j,k);
 				}
 			}
 		}
@@ -211,7 +219,7 @@ int main(int argc, char **argv)
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
 
-		for (int k=0; k<nz; k++){
+/*		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					if (BlockDist(i,j,k) > 0.f)
@@ -221,6 +229,7 @@ int main(int argc, char **argv)
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);
@@ -232,6 +241,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j+ny-2,k);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j+ny-2,k);
 				}
 			}
 		}
@@ -240,7 +251,7 @@ int main(int argc, char **argv)
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
 
-		for (int k=0; k<nz; k++){
+/*		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					if (BlockDist(i,j,k) > 0.f)
@@ -250,6 +261,7 @@ int main(int argc, char **argv)
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);
@@ -260,6 +272,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i,j+ny-2,k);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i,j+ny-2,k);
 				}
 			}
 		}
@@ -267,7 +281,7 @@ int main(int argc, char **argv)
 		REFINEDIST = fopen(LocalRankFilename,"wb");
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
-
+/*
 		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
@@ -278,6 +292,7 @@ int main(int argc, char **argv)
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);
@@ -288,6 +303,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i,j,k+nz-2);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i,j,k+nz-2);
 				}
 			}
 		}
@@ -295,7 +312,7 @@ int main(int argc, char **argv)
 		REFINEDIST = fopen(LocalRankFilename,"wb");
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
-
+/*
 		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
@@ -306,6 +323,7 @@ int main(int argc, char **argv)
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);
@@ -316,6 +334,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j,k+nz-2);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j,k+nz-2);
 				}
 			}
 		}
@@ -324,7 +344,7 @@ int main(int argc, char **argv)
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
 
-		for (int k=0; k<nz; k++){
+/*		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					if (BlockDist(i,j,k) > 0.f)
@@ -334,6 +354,7 @@ int main(int argc, char **argv)
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);
@@ -344,6 +365,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i,j+ny-2,k+nz-2);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i,j+ny-2,k+nz-2);
 				}
 			}
 		}
@@ -351,7 +374,7 @@ int main(int argc, char **argv)
 		REFINEDIST = fopen(LocalRankFilename,"wb");
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
-
+/*
 		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
@@ -362,6 +385,7 @@ int main(int argc, char **argv)
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);
@@ -372,6 +396,8 @@ int main(int argc, char **argv)
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j+ny-2,k+nz-2);
+					if (BlockDist(i,j,k) > 0) 	id[k*nx*ny + j*nx + i]=2;
+					else 						id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j+ny-2,k+nz-2);
 				}
 			}
 		}
@@ -380,7 +406,8 @@ int main(int argc, char **argv)
 		REFINEDIST = fopen(LocalRankFilename,"wb");
 		fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
 		fclose(REFINEDIST);
-		for (int k=0; k<nz; k++){
+		
+/*		for (int k=0; k<nz; k++){
 			for (int j=0; j<ny; j++){
 				for (int i=0; i<nx; i++){
 					if (BlockDist(i,j,k) > 0.f)
@@ -390,6 +417,7 @@ int main(int argc, char **argv)
 				}
 			}
 		}
+		*/
 		sprintf(LocalRankFilename,"RefineID.%05i",writerank);
 		WRITEID = fopen(LocalRankFilename,"wb");
 		fwrite(id,1,nx*ny*nz,WRITEID);