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

common/Domain.cpp

@@ -133,11 +133,6 @@ void Domain::initialize( std::shared_ptr<Database> db )
     d_db = 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" );
@@ -178,6 +173,12 @@ 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;
+    voxel_length = 1.0;
+    if (d_db->keyExists( "voxel_length" )){
+    	auto voxel_length = d_db->getScalar<double>("voxel_length");
+    	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/Workflow/HelperFunctions.R

@@ -0,0 +1,67 @@
+require("ggplot2")
+
+ReadSubphase<-function(PATH){
+  FILE=paste0(PATH,"/subphase.csv")
+  S<-read.csv(FILE,head=TRUE,sep=" ")
+  S$Vw<-S$Vwc+S$Vwd
+  S$Vn<-S$Vnc+S$Vnd
+  S$Sw<-S$Vw/(S$Vn+S$Vw)
+  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$Case<-PATH
+  return(S)
+}
+
+ReadTimelog<-function(PATH){
+	FILE=paste0(PATH,"/timelog.csv")
+	D<-read.csv(file=FILE,head=TRUE,sep=" ")
+	return(D)
+}
+
+ReadRelperm<-function(PATH){
+	FILE=paste0(PATH,"/relperm.csv")
+	D<-read.csv(file=FILE,head=TRUE,sep=" ")
+
+	p<-ggplot(D)+
+	  geom_line(aes(sat.water,eff.perm.oil,color="oil"))+
+	  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)
+}
+
+ReadCase<-function(PATTERN){
+    list<-list.files(pattern=PATTERN)
+    Data=NULL
+    for (k in 1:length(list)){
+        print(list[k])
+        tmp=ReadRelperm(list[k])
+	tmp$Case<-list[k]
+        Data<-rbind(Data,tmp)
+    }
+    return(Data)
+}
+
+
+ReadMorphDrain<-function(PATH){
+	FILE=paste0(PATH,"/morphdrain.csv")
+	B<-read.csv(FILE,head=TRUE,sep=" ")
+	B$Media<-PATH
+
+	p<-ggplot()+
+	  geom_line(data=B,aes(Sw,2/R,colour=Media))+
+	  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+
+	  theme_bw()
+
+	FILE=paste0(PATH,"-morphdrain.png")
+ 	ggsave(FILE,p,height=4.0,width=6.0)
+	return(B)
+}
+

models/ColorModel.cpp

@@ -704,6 +704,7 @@ void ScaLBL_ColorModel::Run(){
 						fclose(kr_log_file);
 
 						printf("  Measured capillary number %f \n ",Ca);
+						CURRENT_STEADY_TIMESTEPS = 0;
 					}
 
 					if (SET_CAPILLARY_NUMBER ){

models/MRTModel.cpp

@@ -261,6 +261,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();
@@ -275,7 +288,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,