merging with morphLBM (may be some problems still)

analysis/morphology.cpp

@@ -131,7 +131,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 
 	// Increase the critical radius until the target saturation is met
 	double deltaR=0.05; // amount to change the radius in voxel units
-	double Rcrit_old;
+	double Rcrit_old=0.0;
 
 	int imin,jmin,kmin,imax,jmax,kmax;
 

common/Domain.cpp

@@ -578,6 +578,51 @@ void Domain::Decomp( const std::string& Filename )
 		Comm.recv(id,N,0,15);
 	}
 	Comm.barrier();
+
+	// Compute the porosity
+	double sum;
+	double sum_local=0.0;
+	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));
+	//.........................................................
+	// If external boundary conditions are applied remove solid
+	if (BoundaryCondition >  0  && kproc() == 0){
+    	if (inlet_layers_z < 4)	inlet_layers_z=4;
+		for (int k=0; k<inlet_layers_z; k++){
+			for (int j=0;j<Ny;j++){
+				for (int i=0;i<Nx;i++){
+					int n = k*Nx*Ny+j*Nx+i;
+					id[n] = 1;
+				}                    
+			}
+ 		}
+ 	}
+    if (BoundaryCondition >  0  && kproc() == nprocz-1){
+    	if (outlet_layers_z < 4)	outlet_layers_z=4;
+ 		for (int k=Nz-outlet_layers_z; k<Nz; k++){
+ 			for (int j=0;j<Ny;j++){
+ 				for (int i=0;i<Nx;i++){
+ 					int n = k*Nx*Ny+j*Nx+i;
+ 					id[n] = 2;
+ 				}                    
+ 			}
+ 		}
+ 	}
+    for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
+        for (int j=1;j<Ny-1;j++){
+            for (int i=1;i<Nx-1;i++){
+                int n = k*Nx*Ny+j*Nx+i;
+                if (id[n] > 0){
+                    sum_local+=1.0;
+                }
+            }
+        }
+    }
+    sum = Comm.sumReduce(sum_local);
+    porosity = sum*iVol_global;
+    if (rank()==0) printf("Media porosity = %f \n",porosity);
+ 	//.........................................................
+
 }
 
 void Domain::AggregateLabels( const std::string& filename ){

cpu/D3Q19.cpp

@@ -687,7 +687,7 @@ extern "C" void ScaLBL_D3Q19_Pressure_BC_z(int *list, double *dist, double din,
 	double f0,f1,f2,f3,f4,f5,f6,f7,f8,f9;
 	double f10,f11,f12,f13,f14,f15,f16,f17,f18;
 	double ux,uy,uz;
-	double Cxz,Cyz;
+	double Cxz,Cyz;	
 
 	for (int idx=0; idx<count; idx++){
 		n = list[idx];

example/Workflow/HelperFunctions.R

@@ -7,11 +7,26 @@ ReadSubphase<-function(PATH){
   S<-read.csv(FILE,head=TRUE,sep=" ")
   S$Vw<-S$Vwc+S$Vwd
   S$Vn<-S$Vnc+S$Vnd
+  S$Aw<-S$Awc+S$Awd
+  S$An<-S$Anc+S$And
+  S$Hw<-S$Hwc+S$Hwd
+  S$Hn<-S$Hnc+S$Hnd
+  S$Xw<-S$Xwc+S$Xwd
+  S$Xn<-S$Xnc+S$Xnd
+
   S$Sw<-S$Vw/(S$Vn+S$Vw)
+  S$pw<-(S$pwc*S$Vwc+S$pwd*S$Vwd) / (S$Vwc+S$Vwd)
+  S$pn<-(S$pnc*S$Vnc+S$pnd*S$Vnd) / (S$Vnc+S$Vnd)
+
   S$Qwx<-S$Vw*(S$Pwc_x+S$Pwd_x)/(S$Mwc+S$Mwd)
   S$Qnx<-S$Vn*(S$Pnc_x+S$Pnd_x)/(S$Mnc+S$Mnd)
-  S$Krn<-S$nun*S$Qnx/S$Fx
-  S$Krw<-S$nuw*S$Qwx/S$Fx
+  S$Qwy<-S$Vw*(S$Pwc_y+S$Pwd_y)/(S$Mwc+S$Mwd)
+  S$Qny<-S$Vn*(S$Pnc_y+S$Pnd_y)/(S$Mnc+S$Mnd)
+  S$Qwz<-S$Vw*(S$Pwc_z+S$Pwd_z)/(S$Mwc+S$Mwd)
+  S$Qnz<-S$Vn*(S$Pnc_z+S$Pnd_z)/(S$Mnc+S$Mnd)
+
+  S$Krn<-S$nun*S$Qnz/S$Fz
+  S$Krw<-S$nuw*S$Qwz/S$Fz
   S$Case<-PATH
   return(S)
 }

models/ColorModel.cpp

@@ -938,7 +938,7 @@ void ScaLBL_ColorModel::Run(){
 					delta_volume = volA*Dm->Volume - initial_volume;
 					CURRENT_MORPH_TIMESTEPS += analysis_interval;
 					double massChange = SeedPhaseField(seed_water);
-					if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target);
+					if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", massChange, delta_volume, delta_volume_target);
 				}
 				else if (USE_MORPHOPEN_OIL){
 					delta_volume = volA*Dm->Volume - initial_volume;

models/MRTModel.cpp

@@ -93,15 +93,19 @@ void ScaLBL_MRTModel::SetDomain(){
 }
 
 void ScaLBL_MRTModel::ReadInput(){
-    int rank=Dm->rank();
-    //.......................................................................
-    //.......................................................................
-    Mask->ReadIDs();
     
     sprintf(LocalRankString,"%05d",Dm->rank());
     sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
     sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
 
+	if (domain_db->keyExists( "Filename" )){
+		auto Filename = domain_db->getScalar<std::string>( "Filename" );
+		Mask->Decomp(Filename);
+	}
+	else{
+		Mask->ReadIDs();
+	}
+
 	// Generate the signed distance map
 	// Initialize the domain and communication
 	Array<char> id_solid(Nx,Ny,Nz);

tests/GenerateSphereTest.cpp

@@ -145,8 +145,8 @@ inline void MorphOpen(DoubleArray SignDist, char *id, Domain &Dm, int nx, int ny
 
 	// Increase the critical radius until the target saturation is met
 	double deltaR=0.05; // amount to change the radius in voxel units
-	double Rcrit_old;
-	double Rcrit_new;
+	double Rcrit_old=0.0;
+	double Rcrit_new=0.0;
 
 	int imin,jmin,kmin,imax,jmax,kmax;