From 4090aa975a361783a87240fe0b66a9c29b380c36 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 12:36:57 -0400 Subject: [PATCH 01/23] print voxel length --- common/Domain.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/common/Domain.cpp b/common/Domain.cpp index 3d122f04..788f0ae5 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -122,6 +122,7 @@ void Domain::initialize( std::shared_ptr db ) voxel_length = 1.0; if (d_db->keyExists( "voxel_length" )){ auto voxel_length = d_db->getScalar("voxel_length"); + if (rank==0) printf("voxel length = %f micron \n", voxel_length); } if (d_db->keyExists( "InletLayers" )){ From ef968daedee42dc7d203d82c16303ccb6224b138 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 12:37:54 -0400 Subject: [PATCH 02/23] fix likely bug in steady timestep reset --- models/ColorModel.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index daf26b50..6379129a 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -689,6 +689,7 @@ void ScaLBL_ColorModel::Run(){ fclose(kr_log_file); printf(" Measured capillary number %f \n ",Ca); + CURRENT_STEADY_TIMESTEPS = 0; } if (SET_CAPILLARY_NUMBER ){ From cd904c9c91067ba4af264eea07a51178faad8fd9 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 13:14:44 -0400 Subject: [PATCH 03/23] helper function updates --- example/Workflow/HelperFunctions.R | 51 ++++++++++++++++++------------ models/MRTModel.cpp | 15 ++++++++- 2 files changed, 45 insertions(+), 21 deletions(-) diff --git a/example/Workflow/HelperFunctions.R b/example/Workflow/HelperFunctions.R index 7dd0170b..f8f5480b 100644 --- a/example/Workflow/HelperFunctions.R +++ b/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,26 @@ IngestRelperm<-function(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=FALSE,sep=" ") - colnames(D)<-c("time","nun","nuw","ift","Fx","Fy","Fz","Vn","Vw","unx","uny","unz","uwx","uwy","uwz") - D$Sw<-D$Vw/(D$Vn+D$Vw) - D$Krn<-D$Vn*D$nun*D$unx/D$Fx - D$Krw<-D$Vw*D$nuw*D$uwx/D$Fx - subset(D,D$time>100000) + 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) } @@ -39,18 +50,18 @@ ReadCase<-function(PATTERN){ } -D<-ReadCase("benth_w") -require("ggplot2") +ReadMorphDrain<-function(PATH){ + FILE=paste0(PATH,"/morphdrain.csv") + B<-read.csv(FILE,head=TRUE,sep=" ") + B$Media<-PATH -B<-read.csv("bentheimer/drain.csv",head=TRUE,sep=" ") -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))+ + theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+ + theme_bw() -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() + FILE=paste0(PATH,"-morphdrain.png") + ggsave(FILE,p,height=4.0,width=6.0) + return(B) +} -ggsave("morph-drain.png",p,height=4.0,width=6.0) \ No newline at end of file diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index d04009c1..5b0cdff3 100644 --- a/models/MRTModel.cpp +++ b/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, From 18d0affaaabab37d67f9245f2be429f2a80634c4 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 13:15:53 -0400 Subject: [PATCH 04/23] print voxel length when specified --- common/Domain.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 788f0ae5..723d8212 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -118,12 +118,6 @@ void Domain::initialize( std::shared_ptr db ) d_db = db; auto nproc = d_db->getVector("nproc"); auto n = d_db->getVector("n"); - - voxel_length = 1.0; - if (d_db->keyExists( "voxel_length" )){ - auto voxel_length = d_db->getScalar("voxel_length"); - if (rank==0) printf("voxel length = %f micron \n", voxel_length); - } if (d_db->keyExists( "InletLayers" )){ auto InletCount = d_db->getVector( "InletLayers" ); @@ -164,6 +158,12 @@ void Domain::initialize( std::shared_ptr 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("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("BC"); From 521520ecdcb906d694ca5d9f36fd3ae3c82e119b Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 13:16:03 -0400 Subject: [PATCH 05/23] align perm with force direction --- models/MRTModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 5b0cdff3..72b03925 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -273,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()*flow_rate / force_mag); + 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, From 82d0e18845ad454279db39eb931fc90d5d5c992f Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 16:32:51 -0400 Subject: [PATCH 06/23] working on morphopen connected --- analysis/SubPhase.cpp | 8 +-- models/ColorModel.cpp | 118 +++++++++++++++++++++++++++++++++++++++--- models/ColorModel.h | 2 +- 3 files changed, 116 insertions(+), 12 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 9c2f4529..b14c17ad 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -58,7 +58,7 @@ SubPhase::SubPhase(std::shared_ptr dm): } } - else{ +/* else{ char LocalRankString[8]; sprintf(LocalRankString,"%05d",Dm->rank()); char LocalRankFilename[40]; @@ -79,7 +79,7 @@ SubPhase::SubPhase(std::shared_ptr dm): fprintf(SUBPHASE,"Vi Ai Hi Xi "); // interface region fprintf(SUBPHASE,"Vic Aic Hic Xic Nic\n"); // interface region } - + */ if (Dm->rank()==0){ bool WriteHeader=false; TIMELOG = fopen("timelog.csv","r"); @@ -93,7 +93,7 @@ SubPhase::SubPhase(std::shared_ptr dm): { // If timelog is empty, write a short header to list the averages //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"sw krw krn qw qn pw pn\n"); + fprintf(TIMELOG,"sw krw krn qw qn pw pn force\n"); } } } @@ -254,7 +254,7 @@ void SubPhase::Basic(){ double krn = h*h*nu_n*not_water_flow_rate / force_mag ; double krw = h*h*nu_w*water_flow_rate / force_mag; //printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow); - fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*h*h*water_flow_rate,h*h*h*not_water_flow_rate, gwb.p, gnb.p); + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*h*h*water_flow_rate,h*h*h*not_water_flow_rate, gwb.p, gnb.p,force_mag); fflush(TIMELOG); } diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 6379129a..7f8d42df 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -727,7 +727,7 @@ void ScaLBL_ColorModel::Run(){ delta_volume = volA - 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 (rank==0) printf("***Morphological step with target volume change %f ***\n", delta_volume_target); @@ -782,6 +782,116 @@ 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; + Array id_solid(nx,ny,nz); + Array 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 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; kLastExterior(), 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 +922,6 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){ mass_loss= sumReduce( Dm->Comm, mass_loss); 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); diff --git a/models/ColorModel.h b/models/ColorModel.h index 25addc5e..c80a97a0 100644 --- a/models/ColorModel.h +++ b/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); }; From e2a7cca884322e275b7577ce7a7818d99864c129 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 16:37:17 -0400 Subject: [PATCH 07/23] fix domain voxel length bug --- common/Domain.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 723d8212..24c0c396 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -131,7 +131,11 @@ void Domain::initialize( std::shared_ptr db ) outlet_layers_y = OutletCount[1]; outlet_layers_z = OutletCount[2]; } - + voxel_length = 1.0; + if (d_db->keyExists( "voxel_length" )){ + auto voxel_length = d_db->getScalar("voxel_length"); + } + ASSERT( n.size() == 3u ); ASSERT( nproc.size() == 3u ); int nx = n[0]; @@ -158,11 +162,8 @@ void Domain::initialize( std::shared_ptr 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("voxel_length"); - if (myrank==0) printf("voxel length = %f micron \n", voxel_length); - } + + if (myrank==0) printf("voxel length = %f micron \n", voxel_length); id = new signed char[N]; memset(id,0,N); From 28b1e580c2731bda70df688bff0ca5cf110ac6ae Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 16:47:39 -0400 Subject: [PATCH 08/23] fix IO bug in subphase --- analysis/SubPhase.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index b14c17ad..d97fc724 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -93,7 +93,7 @@ SubPhase::SubPhase(std::shared_ptr dm): { // If timelog is empty, write a short header to list the averages //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"sw krw krn qw qn pw pn force\n"); + fprintf(TIMELOG,"sw krw krn qw qn pw pn\n"); } } } @@ -254,7 +254,7 @@ void SubPhase::Basic(){ double krn = h*h*nu_n*not_water_flow_rate / force_mag ; double krw = h*h*nu_w*water_flow_rate / force_mag; //printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow); - fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*h*h*water_flow_rate,h*h*h*not_water_flow_rate, gwb.p, gnb.p,force_mag); + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*h*h*water_flow_rate,h*h*h*not_water_flow_rate, gwb.p, gnb.p); fflush(TIMELOG); } From 3cdc9d84b67241c3a99fc55c47063e3240b7d540 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 1 May 2019 16:58:47 -0400 Subject: [PATCH 09/23] fix subphase bug --- analysis/SubPhase.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index d97fc724..9c2f4529 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -58,7 +58,7 @@ SubPhase::SubPhase(std::shared_ptr dm): } } -/* else{ + else{ char LocalRankString[8]; sprintf(LocalRankString,"%05d",Dm->rank()); char LocalRankFilename[40]; @@ -79,7 +79,7 @@ SubPhase::SubPhase(std::shared_ptr dm): fprintf(SUBPHASE,"Vi Ai Hi Xi "); // interface region fprintf(SUBPHASE,"Vic Aic Hic Xic Nic\n"); // interface region } - */ + if (Dm->rank()==0){ bool WriteHeader=false; TIMELOG = fopen("timelog.csv","r"); From 3957ea376fda693b05f09550d78703da71c6d95b Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 1 May 2019 21:16:13 -0400 Subject: [PATCH 10/23] fix likely mpi hang --- models/ColorModel.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 7f8d42df..1afeb161 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -689,9 +689,7 @@ void ScaLBL_ColorModel::Run(){ fclose(kr_log_file); printf(" Measured capillary number %f \n ",Ca); - CURRENT_STEADY_TIMESTEPS = 0; } - if (SET_CAPILLARY_NUMBER ){ Fx *= capillary_number / Ca; Fy *= capillary_number / Ca; @@ -710,6 +708,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){ From 0e8a36c45caae22eb29f446a6fe5376d5ca9a6f4 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 1 May 2019 21:44:49 -0400 Subject: [PATCH 11/23] clean up examples --- example/Bubble/input.db | 8 -------- example/Piston/Color.in | 6 ------ example/Piston/ComponentLabels.csv | 3 --- example/Piston/Domain.in | 4 ---- example/Piston/input.db | 8 -------- example/Plates/input.db | 8 -------- example/Poiseuille/Domain.in | 3 --- example/Poiseuille/Permeability.in | 4 ---- example/Sph125/Domain.in | 3 --- 9 files changed, 47 deletions(-) delete mode 100644 example/Piston/Color.in delete mode 100644 example/Piston/ComponentLabels.csv delete mode 100644 example/Piston/Domain.in delete mode 100644 example/Poiseuille/Domain.in delete mode 100644 example/Poiseuille/Permeability.in delete mode 100644 example/Sph125/Domain.in diff --git a/example/Bubble/input.db b/example/Bubble/input.db index 861788f0..3d07b2a5 100644 --- a/example/Bubble/input.db +++ b/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 } diff --git a/example/Piston/Color.in b/example/Piston/Color.in deleted file mode 100644 index e5a783f3..00000000 --- a/example/Piston/Color.in +++ /dev/null @@ -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 diff --git a/example/Piston/ComponentLabels.csv b/example/Piston/ComponentLabels.csv deleted file mode 100644 index a9bba3e2..00000000 --- a/example/Piston/ComponentLabels.csv +++ /dev/null @@ -1,3 +0,0 @@ -0 -0.5 -1 1.0 -2 -1.0 \ No newline at end of file diff --git a/example/Piston/Domain.in b/example/Piston/Domain.in deleted file mode 100644 index 641e6262..00000000 --- a/example/Piston/Domain.in +++ /dev/null @@ -1,4 +0,0 @@ -1 1 10 -40 40 40 -229 -1.0 1.0 1.0 diff --git a/example/Piston/input.db b/example/Piston/input.db index 82dbf077..fab67cc5 100644 --- a/example/Piston/input.db +++ b/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" diff --git a/example/Plates/input.db b/example/Plates/input.db index 3253c993..2e74a43f 100644 --- a/example/Plates/input.db +++ b/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" diff --git a/example/Poiseuille/Domain.in b/example/Poiseuille/Domain.in deleted file mode 100644 index 8a657dff..00000000 --- a/example/Poiseuille/Domain.in +++ /dev/null @@ -1,3 +0,0 @@ -1 1 1 -100 100 100 -1.0 1.0 1.0 diff --git a/example/Poiseuille/Permeability.in b/example/Poiseuille/Permeability.in deleted file mode 100644 index 88640c05..00000000 --- a/example/Poiseuille/Permeability.in +++ /dev/null @@ -1,4 +0,0 @@ -1.0 -0.0 0.0 1.0e-5 -0 0 1.0 1.0 -1000 1000 1.0e-5 diff --git a/example/Sph125/Domain.in b/example/Sph125/Domain.in deleted file mode 100644 index 55a65a55..00000000 --- a/example/Sph125/Domain.in +++ /dev/null @@ -1,3 +0,0 @@ -2 2 2 -100 100 100 -1.0 1.0 1.0 From fd713ba3af9b218eb14ac951c088cab82762f72a Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Thu, 2 May 2019 15:25:43 -0400 Subject: [PATCH 12/23] fix voxel length --- common/Domain.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 24c0c396..75f2082b 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -133,7 +133,7 @@ void Domain::initialize( std::shared_ptr db ) } voxel_length = 1.0; if (d_db->keyExists( "voxel_length" )){ - auto voxel_length = d_db->getScalar("voxel_length"); + voxel_length = d_db->getScalar("voxel_length"); } ASSERT( n.size() == 3u ); From 89461a265c966bb3feb17a41a2370396bf3c20bb Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 2 May 2019 21:16:12 -0400 Subject: [PATCH 13/23] adjust kr computation --- analysis/SubPhase.cpp | 4 ++-- models/ColorModel.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 9c2f4529..3047bd0a 100644 --- a/analysis/SubPhase.cpp +++ b/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 not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.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*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; diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 1afeb161..0e39ee73 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -641,8 +641,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_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_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*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 ){ From 8822d5036d23c40ed7870aafa12d5ee1c9871d1b Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 2 May 2019 21:51:02 -0400 Subject: [PATCH 14/23] added morphopen connected oil option --- models/ColorModel.cpp | 199 ++++++++++++++++++++++-------------------- 1 file changed, 103 insertions(+), 96 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 0e39ee73..27f4b8c0 100644 --- a/models/ColorModel.cpp +++ 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; @@ -492,6 +493,10 @@ void ScaLBL_ColorModel::Run(){ morph_interval = analysis_db->getScalar( "morph_interval" ); USE_MORPH = true; } + if (analysis_db->keyExists( "use_morphopen_oil" )){ + USE_MORPHOPEN_OIL = analysis_db->getScalar( "use_morphopen_oil" ); + USE_MORPH = true; + } if (analysis_db->keyExists( "tolerance" )){ tolerance = analysis_db->getScalar( "tolerance" ); } @@ -656,7 +661,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 ); @@ -723,11 +728,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, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target); } + else if (USE_MORPHOPEN_OIL){ + if (rank==0) printf("***Morphological opening of connected oil, with target volume change ***\n", delta_volume_target); + MorphologicalOpening(delta_volume_target); + } else{ if (rank==0) printf("***Morphological step with target volume change %f ***\n", delta_volume_target); //double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A @@ -735,13 +744,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; } @@ -789,105 +798,103 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ int nz = Nz; int n; int N = nx*ny*nz; - double volume_change; - Array id_solid(nx,ny,nz); - Array 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]; + double volume_change=0.0; - ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); + if (target_volume_change > 0.0){ + Array id_solid(nx,ny,nz); + Array 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]; - // 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); + 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 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 0){ - count_porespace++; - } - if (id[n] == 2){ - count_water++; - } + 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); } - } - } - 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; kkproc() == 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); } } } - CalcDist(distance,id_solid,*Dm); - - double SW=0.5; - // target water increase in voxels, normalized by connected volume - double St = (SW*count_porespace - count_water)/count_porespace; - - signed char water=2; - signed char notwater=1; - MorphOpen(distance, id_connected, Dm, St, water, notwater); - - for (int k=0; kLastExterior(), 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); } From 1b8e6cbbd5177ac547c431342baff1f118c3cb03 Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Thu, 2 May 2019 21:52:51 -0400 Subject: [PATCH 15/23] fix open bug --- example/Workflow/HelperFunctions.R | 2 ++ models/ColorModel.cpp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/example/Workflow/HelperFunctions.R b/example/Workflow/HelperFunctions.R index f8f5480b..57a7db11 100644 --- a/example/Workflow/HelperFunctions.R +++ b/example/Workflow/HelperFunctions.R @@ -17,12 +17,14 @@ ReadSubphase<-function(PATH){ 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=TRUE,sep=" ") + D$Case<-PATH p<-ggplot(D)+ geom_line(aes(sat.water,eff.perm.oil,color="oil"))+ diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 27f4b8c0..d2a1cc7f 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -735,7 +735,7 @@ void ScaLBL_ColorModel::Run(){ } else if (USE_MORPHOPEN_OIL){ if (rank==0) printf("***Morphological opening of connected oil, with target volume change ***\n", delta_volume_target); - MorphologicalOpening(delta_volume_target); + MorphOpenConnected(delta_volume_target); } else{ if (rank==0) printf("***Morphological step with target volume change %f ***\n", delta_volume_target); From b5225efb92d522a9b81de15f6f9f89507808f530 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 3 May 2019 14:26:36 -0400 Subject: [PATCH 16/23] morph open connected --- models/ColorModel.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index d2a1cc7f..a1226741 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -867,6 +867,7 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ double SW=-(target_volume_change)/count_connected; MorphOpen(distance, id_connected, Dm, SW, water, notwater); + int count_morphopen=0.0; for (int k=0; kComm, count_morphopen); + volume_change = double(count_morphopen - count_connected); + if (rank==0) printf(" opening of connected oil %f \n",volume_change); + 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); From 9c6e853be44cd8a638a25288387dcc08c0d20343 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 3 May 2019 16:31:12 -0400 Subject: [PATCH 17/23] volume change by morphological opening --- models/ColorModel.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index a1226741..f18e7b28 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -484,10 +484,12 @@ void ScaLBL_ColorModel::Run(){ } if (analysis_db->keyExists( "seed_water" )){ seed_water = analysis_db->getScalar( "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( "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( "morph_interval" ); @@ -495,6 +497,7 @@ void ScaLBL_ColorModel::Run(){ } if (analysis_db->keyExists( "use_morphopen_oil" )){ USE_MORPHOPEN_OIL = analysis_db->getScalar( "use_morphopen_oil" ); + if (rank == 0 && USE_MORPHOPEN_OIL) printf("Volume change by morphological opening \n"); USE_MORPH = true; } if (analysis_db->keyExists( "tolerance" )){ @@ -734,6 +737,7 @@ void ScaLBL_ColorModel::Run(){ 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 ***\n", delta_volume_target); MorphOpenConnected(delta_volume_target); } @@ -800,7 +804,7 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ int N = nx*ny*nz; double volume_change=0.0; - if (target_volume_change > 0.0){ + if (target_volume_change < 0.0){ Array id_solid(nx,ny,nz); Array phase_label(nx,ny,nz); DoubleArray distance(Nx,Ny,Nz); From 5bb2bcb62a3ae243de0363358da864de2a049cf7 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 3 May 2019 17:43:13 -0400 Subject: [PATCH 18/23] morph open connected --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index f18e7b28..8b29785a 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -738,7 +738,7 @@ void ScaLBL_ColorModel::Run(){ } 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 ***\n", delta_volume_target); + if (rank==0) printf("***Morphological opening of connected oil, with target volume change %f ***\n", delta_volume_target); MorphOpenConnected(delta_volume_target); } else{ @@ -869,7 +869,7 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ signed char water=2; signed char notwater=1; double SW=-(target_volume_change)/count_connected; - MorphOpen(distance, id_connected, Dm, SW, water, notwater); + MorphOpen(distance, id_connected, Dm, SW, notwater, water); int count_morphopen=0.0; for (int k=0; k Date: Fri, 3 May 2019 19:23:54 -0400 Subject: [PATCH 19/23] morph open connected --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8b29785a..4cbc3e2a 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -869,7 +869,7 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ signed char water=2; signed char notwater=1; double SW=-(target_volume_change)/count_connected; - MorphOpen(distance, id_connected, Dm, SW, notwater, water); + MorphOpen(distance, id_connected, Dm, SW, water, notwater); int count_morphopen=0.0; for (int k=0; kComm, count_morphopen); volume_change = double(count_morphopen - count_connected); - if (rank==0) printf(" opening of connected oil %f \n",volume_change); + 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); From b28f0e642932e8085708b31bc10315ff0e82d481 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 3 May 2019 20:39:20 -0400 Subject: [PATCH 20/23] remove extra printing --- analysis/morphology.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 4b4e9a82..ebc7fc0c 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -279,10 +279,10 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr 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 Date: Fri, 3 May 2019 20:47:48 -0400 Subject: [PATCH 21/23] morph open connected --- models/ColorModel.cpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 4cbc3e2a..fb639bc1 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -870,15 +870,26 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ signed char notwater=1; double SW=-(target_volume_change)/count_connected; MorphOpen(distance, id_connected, Dm, SW, water, notwater); - - int count_morphopen=0.0; + for (int k=0; k Date: Fri, 3 May 2019 20:55:37 -0400 Subject: [PATCH 22/23] morph open connected --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index fb639bc1..fd7bc935 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -877,7 +877,7 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ n=k*nx*ny+j*nx+i; // only apply opening to connected component if ( id_connected[n] == 1){ - count_morphopen++; + phase(i,j,k) = 1.0; } } } From 6ff2f6ce1d60b7538f7935f0a2f2a98a2c9de3b7 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 3 May 2019 21:54:50 -0400 Subject: [PATCH 23/23] use distance for re-init --- models/ColorModel.cpp | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index fd7bc935..d73cbee8 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -877,12 +877,32 @@ double ScaLBL_ColorModel::MorphOpenConnected(double target_volume_change){ n=k*nx*ny+j*nx+i; // only apply opening to connected component if ( id_connected[n] == 1){ - phase(i,j,k) = 1.0; + 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; kSDs(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