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

analysis/SubPhase.cpp

@@ -245,8 +245,8 @@ void SubPhase::Basic(){
 			force_mag = 1.0;
 		}
 		double saturation=gwb.V/(gwb.V + gnb.V);
-		double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume;
+		double water_flow_rate=gwb.V*sqrt(gwb.Px*gwb.Px + gwb.Py*gwb.Py + gwb.Pz*gwb.Pz)/gwb.M / Dm->Volume;
-		double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume;
+		double not_water_flow_rate=gnb.V*sqrt(gnb.Px*gnb.Px + gnb.Py*gnb.Py + gnb.Pz*gnb.Pz)/gnb.M/ Dm->Volume;
 		double total_flow_rate = water_flow_rate + not_water_flow_rate;
 		double fractional_flow= water_flow_rate / total_flow_rate;
 

analysis/morphology.cpp

@@ -279,10 +279,10 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain>
 		MPI_Allreduce(&count,&countGlobal,1,MPI_DOUBLE,MPI_SUM,Dm->Comm);
 		void_fraction_new = countGlobal/totalGlobal;
 		void_fraction_diff_new = abs(void_fraction_new-VoidFraction);
-		if (rank==0){
+	/*	if (rank==0){
 			printf("     %f ",void_fraction_new);
 			printf("     %f\n",Rcrit_new);
-		}
+		} */
 	}
 
 	if (void_fraction_diff_new<void_fraction_diff_old){

common/Domain.cpp

@@ -119,11 +119,6 @@ void Domain::initialize( std::shared_ptr<Database> db )
     auto nproc = d_db->getVector<int>("nproc");
     auto n = d_db->getVector<int>("n");
     
-    voxel_length = 1.0;
-    if (d_db->keyExists( "voxel_length" )){
-    	auto voxel_length = d_db->getScalar<double>("voxel_length");
-    }
-    
 	if (d_db->keyExists( "InletLayers" )){
 		auto InletCount = d_db->getVector<int>( "InletLayers" );
 		inlet_layers_x = InletCount[0];
@@ -136,6 +131,10 @@ void Domain::initialize( std::shared_ptr<Database> db )
 		outlet_layers_y = OutletCount[1];
 		outlet_layers_z = OutletCount[2];
 	}
+    voxel_length = 1.0;
+    if (d_db->keyExists( "voxel_length" )){
+    	voxel_length = d_db->getScalar<double>("voxel_length");
+    }
 
     ASSERT( n.size() == 3u );
     ASSERT( nproc.size() == 3u );
@@ -163,6 +162,9 @@ void Domain::initialize( std::shared_ptr<Database> db )
     // Fill remaining variables
 	N = Nx*Ny*Nz;
 	Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
+
+	if (myrank==0) printf("voxel length = %f micron \n", voxel_length);
+
 	id = new signed char[N];
 	memset(id,0,N);
 	BoundaryCondition = d_db->getScalar<int>("BC");

example/Bubble/input.db

@@ -7,15 +7,7 @@ Color {
     beta  = 0.95;
     F = 0, 0, 0
     Restart = false
-    pBC = 0
-    din = 1.0
-    dout = 1.0
-    flux = 1.0e-3
     timestepMax = 200
-    interval = 1000
-    tol = 1e-5;
-    das = 0.1
-    dbs = 0.9
     ComponentLabels = 0, 1, 2
     ComponentAffinity = -1.0, 1.0, -1.0
 }

example/Piston/Color.in

@@ -1,6 +0,0 @@
-1.0 1.0
-1.0 1.0
-1.0e-3 0.95
-0.0 0.0 0.0
-0 4 10.0 1.0
-10000 2000 1e-5

example/Piston/ComponentLabels.csv

@@ -1,3 +0,0 @@
-0 -0.5
-1 1.0
-2 -1.0

example/Piston/Domain.in

@@ -1,4 +0,0 @@
-1 1 10
-40 40 40
-229
-1.0 1.0 1.0

example/Piston/input.db

@@ -7,14 +7,7 @@ Color {
     beta  = 0.95;
     F = 0, 0, 0
     Restart = false
-    pBC = 0
-    din = 1.0
-    dout = 1.0
     timestepMax = 3000
-    interval = 1000
-    tol = 1e-5;
-    das = 0.1
-    dbs = 0.9
     flux = 2.0
     ComponentLabels = 0
     ComponentAffinity = -1.0;
@@ -25,7 +18,6 @@ Domain {
     nproc = 1, 1, 4     // Number of processors (Npx,Npy,Npz)
     n = 24, 24, 24      // Size of local domain (Nx,Ny,Nz)
     N = 24, 24, 96      // size of the input image
-    n_spheres = 1       // Number of spheres
     L = 1, 1, 1         // Length of domain (x,y,z)
     BC = 4              // Boundary condition type
     ReadType = "8bit"   

example/Plates/input.db

@@ -7,14 +7,7 @@ Color {
     beta  = 0.95;
     F = 0, 0, 0
     Restart = false
-    pBC = 0
-    din = 1.0
-    dout = 1.0
     timestepMax = 3000
-    interval = 1000
-    tol = 1e-5;
-    das = 0.1
-    dbs = 0.9
     flux = 0.0
     ComponentLabels = -2, -1
     ComponentAffinity = -1.0, -0.5;
@@ -25,7 +18,6 @@ Domain {
     nproc = 1, 1, 4     // Number of processors (Npx,Npy,Npz)
     n = 24, 24, 24      // Size of local domain (Nx,Ny,Nz)
     N = 24, 24, 96      // size of the input image
-    n_spheres = 1       // Number of spheres
     L = 1, 1, 1         // Length of domain (x,y,z)
     BC = 0              // Boundary condition type
     ReadType = "8bit"

example/Poiseuille/Domain.in

@@ -1,3 +0,0 @@
-1 1 1
-100 100 100
-1.0 1.0 1.0

example/Poiseuille/Permeability.in

@@ -1,4 +0,0 @@
-1.0
-0.0 0.0 1.0e-5
-0 0 1.0 1.0
-1000 1000 1.0e-5

example/Sph125/Domain.in

@@ -1,3 +0,0 @@
-2 2 2
-100 100 100
-1.0 1.0 1.0

example/Workflow/HelperFunctions.R

@@ -1,7 +1,6 @@
 require("ggplot2")
 
-IngestRelperm<-function(PATH){
+ReadSubphase<-function(PATH){
-
   FILE=paste0(PATH,"/subphase.csv")
   S<-read.csv(FILE,head=TRUE,sep=" ")
   S$Vw<-S$Vwc+S$Vwd
@@ -15,14 +14,28 @@ IngestRelperm<-function(PATH){
   return(S)
 }
 
+ReadTimelog<-function(PATH){
+	FILE=paste0(PATH,"/timelog.csv")
+	D<-read.csv(file=FILE,head=TRUE,sep=" ")
+	D$time<-seq(0,nrow(D))
+	return(D)
+}
+
 ReadRelperm<-function(PATH){
 	FILE=paste0(PATH,"/relperm.csv")
-	D<-read.csv(file=FILE,head=FALSE,sep=" ")
+	D<-read.csv(file=FILE,head=TRUE,sep=" ")
-	colnames(D)<-c("time","nun","nuw","ift","Fx","Fy","Fz","Vn","Vw","unx","uny","unz","uwx","uwy","uwz")
+	D$Case<-PATH
-	D$Sw<-D$Vw/(D$Vn+D$Vw)
+
-	D$Krn<-D$Vn*D$nun*D$unx/D$Fx
+	p<-ggplot(D)+
-	D$Krw<-D$Vw*D$nuw*D$uwx/D$Fx
+	  geom_line(aes(sat.water,eff.perm.oil,color="oil"))+
-	subset(D,D$time>100000)
+	  geom_line(aes(sat.water,eff.perm.water,color="water"))+
+	  xlab("Water saturation")+ylab("Effective permeability")+
+	  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+
+	  theme_bw()
+
+	FILE=paste0(PATH,"-relperm.png")
+ 	ggsave(FILE,p,height=4.0,width=6.0)
+
 	return(D)
 }
 
@@ -39,18 +52,18 @@ ReadCase<-function(PATTERN){
 }
 
 
-D<-ReadCase("benth_w")
+ReadMorphDrain<-function(PATH){
-require("ggplot2")
+	FILE=paste0(PATH,"/morphdrain.csv")
-
+	B<-read.csv(FILE,head=TRUE,sep=" ")
-B<-read.csv("bentheimer/drain.csv",head=TRUE,sep=" ")
+	B$Media<-PATH
-B$Media<-"Bentheimer"
-M1<-read.csv("mineral_model_1/drain.csv",head=TRUE,sep=" ")
-M1$Media<-"Mineral Sandstone #1"
 
 	p<-ggplot()+
 	  geom_line(data=B,aes(Sw,2/R,colour=Media))+
-	geom_line(data=M1,aes(Sw,2/R,colour=Media))+
 	  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+
 	  theme_bw()
 
-ggsave("morph-drain.png",p,height=4.0,width=6.0)
+	FILE=paste0(PATH,"-morphdrain.png")
+ 	ggsave(FILE,p,height=4.0,width=6.0)
+	return(B)
+}
+

models/ColorModel.cpp

@@ -451,6 +451,7 @@ void ScaLBL_ColorModel::Run(){
 	bool MORPH_ADAPT = false;
 	bool USE_MORPH = false;
 	bool USE_SEED = false;
+	bool USE_MORPHOPEN_OIL = false;
 	int analysis_interval = 1000; 	// number of timesteps in between in situ analysis 
 	int MAX_MORPH_TIMESTEPS = 50000; // maximum number of LBM timesteps to spend in morphological adaptation routine
 	int MIN_STEADY_TIMESTEPS = 100000;
@@ -483,15 +484,22 @@ void ScaLBL_ColorModel::Run(){
 	}
 	if (analysis_db->keyExists( "seed_water" )){
 		seed_water = analysis_db->getScalar<double>( "seed_water" );
+		if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water);
 		USE_SEED = true;
 	}
 	if (analysis_db->keyExists( "morph_delta" )){
 		morph_delta = analysis_db->getScalar<double>( "morph_delta" );
+		if (rank == 0) printf("Target volume change %f (morph_delta) \n",morph_delta);
 	}
 	if (analysis_db->keyExists( "morph_interval" )){
 		morph_interval = analysis_db->getScalar<int>( "morph_interval" );
 		USE_MORPH = true;
 	}
+	if (analysis_db->keyExists( "use_morphopen_oil" )){
+		USE_MORPHOPEN_OIL = analysis_db->getScalar<bool>( "use_morphopen_oil" );
+		if (rank == 0 && USE_MORPHOPEN_OIL) printf("Volume change by morphological opening \n");
+		USE_MORPH = true;
+	}
 	if (analysis_db->keyExists( "tolerance" )){
 		tolerance = analysis_db->getScalar<double>( "tolerance" );
 	}
@@ -641,8 +649,8 @@ void ScaLBL_ColorModel::Run(){
 				force_mag = 1.0;
 			}
 			double current_saturation = volB/(volA+volB);
-			double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
+			double flow_rate_A = volA*sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
-			double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
+			double flow_rate_B = volB*sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
 			double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
 			
 			if ( morph_timesteps > morph_interval ){
@@ -656,7 +664,7 @@ void ScaLBL_ColorModel::Run(){
 				if ( isSteady ){
 					MORPH_ADAPT = true;
 					CURRENT_MORPH_TIMESTEPS=0;
-					delta_volume_target = (volA )*morph_delta; // set target volume change
+					delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
 					Averages->Full();
 					Averages->Write(timestep);
 					analysis.WriteVisData( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
@@ -690,7 +698,6 @@ void ScaLBL_ColorModel::Run(){
 
 						printf("  Measured capillary number %f \n ",Ca);
 					}
-
 					if (SET_CAPILLARY_NUMBER ){
 						Fx *= capillary_number / Ca;
 						Fy *= capillary_number / Ca;
@@ -709,6 +716,7 @@ void ScaLBL_ColorModel::Run(){
 						if (rank == 0) printf("    -- adjust force by factor %f \n ",capillary_number / Ca);
 						Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);
 					}
+					CURRENT_STEADY_TIMESTEPS = 0;
 				}
 				else{
 					if (rank==0){
@@ -723,10 +731,15 @@ void ScaLBL_ColorModel::Run(){
 			if (MORPH_ADAPT ){
 				CURRENT_MORPH_TIMESTEPS += analysis_interval;
 				if (USE_SEED){
-					delta_volume = volA - initial_volume;
+					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, mass change %f ***\n", seed_water, massChange);
+					if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target);
+				}
+				else if (USE_MORPHOPEN_OIL){
+					delta_volume = volA*Dm->Volume - initial_volume;
+					if (rank==0) printf("***Morphological opening of connected oil, with target volume change %f ***\n", delta_volume_target);
+					MorphOpenConnected(delta_volume_target);
 				}
 				else{
 					if (rank==0) printf("***Morphological step with target volume change %f ***\n", delta_volume_target);
@@ -735,13 +748,13 @@ void ScaLBL_ColorModel::Run(){
 				}
 				if ( (delta_volume - delta_volume_target)/delta_volume_target > 0.0 ){
 					MORPH_ADAPT = false;
-					initial_volume = volA;
+					initial_volume = volA*Dm->Volume;
 					delta_volume = 0.0;
 					CURRENT_STEADY_TIMESTEPS=0;
 				}
 				else if (CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
 					delta_volume = 0.0;
-					initial_volume = volA;
+					initial_volume = volA*Dm->Volume;
 					MORPH_ADAPT = false;
 					CURRENT_STEADY_TIMESTEPS=0;
 				}
@@ -781,6 +794,151 @@ void ScaLBL_ColorModel::Run(){
 
 	// ************************************************************************
 }
+
+double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){
+	
+	int nx = Nx;
+	int ny = Ny;
+	int nz = Nz;
+	int n;
+	int N = nx*ny*nz;
+	double volume_change=0.0;
+	
+	if (target_volume_change < 0.0){
+		Array<char> id_solid(nx,ny,nz);
+		Array<int> phase_label(nx,ny,nz);
+		DoubleArray distance(Nx,Ny,Nz);
+		DoubleArray phase(nx,ny,nz);
+		signed char *id_connected;
+		id_connected = new signed char [nx*ny*nz];
+
+		ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
+
+		// Extract only the connected part of NWP
+		BlobIDstruct new_index;
+		double vF=0.0; double vS=0.0;
+		ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm);
+		MPI_Barrier(Dm->Comm);
+
+		int count_oil=0;
+		int count_connected=0;
+		int count_porespace=0;
+		int count_water=0;
+		for (int k=1; k<nz-1; k++){
+			for (int j=1; j<ny-1; j++){
+				for (int i=1; i<nx-1; i++){
+					n=k*nx*ny+j*nx+i;
+					// only apply opening to connected component 
+					if ( phase_label(i,j,k) == 0){
+						count_connected++;
+					}
+					if (id[n] > 0){
+						count_porespace++;
+					}
+					if (id[n] == 2){
+						count_water++;
+					}
+				}
+			}
+		}
+		count_connected=sumReduce( Dm->Comm, count_connected);
+		count_porespace=sumReduce( Dm->Comm, count_porespace);
+		count_water=sumReduce( Dm->Comm, count_water);
+
+		for (int k=0; k<nz; k++){
+			for (int j=0; j<ny; j++){
+				for (int i=0; i<nx; i++){
+					n=k*nx*ny+j*nx+i;
+					// only apply opening to connected component 
+					if ( phase_label(i,j,k) == 0){
+						id_solid(i,j,k) = 1;
+						id_connected[n] = 2;
+						id[n] = 2;
+						/* delete the connected component */
+						phase(i,j,k) = -1.0;
+					}
+					else{
+						id_solid(i,j,k) = 0;
+						id_connected[n] = 0;
+					}
+				}
+			}
+		}
+		CalcDist(distance,id_solid,*Dm);
+
+		signed char water=2;
+		signed char notwater=1;
+		double SW=-(target_volume_change)/count_connected;
+		MorphOpen(distance, id_connected, Dm, SW, water, notwater);
+		
+		for (int k=0; k<nz; k++){
+			for (int j=0; j<ny; j++){
+				for (int i=0; i<nx; i++){
+					n=k*nx*ny+j*nx+i;
+					// only apply opening to connected component 
+					if ( id_connected[n] == 1){
+						id_solid(i,j,k) = 0;
+					}
+					else{
+						id_solid(i,j,k) = 1;
+					}
+				}
+			}
+		}
+		CalcDist(distance,id_solid,*Dm);
+		
+		// re-initialize
+		double beta = 0.95;
+		for (int k=0; k<nz; k++){
+			for (int j=0; j<ny; j++){
+				for (int i=0; i<nx; i++){
+					n=k*nx*ny+j*nx+i;
+					double d = distance(i,j,k);
+					if (Averages->SDs(i,j,k) > 0.f){
+						if (d < 3.f){
+							phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);	
+						}
+					}
+				} 
+			}
+		}
+		
+		int count_morphopen=0.0;
+		for (int k=1; k<nz-1; k++){
+			for (int j=1; j<ny-1; j++){
+				for (int i=1; i<nx-1; i++){
+					n=k*nx*ny+j*nx+i;
+					// only apply opening to connected component 
+					if ( id_connected[n] == 1){
+						count_morphopen++;
+					}
+				}
+			}
+		}
+		count_morphopen=sumReduce( Dm->Comm, count_morphopen);
+		volume_change = double(count_morphopen - count_connected);
+		
+		if (rank==0)  printf("   opening of connected oil %f \n",volume_change/count_connected);
+
+		ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
+		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){
+				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
+				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
+				ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
+			}
+			if (Dm->kproc() == nprocz-1){
+				ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
+				ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
+				ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
+			}
+		}
+	}
+	return(volume_change);
+}
+
 double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
 	srand(time(NULL));
 	double mass_loss =0.f;
@@ -812,12 +970,6 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){
 	if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
 	ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
 
-	FILE *OUTFILE;
-	sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
-	OUTFILE = fopen(LocalRankFilename,"wb");
-	fwrite(phase.data(),8,N,OUTFILE);
-	fclose(OUTFILE);
-
 	// 7. Re-initialize phase field and density
 	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);

models/ColorModel.h

@@ -81,6 +81,6 @@ private:
     void AssignComponentLabels(double *phase);
     double MorphInit(const double beta, const double morph_delta);
     double SeedPhaseField(const double seed_water_in_oil);
-    	
+    double MorphOpenConnected(double target_volume_change);
 };
 

models/MRTModel.cpp

@@ -246,6 +246,19 @@ void ScaLBL_MRTModel::Run(){
 			vay /= count;
 			vaz /= count;
 			
+			double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
+			double dir_x = Fx/force_mag;
+			double dir_y = Fy/force_mag;
+			double dir_z = Fz/force_mag;
+			if (force_mag == 0.0){
+				// default to z direction
+				dir_x = 0.0;
+				dir_y = 0.0;
+				dir_z = 1.0;
+				force_mag = 1.0;
+			}
+			double flow_rate = (vax*dir_x + vay*dir_y + vaz*dir_z);
+			
 			//if (rank==0) printf("Computing Minkowski functionals \n");
 			Morphology.ComputeScalar(Distance,0.f);
 			//Morphology.PrintAll();
@@ -260,7 +273,7 @@ void ScaLBL_MRTModel::Run(){
 			Xs=sumReduce( Dm->Comm, Xs);
 			if (rank==0) {
 				double h = Dm->voxel_length;
-				double absperm = h*h*mu*Mask->Porosity()*sqrt(vax*vax+vay*vay+vaz*vaz)/sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
+				double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
 				printf("     %f\n",absperm);
 				FILE * log_file = fopen("Permeability.csv","a");
 				fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu,