From b82b81dfff471a054ae739c33fa5109bd8b1d5b0 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 14:31:17 -0400 Subject: [PATCH 01/79] skeleton code for morph lbm --- models/ColorModel.cpp | 89 ++++++++++++++++++++++++++++++++++++++++++- models/ColorModel.h | 1 + 2 files changed, 89 insertions(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 3bbd7d91..c2dca2e5 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -326,7 +326,7 @@ void ScaLBL_ColorModel::Initialize(){ cDen = new double[2*Np]; cDist = new double[19*Np]; ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int)); - ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double)); + ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double)); ifstream File(LocalRestartFile,ios::binary); int idx; @@ -520,6 +520,93 @@ void ScaLBL_ColorModel::Run(){ // ************************************************************************ } +void ColorModel::MorphInit(double beta, double morph_delta){ + + double vF = 0.f; + double vS = 0.f; + + DoubleArray phase(Nx,Ny,Nz); + IntArray phase_label(Nx,Ny,Nz);; + DoubleArray phase_distance(Nx,Ny,Nz); + + // Basic algorithm to + // 1. Copy phase field to CPU + ScaLBL_CopyToHost(phase, Phi, N*sizeof(double)); + + // 2. Identify connected components of phase field -> phase_label + BlobIDstruct new_index; + new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); + // only operate on component "0" + for (int k=0; k phase_distance + CalcDist(phase_distance,phase_label,*Dm); + + double factor=0.5/Beta; + for (int k=0; kSDs(i,j,k) > 0.f){ + // only update phase field in immediate proximity of largest component + if (d < 3.f){ + phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d)) - 1.f); + } + } + } + } + } + // 6. copy back to the device + ScaLBL_CopyToDevice(Phi,phase,N*sizeof(double)); + + // 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); + 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); + } + } + +} + void ScaLBL_ColorModel::WriteDebug(){ // Copy back final phase indicator field and convert to regular layout DoubleArray PhaseField(Nx,Ny,Nz); diff --git a/models/ColorModel.h b/models/ColorModel.h index 4a43de66..8c40a22a 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -77,6 +77,7 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr db0); void AssignComponentLabels(double *phase); + void MorphInit(double beta, double morph_delta); }; From 26216b23ad6dbc3720ad40c96cf82a254fcae85c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 14:49:34 -0400 Subject: [PATCH 02/79] morphinit builds --- analysis/morphology.cpp | 2 ++ models/ColorModel.cpp | 18 ++++++++++-------- 2 files changed, 12 insertions(+), 8 deletions(-) create mode 100644 analysis/morphology.cpp diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp new file mode 100644 index 00000000..9302cd9a --- /dev/null +++ b/analysis/morphology.cpp @@ -0,0 +1,2 @@ +#include "analysis/morphology.h" + diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index c2dca2e5..65827cf5 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -520,7 +520,7 @@ void ScaLBL_ColorModel::Run(){ // ************************************************************************ } -void ColorModel::MorphInit(double beta, double morph_delta){ +void ScaLBL_ColorModel::MorphInit(double beta, double morph_delta){ double vF = 0.f; double vS = 0.f; @@ -528,28 +528,30 @@ void ColorModel::MorphInit(double beta, double morph_delta){ DoubleArray phase(Nx,Ny,Nz); IntArray phase_label(Nx,Ny,Nz);; DoubleArray phase_distance(Nx,Ny,Nz); + Array phase_id(Nx,Ny,Nz); // Basic algorithm to // 1. Copy phase field to CPU - ScaLBL_CopyToHost(phase, Phi, N*sizeof(double)); + ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); // 2. Identify connected components of phase field -> phase_label BlobIDstruct new_index; - new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); + new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); // only operate on component "0" for (int k=0; k phase_distance - CalcDist(phase_distance,phase_label,*Dm); + CalcDist(phase_distance,phase_id,*Dm); - double factor=0.5/Beta; + double temp,value; + double factor=0.5/beta; for (int k=0; kLastExterior(), Np); From 59b2667188109d16c4379df7a6a034a3bb684f2e Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 14:52:37 -0400 Subject: [PATCH 03/79] test color model morph --- models/ColorModel.cpp | 8 +++++++- models/ColorModel.h | 2 +- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 65827cf5..8054b8dd 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -495,6 +495,12 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); + + if (timestep%2000 == 0){ + double morph_delta=0.5; + MorphInit(beta,morph_delta); + MPI_Barrier(comm); + } } analysis.finish(); @@ -520,7 +526,7 @@ void ScaLBL_ColorModel::Run(){ // ************************************************************************ } -void ScaLBL_ColorModel::MorphInit(double beta, double morph_delta){ +void ScaLBL_ColorModel::MorphInit(const double beta, const 2,double morph_delta){ double vF = 0.f; double vS = 0.f; diff --git a/models/ColorModel.h b/models/ColorModel.h index 8c40a22a..5b5e9ccf 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -77,7 +77,7 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr db0); void AssignComponentLabels(double *phase); - void MorphInit(double beta, double morph_delta); + void MorphInit(const double beta, const double morph_delta); }; From 30e5ffebc616a312978bbafca23ca9ad4b599446 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 14:54:01 -0400 Subject: [PATCH 04/79] remove typoe --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8054b8dd..63a2d90e 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -526,7 +526,7 @@ void ScaLBL_ColorModel::Run(){ // ************************************************************************ } -void ScaLBL_ColorModel::MorphInit(const double beta, const 2,double morph_delta){ +void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ double vF = 0.f; double vS = 0.f; From 7cdebdd181661f9703033017bac7967c6dfccf52 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 15:44:44 -0400 Subject: [PATCH 05/79] morphological intial condition runs --- models/ColorModel.cpp | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 63a2d90e..08c0e92b 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -496,7 +496,7 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); - if (timestep%2000 == 0){ + if (timestep%100 == 0){ double morph_delta=0.5; MorphInit(beta,morph_delta); MPI_Barrier(comm); @@ -527,7 +527,8 @@ void ScaLBL_ColorModel::Run(){ } void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ - + const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + double vF = 0.f; double vS = 0.f; @@ -538,11 +539,15 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ // Basic algorithm to // 1. Copy phase field to CPU + if (rank==0) printf("MorphInit: copy data \n"); ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); + if (rank==0) printf("MorphInit: get blob ids \n"); // 2. Identify connected components of phase field -> phase_label BlobIDstruct new_index; - new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); + ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); + MPI_Barrier(comm); + if (rank==0) printf("MorphInit: label largest feature \n"); // only operate on component "0" for (int k=0; k phase_distance CalcDist(phase_distance,phase_id,*Dm); @@ -571,6 +577,7 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ } // 4. Apply erosion / dilation operation to phase_distance + if (rank==0) printf("MorphInit: morphological operation \n"); for (int k=0; kLastExterior(), Np); ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); if (BoundaryCondition >0 ){ From 154bb45fee6e010d1b66445511a1c9c6d24ed9f0 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 15:49:16 -0400 Subject: [PATCH 06/79] remove debug write statements --- tests/lbpm_color_simulator.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index eb9d8f62..18ac9a80 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -54,7 +54,7 @@ int main(int argc, char **argv) ColorModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables ColorModel.Initialize(); // initializing the model will set initial conditions for variables ColorModel.Run(); - ColorModel.WriteDebug(); + //ColorModel.WriteDebug(); PROFILE_STOP("Main"); PROFILE_SAVE("lbpm_color_simulator",1); From 5c409e54426c3366958312f9d2d0c6c7e7b0e814 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 16:44:36 -0400 Subject: [PATCH 07/79] morph init getting closer not quite there --- models/ColorModel.cpp | 32 ++++++++++++++++++++++++-------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 08c0e92b..d36771e8 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -553,8 +553,8 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ for (int j=0; j 1.f) value=1.f; + if (value < -1.f) value=-1.f; // temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm. - temp = factor*log((1.0+value)/(1.0-value)); + temp = -factor*log((1.0+value)/(1.0-value)); /// use this approximation close to the object - if (phase_distance(i,j,k) < 1.f) phase_distance(i,j,k) = temp; + if (Averages->SDs(i,j,k) < 0.f ) phase_distance(i,j,k) = (-1.f)*Averages->SDs(i,j,k); + else if (fabs(phase_distance(i,j,k)) < 2.f ) phase_distance(i,j,k) = temp; } } } - + // 4. Apply erosion / dilation operation to phase_distance if (rank==0) printf("MorphInit: morphological operation \n"); for (int k=0; kSDs(i,j,k) > 0.f){ // only update phase field in immediate proximity of largest component - if (d < 3.f){ - phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d)) - 1.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); + } } } } } + FILE *OUTFILE; + sprintf(LocalRankFilename,"Phase.%05i.raw",rank); + OUTFILE = fopen(LocalRankFilename,"wb"); + fwrite(phase.data(),8,N,OUTFILE); + fclose(OUTFILE); + + + FILE *DISTFILE; + sprintf(LocalRankFilename,"Distance.%05i.raw",rank); + DISTFILE = fopen(LocalRankFilename,"wb"); + fwrite(phase_distance.data(),8,N,DISTFILE); + fclose(DISTFILE); + // 6. copy back to the device if (rank==0) printf("MorphInit: copy data back to device\n"); ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double)); From bce4069111b66320c8478d0df2615bcf131078fc Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 16:51:38 -0400 Subject: [PATCH 08/79] fixed indentation --- models/ColorModel.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index d36771e8..9fa20978 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -527,11 +527,11 @@ void ScaLBL_ColorModel::Run(){ } void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ - const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); double vF = 0.f; double vS = 0.f; - + DoubleArray phase(Nx,Ny,Nz); IntArray phase_label(Nx,Ny,Nz);; DoubleArray phase_distance(Nx,Ny,Nz); @@ -544,9 +544,9 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ if (rank==0) printf("MorphInit: get blob ids \n"); // 2. Identify connected components of phase field -> phase_label - BlobIDstruct new_index; - ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); - MPI_Barrier(comm); + BlobIDstruct new_index; + ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); + MPI_Barrier(comm); if (rank==0) printf("MorphInit: label largest feature \n"); // only operate on component "0" for (int k=0; k phase_distance CalcDist(phase_distance,phase_id,*Dm); - + double temp,value; double factor=0.5/beta; for (int k=0; kSDs(i,j,k) > 0.f){ // only update phase field in immediate proximity of largest component - if (d < 3.f){ - phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.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); + } } } } @@ -622,7 +622,7 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ // 6. copy back to the device if (rank==0) printf("MorphInit: copy data back to device\n"); ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double)); - + // 7. Re-initialize phase field and density if (rank==0) printf("MorphInit: re-initialize LBM \n"); From 4933f846328ec9061ba8e327a1e28ce088298285 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 23:31:53 -0400 Subject: [PATCH 09/79] tweak distance in solid --- models/ColorModel.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 9fa20978..73466607 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -573,8 +573,12 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ // temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm. temp = -factor*log((1.0+value)/(1.0-value)); /// use this approximation close to the object - if (Averages->SDs(i,j,k) < 0.f ) phase_distance(i,j,k) = (-1.f)*Averages->SDs(i,j,k); - else if (fabs(phase_distance(i,j,k)) < 2.f ) phase_distance(i,j,k) = temp; + if (fabs(phase_distance(i,j,k)) < 2.f ){ + if (Averages->SDs(i,j,k) < 0.f ) + phase_distance(i,j,k) = (-1.f)*Averages->SDs(i,j,k); + else + phase_distance(i,j,k) = temp; + } } } } From b260e875bb8659a1d3d95075921693eeb2f80d41 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Oct 2018 23:34:40 -0400 Subject: [PATCH 10/79] tweak distance in solid --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 73466607..c5afc893 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -574,7 +574,7 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ temp = -factor*log((1.0+value)/(1.0-value)); /// use this approximation close to the object if (fabs(phase_distance(i,j,k)) < 2.f ){ - if (Averages->SDs(i,j,k) < 0.f ) + if (Averages->SDs(i,j,k) < 0.f && fabs(Averages->SDs(i,j,k)) > phase_distance(i,j,k) ) phase_distance(i,j,k) = (-1.f)*Averages->SDs(i,j,k); else phase_distance(i,j,k) = temp; From e6e3316564e339f9b625909d226bddd076fb0a91 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 23 Oct 2018 00:43:59 -0400 Subject: [PATCH 11/79] morphological operation functioning in tests --- models/ColorModel.cpp | 69 ++++++++++++++++++++++++++++++++----------- 1 file changed, 51 insertions(+), 18 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index c5afc893..72be7f80 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -539,15 +539,34 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ // Basic algorithm to // 1. Copy phase field to CPU - if (rank==0) printf("MorphInit: copy data \n"); + //if (rank==0) printf("MorphInit: copy data \n"); ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); - if (rank==0) printf("MorphInit: get blob ids \n"); + double count,count_global,volume_initial,volume_final; + count = 0.f; + for (int k=0; k 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f; + } + } + } + MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm); + volume_initial = count_global; + + //if (rank==0) printf("MorphInit: get blob ids \n"); // 2. Identify connected components of phase field -> phase_label BlobIDstruct new_index; ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); MPI_Barrier(comm); - if (rank==0) printf("MorphInit: label largest feature \n"); + + /* FILE *IDFILE; + sprintf(LocalRankFilename,"Label.%05i.raw",rank); + IDFILE = fopen(LocalRankFilename,"wb"); + fwrite(phase_label.data(),4,N,IDFILE); + fclose(IDFILE); + */ + //if (rank==0) printf("MorphInit: label largest feature \n"); // only operate on component "0" for (int k=0; k phase_distance CalcDist(phase_distance,phase_id,*Dm); - + double temp,value; double factor=0.5/beta; for (int k=0; k 1.f) value=1.f; if (value < -1.f) value=-1.f; // temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm. temp = -factor*log((1.0+value)/(1.0-value)); /// use this approximation close to the object - if (fabs(phase_distance(i,j,k)) < 2.f ){ - if (Averages->SDs(i,j,k) < 0.f && fabs(Averages->SDs(i,j,k)) > phase_distance(i,j,k) ) - phase_distance(i,j,k) = (-1.f)*Averages->SDs(i,j,k); - else - phase_distance(i,j,k) = temp; + if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){ + phase_distance(i,j,k) = temp; } + } } } } - + // 4. Apply erosion / dilation operation to phase_distance - if (rank==0) printf("MorphInit: morphological operation \n"); + //if (rank==0) printf("MorphInit: morphological operation \n"); for (int k=0; kSDs(i,j,k); + double wallweight = 1.f / (1+exp(-5.f*(walldist-1.f))); + phase_distance(i,j,k) -= wallweight*morph_delta; } } } - if (rank==0) printf("MorphInit: reinitialize phase field \n"); + //if (rank==0) printf("MorphInit: reinitialize phase field \n"); // 5. Update phase indicator field based on new distnace for (int k=0; k 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f; + } + } + } + MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm); + volume_final=count_global; + + if (rank == 0) printf("Morphological operation change volume fraction by %f \n", (volume_final-volume_initial)/volume_initial); + + + /* FILE *OUTFILE; sprintf(LocalRankFilename,"Phase.%05i.raw",rank); OUTFILE = fopen(LocalRankFilename,"wb"); fwrite(phase.data(),8,N,OUTFILE); fclose(OUTFILE); - FILE *DISTFILE; sprintf(LocalRankFilename,"Distance.%05i.raw",rank); DISTFILE = fopen(LocalRankFilename,"wb"); fwrite(phase_distance.data(),8,N,DISTFILE); fclose(DISTFILE); - + */ // 6. copy back to the device - if (rank==0) printf("MorphInit: copy data back to device\n"); + //if (rank==0) printf("MorphInit: copy data back to device\n"); ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double)); // 7. Re-initialize phase field and density From 2b23951bcee4529b8560bbc44345da8dca6d6405 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 23 Oct 2018 00:50:07 -0400 Subject: [PATCH 12/79] clean up code --- models/ColorModel.cpp | 78 ++++++++++++++----------------------------- 1 file changed, 25 insertions(+), 53 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 72be7f80..4e67970d 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -539,7 +539,6 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ // Basic algorithm to // 1. Copy phase field to CPU - //if (rank==0) printf("MorphInit: copy data \n"); ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); double count,count_global,volume_initial,volume_final; @@ -547,26 +546,17 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ for (int k=0; k 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f; + if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f; } } } MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm); volume_initial = count_global; - //if (rank==0) printf("MorphInit: get blob ids \n"); // 2. Identify connected components of phase field -> phase_label BlobIDstruct new_index; ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); MPI_Barrier(comm); - - /* FILE *IDFILE; - sprintf(LocalRankFilename,"Label.%05i.raw",rank); - IDFILE = fopen(LocalRankFilename,"wb"); - fwrite(phase_label.data(),4,N,IDFILE); - fclose(IDFILE); - */ - //if (rank==0) printf("MorphInit: label largest feature \n"); // only operate on component "0" for (int k=0; k phase_distance CalcDist(phase_distance,phase_id,*Dm); - + double temp,value; double factor=0.5/beta; for (int k=0; k 1.f) value=1.f; - if (value < -1.f) value=-1.f; - // temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm. - temp = -factor*log((1.0+value)/(1.0-value)); - /// use this approximation close to the object - if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){ - phase_distance(i,j,k) = temp; + if (phase_distance(i,j,k) < 3.f ){ + value = phase(i,j,k); + if (value > 1.f) value=1.f; + if (value < -1.f) value=-1.f; + // temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm. + temp = -factor*log((1.0+value)/(1.0-value)); + /// use this approximation close to the object + if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){ + phase_distance(i,j,k) = temp; + } } - } - } - } - } - - // 4. Apply erosion / dilation operation to phase_distance - //if (rank==0) printf("MorphInit: morphological operation \n"); - for (int k=0; kSDs(i,j,k); - double wallweight = 1.f / (1+exp(-5.f*(walldist-1.f))); - phase_distance(i,j,k) -= wallweight*morph_delta; } } } - //if (rank==0) printf("MorphInit: reinitialize phase field \n"); + // 4. Apply erosion / dilation operation to phase_distance + for (int k=0; kSDs(i,j,k); + double wallweight = 1.f / (1+exp(-5.f*(walldist-1.f))); + phase_distance(i,j,k) -= wallweight*morph_delta; + } + } + } + // 5. Update phase indicator field based on new distnace for (int k=0; k 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f; + if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f; } } } MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm); volume_final=count_global; - if (rank == 0) printf("Morphological operation change volume fraction by %f \n", (volume_final-volume_initial)/volume_initial); - + if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", (volume_final-volume_initial)/volume_initial); - /* FILE *OUTFILE; - sprintf(LocalRankFilename,"Phase.%05i.raw",rank); - OUTFILE = fopen(LocalRankFilename,"wb"); - fwrite(phase.data(),8,N,OUTFILE); - fclose(OUTFILE); - - FILE *DISTFILE; - sprintf(LocalRankFilename,"Distance.%05i.raw",rank); - DISTFILE = fopen(LocalRankFilename,"wb"); - fwrite(phase_distance.data(),8,N,DISTFILE); - fclose(DISTFILE); - */ // 6. copy back to the device //if (rank==0) printf("MorphInit: copy data back to device\n"); ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double)); // 7. Re-initialize phase field and density - if (rank==0) printf("MorphInit: re-initialize LBM \n"); - 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 ){ From 3e63810065380d3ea52cba47375c8b1616ff8a60 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 23 Oct 2018 15:38:01 -0400 Subject: [PATCH 13/79] morph stuff in input file --- models/ColorModel.cpp | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 4e67970d..ae79158b 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -394,20 +394,35 @@ void ScaLBL_ColorModel::Initialize(){ void ScaLBL_ColorModel::Run(){ int nprocs=nprocx*nprocy*nprocz; const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + + int morph_interval; + double morph_delta; + if (analysis_db->keyExists( "morph_delta" )){ + morph_delta = domain_db->getVector( "morph_delta" ); + } + else{ + morph_delta=0.5; + } + if (analysis_db->keyExists( "morph_interval" )){ + morph_interval = domain_db->getVector( "morph_interval" ); + } + else{ + morph_interval=100000; + } + if (rank==0){ printf("********************************************************\n"); printf("No. of timesteps: %i \n", timestepMax); fflush(stdout); } - //.......create and start timer............ double starttime,stoptime,cputime; ScaLBL_DeviceBarrier(); MPI_Barrier(comm); starttime = MPI_Wtime(); //......................................... - + //************ MAIN ITERATION LOOP ***************************************/ PROFILE_START("Loop"); //std::shared_ptr analysis_db; @@ -496,8 +511,7 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); - if (timestep%100 == 0){ - double morph_delta=0.5; + if (timestep%morph_interval == 0){ MorphInit(beta,morph_delta); MPI_Barrier(comm); } From 48b39580c710d3ff0bd86977e03c8353779bcfb7 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 23 Oct 2018 23:51:45 -0400 Subject: [PATCH 14/79] remove morphology.cpp --- analysis/morphology.cpp | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 analysis/morphology.cpp diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp deleted file mode 100644 index 9302cd9a..00000000 --- a/analysis/morphology.cpp +++ /dev/null @@ -1,2 +0,0 @@ -#include "analysis/morphology.h" - From 760c1bee2e29ef08d6b783e70b6465095fe450b1 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 23 Oct 2018 23:59:38 -0400 Subject: [PATCH 15/79] fix morph input --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index ae79158b..f57b714e 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -398,13 +398,13 @@ void ScaLBL_ColorModel::Run(){ int morph_interval; double morph_delta; if (analysis_db->keyExists( "morph_delta" )){ - morph_delta = domain_db->getVector( "morph_delta" ); + morph_delta = analysis_db->getScalar( "morph_delta" ); } else{ morph_delta=0.5; } if (analysis_db->keyExists( "morph_interval" )){ - morph_interval = domain_db->getVector( "morph_interval" ); + morph_interval = analysis_db->getScalar( "morph_interval" ); } else{ morph_interval=100000; From 2a6c1a88947fed419ca94525a939885ad4809544 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 24 Oct 2018 11:10:45 -0400 Subject: [PATCH 16/79] adding target Ca --- models/ColorModel.cpp | 41 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index f57b714e..29162a04 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -395,6 +395,23 @@ void ScaLBL_ColorModel::Run(){ int nprocs=nprocx*nprocy*nprocz; const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + int lead_timesteps = 0; + double tolerance = 1.f; + + double capillary_number; + bool SET_CAPILLARY_NUMBER=false; + if (color_db->keyExists( "capillary_number" )){ + capillary_number = color_db->getScalar( "capillary_number" ); + SET_CAPILLARY_NUMBER=true; + } + else{ + capillary_number=0; + } + if (BoundaryCondition != 0 && SET_CAPILLARY_NUMBER==true){ + if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 \n"); + SET_CAPILLARY_NUMBER=false; + } + int morph_interval; double morph_delta; if (analysis_db->keyExists( "morph_delta" )){ @@ -409,7 +426,8 @@ void ScaLBL_ColorModel::Run(){ else{ morph_interval=100000; } - + int analysis_interval = analysis_db->getScalar( "analysis_interval" ); + if (rank==0){ printf("********************************************************\n"); @@ -511,9 +529,28 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); - if (timestep%morph_interval == 0){ + if (timestep%morph_interval == 0 || tolerance < 0.01){ MorphInit(beta,morph_delta); MPI_Barrier(comm); + lead_timesteps = 0; + tolerance = 1.f; + } + if (SET_CAPILLARY_NUMBER && timestep%analysis_interval == analysis_interval-20){ + lead_timesteps += analysis_interval; + double vA_z = Averages->van_global(2); + double vB_z = Averages->vaw_global(2); + double volB = Averages->wet_morph->V(); + double volA = Averages->nonwet_morph->V(); + double muA = rhoA*(tauA-0.5)/3.f; + double muB = rhoB*(tauB-0.5)/3.f; + double Ca = fabs(volA*muA*vA_z + volB*muB*vB_z)/(5.796*alpha*double(Nx*Ny*Nz)); + if (rank == 0) printf(" Measured capillary number %f \n ",Ca); + if (lead_timesteps > 5000){ + Fz *= capillary_number / Ca; + if (Fz > 1e-3) Fz = 1e-3; // impose ceiling for stability + tolerance = fabs(capillary_number - Ca) / capillary_number ; + if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); + } } } From 6d5732d66f58eeddad2cc55464b82c29977eee5c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 24 Oct 2018 12:58:47 -0400 Subject: [PATCH 17/79] added ramp up for steady state --- models/ColorModel.cpp | 46 ++++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 22 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 29162a04..12a742e5 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -529,30 +529,32 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); - if (timestep%morph_interval == 0 || tolerance < 0.01){ - MorphInit(beta,morph_delta); - MPI_Barrier(comm); - lead_timesteps = 0; - tolerance = 1.f; - } - if (SET_CAPILLARY_NUMBER && timestep%analysis_interval == analysis_interval-20){ - lead_timesteps += analysis_interval; - double vA_z = Averages->van_global(2); - double vB_z = Averages->vaw_global(2); - double volB = Averages->wet_morph->V(); - double volA = Averages->nonwet_morph->V(); - double muA = rhoA*(tauA-0.5)/3.f; - double muB = rhoB*(tauB-0.5)/3.f; - double Ca = fabs(volA*muA*vA_z + volB*muB*vB_z)/(5.796*alpha*double(Nx*Ny*Nz)); - if (rank == 0) printf(" Measured capillary number %f \n ",Ca); - if (lead_timesteps > 5000){ - Fz *= capillary_number / Ca; - if (Fz > 1e-3) Fz = 1e-3; // impose ceiling for stability - tolerance = fabs(capillary_number - Ca) / capillary_number ; - if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); + if (timestep > 50000){ + // allow initial ramp-up to get closer to steady state + if (timestep%morph_interval == 0 || tolerance < 0.01){ + MorphInit(beta,morph_delta); + MPI_Barrier(comm); + lead_timesteps = 0; + tolerance = 1.f; + } + if (SET_CAPILLARY_NUMBER && timestep%analysis_interval == analysis_interval-20 ){ + lead_timesteps += analysis_interval; + double vA_z = Averages->van_global(2); + double vB_z = Averages->vaw_global(2); + double volB = Averages->wet_morph->V(); + double volA = Averages->nonwet_morph->V(); + double muA = rhoA*(tauA-0.5)/3.f; + double muB = rhoB*(tauB-0.5)/3.f; + double Ca = fabs(volA*muA*vA_z + volB*muB*vB_z)/(5.796*alpha*double(Nx*Ny*Nz)); + if (rank == 0) printf(" Measured capillary number %f \n ",Ca); + if (lead_timesteps > 5000){ + Fz *= capillary_number / Ca; + if (Fz > 1e-3) Fz = 1e-3; // impose ceiling for stability + tolerance = fabs(capillary_number - Ca) / capillary_number ; + if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); + } } } - } analysis.finish(); PROFILE_STOP("Loop"); From 23dd734a40e16b2f1b84f37684a4f70a3402cac7 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 24 Oct 2018 13:26:37 -0400 Subject: [PATCH 18/79] LBM parameters in TwoPhase --- analysis/TwoPhase.cpp | 21 +++++++++++++++++---- analysis/TwoPhase.h | 6 ++++++ models/ColorModel.cpp | 3 +++ 3 files changed, 26 insertions(+), 4 deletions(-) diff --git a/analysis/TwoPhase.cpp b/analysis/TwoPhase.cpp index 7c4ac45d..aeb7b827 100644 --- a/analysis/TwoPhase.cpp +++ b/analysis/TwoPhase.cpp @@ -144,7 +144,7 @@ TwoPhase::TwoPhase(std::shared_ptr dm): { // If timelog is empty, write a short header to list the averages //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"time dEs "); // Timestep, Change in Surface Energy + fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn "); // Timestep, Change in Surface Energy fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz "); @@ -175,7 +175,7 @@ TwoPhase::TwoPhase(std::shared_ptr dm): sprintf(LocalRankFilename,"%s%s","timelog.tcat.",LocalRankString); TIMELOG = fopen(LocalRankFilename,"a+"); //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"time "); // Timestep, Change in Surface Energy + fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn ");; // Timestep, fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz "); @@ -308,6 +308,19 @@ void TwoPhase::Initialize() wwndnw = 0.0; wwnsdnwn = 0.0; Jwnwwndnw=0.0; } +void TwoPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha) +{ + Fx = force_x; + Fy = force_y; + Fz = force_z; + rho_n = rhoA; + rho_w = rhoB; + nu_n = (tauA-0.5)/3.f; + nu_w = (tauB-0.5)/3.f; + gamma_wn = 5.796*alpha; + +} + /* void TwoPhase::SetupCubes(Domain &Dm) { @@ -1162,7 +1175,7 @@ void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) void TwoPhase::PrintAll(int timestep) { if (Dm->rank()==0){ - fprintf(TIMELOG,"%i %.5g ",timestep,dEs); // change in surface energy + fprintf(TIMELOG,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn); fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas fprintf(TIMELOG,"%.5g %.5g ",Jwn_global, Kwn_global); // curvature of wn interface @@ -1189,7 +1202,7 @@ void TwoPhase::PrintAll(int timestep) else{ sat_w = 1.0 - nwp_volume/(nwp_volume+wp_volume); - fprintf(TIMELOG,"%i ",timestep); // change in surface energy + fprintf(TIMELOG,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn); fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw,pan); // saturation and pressure fprintf(TIMELOG,"%.5g %.5g %.5g ",awn,ans,aws); // interfacial areas fprintf(TIMELOG,"%.5g %.5g ",Jwn, Kwn); // curvature of wn interface diff --git a/analysis/TwoPhase.h b/analysis/TwoPhase.h index 74794500..5211f2d9 100644 --- a/analysis/TwoPhase.h +++ b/analysis/TwoPhase.h @@ -88,6 +88,11 @@ public: double efawns,efawns_global; // averaged contact angle double euler,Kn,Jn,An; double euler_global,Kn_global,Jn_global,An_global; + + double rho_n, rho_w; + double nu_n, nu_w; + double gamma_wn; + double Fx, Fy, Fz; double Jwn,Jwn_global; // average mean curavture - wn interface double Kwn,Kwn_global; // average Gaussian curavture - wn interface @@ -169,6 +174,7 @@ public: int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel); void SortBlobs(); void PrintComponents(int timestep); + void SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha); }; diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 12a742e5..480fbb39 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -151,6 +151,8 @@ void ScaLBL_ColorModel::ReadInput(){ CalcDist(Averages->SDs,id_solid,*Mask); if (rank == 0) cout << "Domain set." << endl; + + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); } void ScaLBL_ColorModel::AssignComponentLabels(double *phase) @@ -553,6 +555,7 @@ void ScaLBL_ColorModel::Run(){ tolerance = fabs(capillary_number - Ca) / capillary_number ; if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); } + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); } } } From 757805628cdc8a349d8e9718f5c750504ef02e1a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 27 Oct 2018 18:12:47 -0400 Subject: [PATCH 19/79] added forcing floor for morphLBM --- models/ColorModel.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 480fbb39..9f2fb637 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -552,6 +552,7 @@ void ScaLBL_ColorModel::Run(){ if (lead_timesteps > 5000){ Fz *= capillary_number / Ca; if (Fz > 1e-3) Fz = 1e-3; // impose ceiling for stability + if (Fz < 1e-6) Fz = 1e-6; // impose floor so we don't do something super dumb tolerance = fabs(capillary_number - Ca) / capillary_number ; if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); } From a016ae6b4257d5c2db5caa4e726b1c762609a3a2 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 28 Oct 2018 10:17:47 -0400 Subject: [PATCH 20/79] force tuning for auto relperm --- models/ColorModel.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 9f2fb637..c4011ddd 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -552,7 +552,9 @@ void ScaLBL_ColorModel::Run(){ if (lead_timesteps > 5000){ Fz *= capillary_number / Ca; if (Fz > 1e-3) Fz = 1e-3; // impose ceiling for stability - if (Fz < 1e-6) Fz = 1e-6; // impose floor so we don't do something super dumb + if (Fz < 1e-7) Fz = 1e-6; // impose floor so we don't do something super dumb + if (Fz*vA_z < 0.f || Fz*vB_z < 0.f) + Fz *= 2.f; // bigger forces are needed if flow rate is "too small" tolerance = fabs(capillary_number - Ca) / capillary_number ; if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); } From 52d60f0e01bb8633dc4701ecf7bdd3757a60d1a5 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 28 Oct 2018 10:18:06 -0400 Subject: [PATCH 21/79] force tuning for auto relperm --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index c4011ddd..51fac042 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -552,7 +552,7 @@ void ScaLBL_ColorModel::Run(){ if (lead_timesteps > 5000){ Fz *= capillary_number / Ca; if (Fz > 1e-3) Fz = 1e-3; // impose ceiling for stability - if (Fz < 1e-7) Fz = 1e-6; // impose floor so we don't do something super dumb + if (Fz < 1e-7) Fz = 1e-7; // impose floor so we don't do something super dumb if (Fz*vA_z < 0.f || Fz*vB_z < 0.f) Fz *= 2.f; // bigger forces are needed if flow rate is "too small" tolerance = fabs(capillary_number - Ca) / capillary_number ; From 66299f76cb8c959b98697bc13f4142d37e74c676 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 1 Nov 2018 15:18:14 -0400 Subject: [PATCH 22/79] Force can change in any direction for morph --- models/ColorModel.cpp | 39 +++++++++++++++++++++++++++++++++------ 1 file changed, 33 insertions(+), 6 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 51fac042..7e6aa0e3 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -541,20 +541,47 @@ void ScaLBL_ColorModel::Run(){ } if (SET_CAPILLARY_NUMBER && timestep%analysis_interval == analysis_interval-20 ){ lead_timesteps += analysis_interval; + double vA_x = Averages->van_global(0); + double vA_y = Averages->van_global(1); double vA_z = Averages->van_global(2); - double vB_z = Averages->vaw_global(2); + + double vB_x = Averages->vaw_global(0); + double vB_y = Averages->vaw_global(1); + double vB_z = Averages->vaw_global(2); + double volB = Averages->wet_morph->V(); double volA = Averages->nonwet_morph->V(); double muA = rhoA*(tauA-0.5)/3.f; double muB = rhoB*(tauB-0.5)/3.f; - double Ca = fabs(volA*muA*vA_z + volB*muB*vB_z)/(5.796*alpha*double(Nx*Ny*Nz)); + + double flow_rate_A = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); + double flow_rate_B = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); + + double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz)); + if (rank == 0) printf(" Measured capillary number %f \n ",Ca); if (lead_timesteps > 5000){ + Fx *= capillary_number / Ca; + Fy *= capillary_number / Ca; Fz *= capillary_number / Ca; - if (Fz > 1e-3) Fz = 1e-3; // impose ceiling for stability - if (Fz < 1e-7) Fz = 1e-7; // impose floor so we don't do something super dumb - if (Fz*vA_z < 0.f || Fz*vB_z < 0.f) - Fz *= 2.f; // bigger forces are needed if flow rate is "too small" + + double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); + + if (force_magnitude > 1e-3){ + Fx *= 1e-3/force_magnitude; // impose ceiling for stability + Fy *= 1e-3/force_magnitude; + Fz *= 1e-3/force_magnitude; + } + if (force_magnitude < 1e-7){ + Fx *= 1e-7/force_magnitude; // impose floor + Fy *= 1e-7/force_magnitude; + Fz *= 1e-7/force_magnitude; + } + if (Fx*vA_x + Fy*vA_y + Fz*vA_z < 0.f || Fx*vB_x + Fy*vB_y +Fz*vB_z < 0.f){ + Fx *= 2.f; // bigger forces are needed if flow rate is "too small" + Fy *= 2.f; + Fz *= 2.f; + } tolerance = fabs(capillary_number - Ca) / capillary_number ; if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); } From 4e2bc18d8a516fd0aebe42d7c302aa61f97d9517 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 07:38:09 -0400 Subject: [PATCH 23/79] compute signdist in morphdrain_pp --- tests/lbpm_morphdrain_pp.cpp | 36 +++++++++++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 45378811..c437b82d 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -121,16 +121,46 @@ int main(int argc, char **argv) ydim=Dm.Ny-2; zdim=Dm.Nz-2; fillHalo fillData(Dm.Comm, Dm.rank_info,{xdim,ydim,zdim},{1,1,1},0,1); + + // Generate the signed distance map + // Initialize the domain and communication + Array id_solid(Nx,Ny,Nz); + DoubleArray SignDist(Nx,Ny,Nz); - DoubleArray SignDist(nx,ny,nz); - // Read the signed distance from file + int count = 0; + // Solve for the position of the solid phase + for (int k=0;kid[n] > 0) id_solid(i,j,k) = 1; + else id_solid(i,j,k) = 0; + } + } + } + // Initialize the signed distance function + for (int k=0;k Date: Fri, 2 Nov 2018 07:46:15 -0400 Subject: [PATCH 24/79] compute signdist in morphdrain_pp --- tests/lbpm_morphdrain_pp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index c437b82d..905172d7 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -11,7 +11,7 @@ #include #include "common/Array.h" #include "common/Domain.h" - +#include "analysis/distance.h" //************************************************************************* // Morpohological drainage pre-processor From 57fe1ee0ba6680584b8f3f06b42e415bde438702 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 07:50:46 -0400 Subject: [PATCH 25/79] compute signdist in morphdrain_pp --- tests/lbpm_morphdrain_pp.cpp | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 905172d7..5d2da8dd 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -124,26 +124,25 @@ int main(int argc, char **argv) // Generate the signed distance map // Initialize the domain and communication - Array id_solid(Nx,Ny,Nz); - DoubleArray SignDist(Nx,Ny,Nz); + Array id_solid(nx,ny,nz); + DoubleArray SignDist(nx,ny,nz); - int count = 0; // Solve for the position of the solid phase - for (int k=0;kid[n] > 0) id_solid(i,j,k) = 1; + if (Dm.id[n] > 0) id_solid(i,j,k) = 1; else id_solid(i,j,k) = 0; } } } // Initialize the signed distance function - for (int k=0;k Date: Fri, 2 Nov 2018 07:59:53 -0400 Subject: [PATCH 26/79] compute signdist in morphdrain_pp --- tests/lbpm_morphdrain_pp.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 5d2da8dd..af481bc4 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -96,8 +96,8 @@ int main(int argc, char **argv) char LocalRankFilename[40]; int BoundaryCondition=1; - Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition); - + Domain(domain_db,comm); + nx+=2; ny+=2; nz+=2; int N = nx*ny*nz; char *id; From e6be6da1ed468b3500c05928b26c5519cdb0fc4a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 08:05:10 -0400 Subject: [PATCH 27/79] compute signdist in morphdrain_pp --- tests/lbpm_morphdrain_pp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index af481bc4..4d96f1db 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -96,7 +96,7 @@ int main(int argc, char **argv) char LocalRankFilename[40]; int BoundaryCondition=1; - Domain(domain_db,comm); + Domain Dm(domain_db,comm); nx+=2; ny+=2; nz+=2; int N = nx*ny*nz; From 19bf87a33553e187bd13ad0be5a3962b2804240e Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 08:48:55 -0400 Subject: [PATCH 28/79] compute signdist in morphdrain_pp --- tests/lbpm_morphdrain_pp.cpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index 4d96f1db..f92c5143 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -116,6 +116,14 @@ int main(int argc, char **argv) Dm.CommInit(); + sprintf(LocalRankFilename,"ID.%05i",rank); + size_t readID; + FILE *IDFILE = fopen(LocalRankFilename,"rb"); + if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); + readID=fread(id,1,N,IDFILE); + if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank); + fclose(IDFILE); + int xdim,ydim,zdim; xdim=Dm.Nx-2; ydim=Dm.Ny-2; @@ -133,8 +141,8 @@ int main(int argc, char **argv) for (int i=0;i 0) id_solid(i,j,k) = 1; - else id_solid(i,j,k) = 0; + if (id[n] > 0) id_solid(i,j,k) = 1; + else id_solid(i,j,k) = 0; } } } @@ -161,15 +169,6 @@ int main(int argc, char **argv) fclose(DIST); */ fillData.fill(SignDist); - - sprintf(LocalRankFilename,"ID.%05i",rank); - size_t readID; - FILE *IDFILE = fopen(LocalRankFilename,"rb"); - if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); - readID=fread(id,1,N,IDFILE); - if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank); - fclose(IDFILE); - Dm.CommInit(); int iproc = Dm.iproc(); From b922bab34e12118d9043379b77169462d924cd42 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 11:33:12 -0400 Subject: [PATCH 29/79] update morph LBM to use target saturation list --- models/ColorModel.cpp | 51 ++++++++++++++++++++++++++++++++++++------- models/ColorModel.h | 2 +- 2 files changed, 44 insertions(+), 9 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 7e6aa0e3..c1588dad 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -80,8 +80,6 @@ void ScaLBL_ColorModel::ReadParams(string filename){ outletA=0.f; outletB=1.f; - if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details) - // Read domain parameters auto L = domain_db->getVector( "L" ); auto size = domain_db->getVector( "n" ); @@ -96,6 +94,8 @@ void ScaLBL_ColorModel::ReadParams(string filename){ nprocx = nproc[0]; nprocy = nproc[1]; nprocz = nproc[2]; + + if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details) } void ScaLBL_ColorModel::SetDomain(){ @@ -397,11 +397,20 @@ void ScaLBL_ColorModel::Run(){ int nprocs=nprocx*nprocy*nprocz; const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); + bool SET_CAPILLARY_NUMBER = false; + bool MORPH_ADAPT = false; int lead_timesteps = 0; + int ramp_timesteps = 50000; double tolerance = 1.f; + int target_saturation_index=0; + std::vector target_saturation; + double TARGET_SATURATION = 0.f; + if (domain_db->keyExists( "Sw" )){ + target_saturation = domain_db->getVector( "Sw" ); + } + double capillary_number; - bool SET_CAPILLARY_NUMBER=false; if (color_db->keyExists( "capillary_number" )){ capillary_number = color_db->getScalar( "capillary_number" ); SET_CAPILLARY_NUMBER=true; @@ -531,10 +540,35 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); - if (timestep > 50000){ + if (timestep > ramp_timesteps){ // allow initial ramp-up to get closer to steady state if (timestep%morph_interval == 0 || tolerance < 0.01){ - MorphInit(beta,morph_delta); + tolerance = 1.f; + MORPH_ADAPT = true; + TARGET_SATURATION = target_saturation[target_saturation_index++]; + double volB = Averages->wet_morph->V(); + double volA = Averages->nonwet_morph->V(); + double current_saturation = volB/(volA+volB); + if (morph_delta > 0.f){ + // wetting phase saturation will decrease + while (current_saturation > TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ + TARGET_SATURATION = target_saturation[target_saturation_index++]; + } + } + else{ + // wetting phase saturation will increase + while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ + TARGET_SATURATION = target_saturation[target_saturation_index++]; + } + } + } + if (MORPH_ADAPT){ + double volB = Averages->wet_morph->V(); + double volA = Averages->nonwet_morph->V(); + double delta_volume = MorphInit(beta,morph_delta); + if ((1.f+delta_volume)*volA/(volA + volB) > 1.f - TARGET_SATURATION){ + MORPH_ADAPT = false; + } MPI_Barrier(comm); lead_timesteps = 0; tolerance = 1.f; @@ -612,7 +646,7 @@ void ScaLBL_ColorModel::Run(){ // ************************************************************************ } -void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ +double ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); double vF = 0.f; @@ -714,7 +748,8 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm); volume_final=count_global; - if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", (volume_final-volume_initial)/volume_initial); + double delta_volume = (volume_final-volume_initial)/volume_initial; + if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume); // 6. copy back to the device //if (rank==0) printf("MorphInit: copy data back to device\n"); @@ -735,7 +770,7 @@ void ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){ ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); } } - + return delta_volume; } void ScaLBL_ColorModel::WriteDebug(){ diff --git a/models/ColorModel.h b/models/ColorModel.h index 5b5e9ccf..67675207 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -77,7 +77,7 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr db0); void AssignComponentLabels(double *phase); - void MorphInit(const double beta, const double morph_delta); + double MorphInit(const double beta, const double morph_delta); }; From e8c1db93dcbdaa262cb135b3fd8339bcc0b8dc4a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 12:47:42 -0400 Subject: [PATCH 30/79] read Sw from color db --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index c1588dad..42f1eb4c 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -406,8 +406,8 @@ void ScaLBL_ColorModel::Run(){ int target_saturation_index=0; std::vector target_saturation; double TARGET_SATURATION = 0.f; - if (domain_db->keyExists( "Sw" )){ - target_saturation = domain_db->getVector( "Sw" ); + if (color_db->keyExists( "Sw" )){ + target_saturation = color_db->getVector( "Sw" ); } double capillary_number; From 33e3d4395ebaff592e058dd3781bb6949778dade Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 12:50:10 -0400 Subject: [PATCH 31/79] updating morph tools compute dist inside --- tests/lbpm_morphdrain_pp.cpp | 2 +- tests/lbpm_morphopen_pp.cpp | 43 ++++++++++++++++++++++++++++++------ 2 files changed, 37 insertions(+), 8 deletions(-) diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index f92c5143..70798e88 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -121,7 +121,7 @@ int main(int argc, char **argv) FILE *IDFILE = fopen(LocalRankFilename,"rb"); if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); readID=fread(id,1,N,IDFILE); - if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank); + if (readID != size_t(N)) printf("lbpm_morphdrain_pp: Error reading ID (rank=%i) \n",rank); fclose(IDFILE); int xdim,ydim,zdim; diff --git a/tests/lbpm_morphopen_pp.cpp b/tests/lbpm_morphopen_pp.cpp index a657b4f1..07608306 100644 --- a/tests/lbpm_morphopen_pp.cpp +++ b/tests/lbpm_morphopen_pp.cpp @@ -106,16 +106,45 @@ int main(int argc, char **argv) char *id; id = new char [N]; + sprintf(LocalRankFilename,"ID.%05i",rank); + size_t readID; + FILE *IDFILE = fopen(LocalRankFilename,"rb"); + if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx"); + readID=fread(id,1,N,IDFILE); + if (readID != size_t(N)) printf("lbpm_morphopen_pp: Error reading ID (rank=%i) \n",rank); + fclose(IDFILE); + nx+=2; ny+=2; nz+=2; + // Generate the signed distance map + // Initialize the domain and communication + Array id_solid(nx,ny,nz); DoubleArray SignDist(nx,ny,nz); - // Read the signed distance from file - sprintf(LocalRankFilename,"SignDist.%05i",rank); - FILE *DIST = fopen(LocalRankFilename,"rb"); - size_t ReadSignDist; - ReadSignDist=fread(SignDist.data(),8,N,DIST); - if (ReadSignDist != size_t(N)) printf("lbpm_morphopen_pp: Error reading signed distance function (rank=%i)\n",rank); - fclose(DIST); + + // Solve for the position of the solid phase + for (int k=0;k 0) id_solid(i,j,k) = 1; + else id_solid(i,j,k) = 0; + } + } + } + // Initialize the signed distance function + for (int k=0;k Date: Fri, 2 Nov 2018 12:53:09 -0400 Subject: [PATCH 32/79] updating morph tools compute dist inside --- tests/lbpm_morphopen_pp.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/lbpm_morphopen_pp.cpp b/tests/lbpm_morphopen_pp.cpp index 07608306..1c9ce2ea 100644 --- a/tests/lbpm_morphopen_pp.cpp +++ b/tests/lbpm_morphopen_pp.cpp @@ -11,6 +11,7 @@ #include #include "common/Array.h" #include "common/Domain.h" +#include "analysis/distance.h" //************************************************************************* From 2e7d086b31ce6d6f944901a174930a2449e9adab Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 12:57:03 -0400 Subject: [PATCH 33/79] updating morph tools compute dist inside --- tests/lbpm_morphopen_pp.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/lbpm_morphopen_pp.cpp b/tests/lbpm_morphopen_pp.cpp index 1c9ce2ea..3dbc9c74 100644 --- a/tests/lbpm_morphopen_pp.cpp +++ b/tests/lbpm_morphopen_pp.cpp @@ -145,7 +145,7 @@ int main(int argc, char **argv) } if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); - CalcDist(SignDist,id_solid,Dm); + CalcDist(SignDist,id_solid,*Dm); MPI_Barrier(comm); double count,countGlobal,totalGlobal; @@ -422,11 +422,7 @@ int main(int argc, char **argv) if (rank==0) printf("Writing ID file \n"); sprintf(LocalRankFilename,"ID.%05i",rank); - size_t readID; - FILE *ID = fopen(LocalRankFilename,"rb"); - readID=fread(Dm->id,1,N,ID); - if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID \n"); - fclose(ID); + // Preserve mineral labels for (int k=0; k Date: Fri, 2 Nov 2018 12:59:01 -0400 Subject: [PATCH 34/79] added target_saturation key to color --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 42f1eb4c..66b6d1fe 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -406,8 +406,8 @@ void ScaLBL_ColorModel::Run(){ int target_saturation_index=0; std::vector target_saturation; double TARGET_SATURATION = 0.f; - if (color_db->keyExists( "Sw" )){ - target_saturation = color_db->getVector( "Sw" ); + if (color_db->keyExists( "target_saturation" )){ + target_saturation = color_db->getVector( "target_saturation" ); } double capillary_number; From ca78ab25116cfa3543244354cc4ec7808778874d Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 13:00:26 -0400 Subject: [PATCH 35/79] updating morph tools compute dist inside --- tests/lbpm_morphopen_pp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/lbpm_morphopen_pp.cpp b/tests/lbpm_morphopen_pp.cpp index 3dbc9c74..bc591d70 100644 --- a/tests/lbpm_morphopen_pp.cpp +++ b/tests/lbpm_morphopen_pp.cpp @@ -432,7 +432,7 @@ int main(int argc, char **argv) } } } - ID = fopen(LocalRankFilename,"wb"); + FILE *ID = fopen(LocalRankFilename,"wb"); fwrite(id,1,N,ID); fclose(ID); } From 9f7b530587d030adcef40d9973b3d9eb472e0872 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 13:51:50 -0400 Subject: [PATCH 36/79] tweak to morph --- models/ColorModel.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 66b6d1fe..73feec3a 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -559,10 +559,12 @@ void ScaLBL_ColorModel::Run(){ // wetting phase saturation will increase while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ TARGET_SATURATION = target_saturation[target_saturation_index++]; + if (rank==0) printf(" Set target saturation as %f (currently %f)\n",TARGET_SATURATION,current_saturation); } } } if (MORPH_ADAPT){ + if (rank==0) printf("***Morphological step***\n"); double volB = Averages->wet_morph->V(); double volA = Averages->nonwet_morph->V(); double delta_volume = MorphInit(beta,morph_delta); @@ -591,7 +593,7 @@ void ScaLBL_ColorModel::Run(){ double flow_rate_A = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); double flow_rate_B = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); - double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz)); + double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); if (rank == 0) printf(" Measured capillary number %f \n ",Ca); if (lead_timesteps > 5000){ From f9d5c2965af4bc1a137b4f5357de4b3586aa2e45 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 14:17:17 -0400 Subject: [PATCH 37/79] tweak to morph --- models/ColorModel.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 73feec3a..6e3ac910 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -568,7 +568,9 @@ void ScaLBL_ColorModel::Run(){ double volB = Averages->wet_morph->V(); double volA = Averages->nonwet_morph->V(); double delta_volume = MorphInit(beta,morph_delta); - if ((1.f+delta_volume)*volA/(volA + volB) > 1.f - TARGET_SATURATION){ + volA += delta_volume; + volB -= delta_volume; + if (volB/(volA + volB) > TARGET_SATURATION){ MORPH_ADAPT = false; } MPI_Barrier(comm); @@ -750,8 +752,8 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta) MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm); volume_final=count_global; - double delta_volume = (volume_final-volume_initial)/volume_initial; - if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume); + double delta_volume = (volume_final-volume_initial); + if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial); // 6. copy back to the device //if (rank==0) printf("MorphInit: copy data back to device\n"); From d7ada6f10c4140e73e53717c8c8ce2b0b5abac07 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 14:33:48 -0400 Subject: [PATCH 38/79] fixing morphopen --- tests/lbpm_morphopen_pp.cpp | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/tests/lbpm_morphopen_pp.cpp b/tests/lbpm_morphopen_pp.cpp index bc591d70..0d280fab 100644 --- a/tests/lbpm_morphopen_pp.cpp +++ b/tests/lbpm_morphopen_pp.cpp @@ -165,7 +165,9 @@ int main(int argc, char **argv) for (int j=0; jid[n]; - } - } - } FILE *ID = fopen(LocalRankFilename,"wb"); fwrite(id,1,N,ID); fclose(ID); From c062809e01c5f8e5d321a61648a0bbf1c7cdc66f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 15:31:11 -0400 Subject: [PATCH 39/79] fix morph too often --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 6e3ac910..ef252ca8 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -542,7 +542,7 @@ void ScaLBL_ColorModel::Run(){ if (timestep > ramp_timesteps){ // allow initial ramp-up to get closer to steady state - if (timestep%morph_interval == 0 || tolerance < 0.01){ + if (timestep%morph_interval-20 == 0 || tolerance < 0.01){ tolerance = 1.f; MORPH_ADAPT = true; TARGET_SATURATION = target_saturation[target_saturation_index++]; @@ -563,7 +563,7 @@ void ScaLBL_ColorModel::Run(){ } } } - if (MORPH_ADAPT){ + if (MORPH_ADAPT && timestep%analysis_interval == analysis_interval-20 ){ if (rank==0) printf("***Morphological step***\n"); double volB = Averages->wet_morph->V(); double volA = Averages->nonwet_morph->V(); From 76797aa005900ce261dc9144049f44acdcef9cda Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 16:01:09 -0400 Subject: [PATCH 40/79] fixed morph volume calc --- analysis/TwoPhase.h | 6 ++++++ models/ColorModel.cpp | 12 ++++++------ 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/analysis/TwoPhase.h b/analysis/TwoPhase.h index 5211f2d9..01df349d 100644 --- a/analysis/TwoPhase.h +++ b/analysis/TwoPhase.h @@ -175,6 +175,12 @@ public: void SortBlobs(); void PrintComponents(int timestep); void SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha); + double Volume_w(){ + return wp_volume_global; + } + double Volume_n(){ + return nwp_volume_global; + } }; diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index ef252ca8..c2e52b7a 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -546,8 +546,8 @@ void ScaLBL_ColorModel::Run(){ tolerance = 1.f; MORPH_ADAPT = true; TARGET_SATURATION = target_saturation[target_saturation_index++]; - double volB = Averages->wet_morph->V(); - double volA = Averages->nonwet_morph->V(); + double volB = Averages->Volume_w(); + double volA = Averages->Volume_n(); double current_saturation = volB/(volA+volB); if (morph_delta > 0.f){ // wetting phase saturation will decrease @@ -565,8 +565,8 @@ void ScaLBL_ColorModel::Run(){ } if (MORPH_ADAPT && timestep%analysis_interval == analysis_interval-20 ){ if (rank==0) printf("***Morphological step***\n"); - double volB = Averages->wet_morph->V(); - double volA = Averages->nonwet_morph->V(); + double volB = Averages->Volume_w(); + double volA = Averages->Volume_n(); double delta_volume = MorphInit(beta,morph_delta); volA += delta_volume; volB -= delta_volume; @@ -587,8 +587,8 @@ void ScaLBL_ColorModel::Run(){ double vB_y = Averages->vaw_global(1); double vB_z = Averages->vaw_global(2); - double volB = Averages->wet_morph->V(); - double volA = Averages->nonwet_morph->V(); + double volB = Averages->Volume_w(); + double volA = Averages->Volume_n(); double muA = rhoA*(tauA-0.5)/3.f; double muB = rhoB*(tauB-0.5)/3.f; From 6d3e7ebc3dc21306a9a29b43936e1d3a8107cc59 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 16:07:39 -0400 Subject: [PATCH 41/79] fixed convention for morphological analysis --- analysis/TwoPhase.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/analysis/TwoPhase.cpp b/analysis/TwoPhase.cpp index aeb7b827..14033f98 100644 --- a/analysis/TwoPhase.cpp +++ b/analysis/TwoPhase.cpp @@ -572,15 +572,15 @@ void TwoPhase::ComputeLocal() n = k*Nx*Ny+j*Nx+i; if (!(Dm->id[n] > 0)){ // Solid phase - phase_label(i,j,k) = 0; + phase_label(i,j,k) = 1; } else if (SDn(i,j,k) < 0.0){ // wetting phase - phase_label(i,j,k) = 1; + phase_label(i,j,k) = 0; } else { // non-wetting phase - phase_label(i,j,k) = 0; + phase_label(i,j,k) = 1; } phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0; } @@ -596,15 +596,15 @@ void TwoPhase::ComputeLocal() n = k*Nx*Ny+j*Nx+i; if (!(Dm->id[n] > 0)){ // Solid phase - phase_label(i,j,k) = 0; + phase_label(i,j,k) = 1; } else if (SDn(i,j,k) < 0.0){ // wetting phase - phase_label(i,j,k) = 0; + phase_label(i,j,k) = 1; } else { // non-wetting phase - phase_label(i,j,k) = 1; + phase_label(i,j,k) = 0; } phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0; } From 6f29a1461af6fa75d567915d1d2b60ddabc931e8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 16:25:55 -0400 Subject: [PATCH 42/79] don't run amok with MorphAdapt --- models/ColorModel.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index c2e52b7a..8a819b90 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -570,6 +570,7 @@ void ScaLBL_ColorModel::Run(){ double delta_volume = MorphInit(beta,morph_delta); volA += delta_volume; volB -= delta_volume; + MORPH_ADAPT = false; if (volB/(volA + volB) > TARGET_SATURATION){ MORPH_ADAPT = false; } From 53004a10e03cf5254bdf10b00a7ae5043e0b71a5 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 16:30:44 -0400 Subject: [PATCH 43/79] fix morph adapt --- models/ColorModel.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8a819b90..3914964d 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -551,13 +551,13 @@ void ScaLBL_ColorModel::Run(){ double current_saturation = volB/(volA+volB); if (morph_delta > 0.f){ // wetting phase saturation will decrease - while (current_saturation > TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ + while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ TARGET_SATURATION = target_saturation[target_saturation_index++]; } } else{ // wetting phase saturation will increase - while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ + while (current_saturation > TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ TARGET_SATURATION = target_saturation[target_saturation_index++]; if (rank==0) printf(" Set target saturation as %f (currently %f)\n",TARGET_SATURATION,current_saturation); } @@ -570,7 +570,7 @@ void ScaLBL_ColorModel::Run(){ double delta_volume = MorphInit(beta,morph_delta); volA += delta_volume; volB -= delta_volume; - MORPH_ADAPT = false; + //MORPH_ADAPT = false; if (volB/(volA + volB) > TARGET_SATURATION){ MORPH_ADAPT = false; } From 30df673fc0503b281298aa511dbd3206f2f4ff6a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 16:57:12 -0400 Subject: [PATCH 44/79] tweak morph conditions --- models/ColorModel.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 3914964d..3bfcd798 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -408,6 +408,7 @@ void ScaLBL_ColorModel::Run(){ double TARGET_SATURATION = 0.f; if (color_db->keyExists( "target_saturation" )){ target_saturation = color_db->getVector( "target_saturation" ); + TARGET_SATURATION = target_saturation[0]; } double capillary_number; From 3e5b7c6e53386155d49a90393911fada2cada5c7 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 2 Nov 2018 17:48:47 -0400 Subject: [PATCH 45/79] debugging --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 3bfcd798..1ae67e84 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -565,7 +565,7 @@ void ScaLBL_ColorModel::Run(){ } } if (MORPH_ADAPT && timestep%analysis_interval == analysis_interval-20 ){ - if (rank==0) printf("***Morphological step***\n"); + if (rank==0) printf("***Morphological step with target saturation %f ***\n",TARGET_SATURATION); double volB = Averages->Volume_w(); double volA = Averages->Volume_n(); double delta_volume = MorphInit(beta,morph_delta); From 398fe3b1869befe700df9dfa06cfcf545ed191c1 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 3 Nov 2018 06:57:38 -0400 Subject: [PATCH 46/79] debugging --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 1ae67e84..84bd4430 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -543,7 +543,7 @@ void ScaLBL_ColorModel::Run(){ if (timestep > ramp_timesteps){ // allow initial ramp-up to get closer to steady state - if (timestep%morph_interval-20 == 0 || tolerance < 0.01){ + if (!MORPH_ADAPT && (timestep%morph_interval-20 == 0 || tolerance < 0.01)){ tolerance = 1.f; MORPH_ADAPT = true; TARGET_SATURATION = target_saturation[target_saturation_index++]; From 0e22b434659b0d6b2a83007a8e8f92f912b6d1c8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 4 Nov 2018 10:51:45 -0500 Subject: [PATCH 47/79] separate morph timesteps to avoid problems --- models/ColorModel.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 84bd4430..5d60d7bd 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -399,7 +399,7 @@ void ScaLBL_ColorModel::Run(){ bool SET_CAPILLARY_NUMBER = false; bool MORPH_ADAPT = false; - int lead_timesteps = 0; + int morph_timesteps = 0; int ramp_timesteps = 50000; double tolerance = 1.f; @@ -543,7 +543,7 @@ void ScaLBL_ColorModel::Run(){ if (timestep > ramp_timesteps){ // allow initial ramp-up to get closer to steady state - if (!MORPH_ADAPT && (timestep%morph_interval-20 == 0 || tolerance < 0.01)){ + if (!MORPH_ADAPT && (morph_timesteps%morph_interval == 0 || tolerance < 0.01)){ tolerance = 1.f; MORPH_ADAPT = true; TARGET_SATURATION = target_saturation[target_saturation_index++]; @@ -576,11 +576,11 @@ void ScaLBL_ColorModel::Run(){ MORPH_ADAPT = false; } MPI_Barrier(comm); - lead_timesteps = 0; + morph_timesteps = 0; tolerance = 1.f; } if (SET_CAPILLARY_NUMBER && timestep%analysis_interval == analysis_interval-20 ){ - lead_timesteps += analysis_interval; + morph_timesteps += analysis_interval; double vA_x = Averages->van_global(0); double vA_y = Averages->van_global(1); double vA_z = Averages->van_global(2); @@ -600,7 +600,7 @@ void ScaLBL_ColorModel::Run(){ double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); if (rank == 0) printf(" Measured capillary number %f \n ",Ca); - if (lead_timesteps > 5000){ + if (morph_timesteps > 5000){ Fx *= capillary_number / Ca; Fy *= capillary_number / Ca; Fz *= capillary_number / Ca; From edbc3d50dfe0ca8467fea3f897c7bd08a7eeea14 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 4 Nov 2018 11:05:16 -0500 Subject: [PATCH 48/79] separate morph timesteps to avoid problems --- models/ColorModel.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 5d60d7bd..55351876 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -601,9 +601,9 @@ void ScaLBL_ColorModel::Run(){ if (rank == 0) printf(" Measured capillary number %f \n ",Ca); if (morph_timesteps > 5000){ - Fx *= capillary_number / Ca; - Fy *= capillary_number / Ca; - Fz *= capillary_number / Ca; + Fx *= (0.95 + 0.05*capillary_number / Ca); + Fy *= (0.95 + 0.05*capillary_number / Ca); + Fz *= (0.95 + 0.05*capillary_number / Ca); double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); @@ -622,7 +622,7 @@ void ScaLBL_ColorModel::Run(){ Fy *= 2.f; Fz *= 2.f; } - tolerance = fabs(capillary_number - Ca) / capillary_number ; + tolerance = fabs(1.f - (0.95 + 0.05*capillary_number / Ca)); if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); } Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); From 11f9fe00cdd20e525280b660093c196a3042e6cf Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 4 Nov 2018 13:18:56 -0500 Subject: [PATCH 49/79] tweaking morph adapt --- models/ColorModel.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 55351876..23ebe3fe 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -541,12 +541,11 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); - if (timestep > ramp_timesteps){ - // allow initial ramp-up to get closer to steady state - if (!MORPH_ADAPT && (morph_timesteps%morph_interval == 0 || tolerance < 0.01)){ + // allow initial ramp-up to get closer to steady state + if (timestep > ramp_timesteps && timestep%analysis_interval == analysis_interval-20 ){ + if ( morph_timesteps > morph_interval ){ tolerance = 1.f; MORPH_ADAPT = true; - TARGET_SATURATION = target_saturation[target_saturation_index++]; double volB = Averages->Volume_w(); double volA = Averages->Volume_n(); double current_saturation = volB/(volA+volB); @@ -564,7 +563,7 @@ void ScaLBL_ColorModel::Run(){ } } } - if (MORPH_ADAPT && timestep%analysis_interval == analysis_interval-20 ){ + if (MORPH_ADAPT ){ if (rank==0) printf("***Morphological step with target saturation %f ***\n",TARGET_SATURATION); double volB = Averages->Volume_w(); double volA = Averages->Volume_n(); @@ -574,12 +573,13 @@ void ScaLBL_ColorModel::Run(){ //MORPH_ADAPT = false; if (volB/(volA + volB) > TARGET_SATURATION){ MORPH_ADAPT = false; + TARGET_SATURATION = target_saturation[target_saturation_index++]; } MPI_Barrier(comm); morph_timesteps = 0; tolerance = 1.f; } - if (SET_CAPILLARY_NUMBER && timestep%analysis_interval == analysis_interval-20 ){ + if (SET_CAPILLARY_NUMBER ){ morph_timesteps += analysis_interval; double vA_x = Averages->van_global(0); double vA_y = Averages->van_global(1); From 37702347ccae57769ad445a742a9bf9f6818af58 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 4 Nov 2018 15:03:19 -0500 Subject: [PATCH 50/79] tweaking morph adapt --- models/ColorModel.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 23ebe3fe..8729bc3a 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -612,16 +612,17 @@ void ScaLBL_ColorModel::Run(){ Fy *= 1e-3/force_magnitude; Fz *= 1e-3/force_magnitude; } - if (force_magnitude < 1e-7){ - Fx *= 1e-7/force_magnitude; // impose floor - Fy *= 1e-7/force_magnitude; - Fz *= 1e-7/force_magnitude; + if (force_magnitude < 1e-6){ + Fx *= 1e-6/force_magnitude; // impose floor + Fy *= 1e-6/force_magnitude; + Fz *= 1e-6/force_magnitude; } - if (Fx*vA_x + Fy*vA_y + Fz*vA_z < 0.f || Fx*vB_x + Fy*vB_y +Fz*vB_z < 0.f){ + /* if (Fx*vA_x + Fy*vA_y + Fz*vA_z < 0.f || Fx*vB_x + Fy*vB_y +Fz*vB_z < 0.f){ Fx *= 2.f; // bigger forces are needed if flow rate is "too small" Fy *= 2.f; Fz *= 2.f; } + */ tolerance = fabs(1.f - (0.95 + 0.05*capillary_number / Ca)); if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); } From d6bcf5b67707989e06e556b9d2004568f760356b Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 6 Nov 2018 15:48:54 -0500 Subject: [PATCH 51/79] removed force regulator --- models/ColorModel.cpp | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8729bc3a..10ba6a1c 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -600,7 +600,7 @@ void ScaLBL_ColorModel::Run(){ double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); if (rank == 0) printf(" Measured capillary number %f \n ",Ca); - if (morph_timesteps > 5000){ + /*if (morph_timesteps > 5000){ Fx *= (0.95 + 0.05*capillary_number / Ca); Fy *= (0.95 + 0.05*capillary_number / Ca); Fz *= (0.95 + 0.05*capillary_number / Ca); @@ -617,16 +617,11 @@ void ScaLBL_ColorModel::Run(){ Fy *= 1e-6/force_magnitude; Fz *= 1e-6/force_magnitude; } - /* if (Fx*vA_x + Fy*vA_y + Fz*vA_z < 0.f || Fx*vB_x + Fy*vB_y +Fz*vB_z < 0.f){ - Fx *= 2.f; // bigger forces are needed if flow rate is "too small" - Fy *= 2.f; - Fz *= 2.f; - } - */ tolerance = fabs(1.f - (0.95 + 0.05*capillary_number / Ca)); if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); } Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); + */ } } } From c4b8f3831eece1aac9bbf131c598ead44de1a09a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 6 Nov 2018 16:06:27 -0500 Subject: [PATCH 52/79] moved force regulator --- models/ColorModel.cpp | 85 +++++++++++++++++++++---------------------- 1 file changed, 41 insertions(+), 44 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 10ba6a1c..6c34397c 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -546,9 +546,48 @@ void ScaLBL_ColorModel::Run(){ if ( morph_timesteps > morph_interval ){ tolerance = 1.f; MORPH_ADAPT = true; + double volB = Averages->Volume_w(); double volA = Averages->Volume_n(); + double vA_x = Averages->van_global(0); + double vA_y = Averages->van_global(1); + double vA_z = Averages->van_global(2); + double vB_x = Averages->vaw_global(0); + double vB_y = Averages->vaw_global(1); + double vB_z = Averages->vaw_global(2); + double muA = rhoA*(tauA-0.5)/3.f; + double muB = rhoB*(tauB-0.5)/3.f; + + double flow_rate_A = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); + double flow_rate_B = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); double current_saturation = volB/(volA+volB); + double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); + + if (rank == 0) printf(" Measured capillary number %f \n ",Ca); + + if (SET_CAPILLARY_NUMBER ){ + Fx *= capillary_number / Ca; + Fy *= capillary_number / Ca; + Fz *= capillary_number / Ca; + + double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); + + if (force_magnitude > 1e-3){ + Fx *= 1e-3/force_magnitude; // impose ceiling for stability + Fy *= 1e-3/force_magnitude; + Fz *= 1e-3/force_magnitude; + } + if (force_magnitude < 1e-6){ + Fx *= 1e-6/force_magnitude; // impose floor + Fy *= 1e-6/force_magnitude; + Fz *= 1e-6/force_magnitude; + } + tolerance = fabs(Ca-capillary_number)/ Ca ; + if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); + + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); + } + if (morph_delta > 0.f){ // wetting phase saturation will decrease while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ @@ -562,6 +601,7 @@ void ScaLBL_ColorModel::Run(){ if (rank==0) printf(" Set target saturation as %f (currently %f)\n",TARGET_SATURATION,current_saturation); } } + } if (MORPH_ADAPT ){ if (rank==0) printf("***Morphological step with target saturation %f ***\n",TARGET_SATURATION); @@ -579,50 +619,7 @@ void ScaLBL_ColorModel::Run(){ morph_timesteps = 0; tolerance = 1.f; } - if (SET_CAPILLARY_NUMBER ){ - morph_timesteps += analysis_interval; - double vA_x = Averages->van_global(0); - double vA_y = Averages->van_global(1); - double vA_z = Averages->van_global(2); - - double vB_x = Averages->vaw_global(0); - double vB_y = Averages->vaw_global(1); - double vB_z = Averages->vaw_global(2); - - double volB = Averages->Volume_w(); - double volA = Averages->Volume_n(); - double muA = rhoA*(tauA-0.5)/3.f; - double muB = rhoB*(tauB-0.5)/3.f; - - double flow_rate_A = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); - double flow_rate_B = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); - - double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); - - if (rank == 0) printf(" Measured capillary number %f \n ",Ca); - /*if (morph_timesteps > 5000){ - Fx *= (0.95 + 0.05*capillary_number / Ca); - Fy *= (0.95 + 0.05*capillary_number / Ca); - Fz *= (0.95 + 0.05*capillary_number / Ca); - - double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); - - if (force_magnitude > 1e-3){ - Fx *= 1e-3/force_magnitude; // impose ceiling for stability - Fy *= 1e-3/force_magnitude; - Fz *= 1e-3/force_magnitude; - } - if (force_magnitude < 1e-6){ - Fx *= 1e-6/force_magnitude; // impose floor - Fy *= 1e-6/force_magnitude; - Fz *= 1e-6/force_magnitude; - } - tolerance = fabs(1.f - (0.95 + 0.05*capillary_number / Ca)); - if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); - } - Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); - */ - } + morph_timesteps += analysis_interval; } } analysis.finish(); From 17fb7ec73681759d4cc421745e3e2b1178ea37f7 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 9 Nov 2018 09:37:21 -0500 Subject: [PATCH 53/79] added printout for steady kr --- models/ColorModel.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 6c34397c..0ad3edc8 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -563,8 +563,11 @@ void ScaLBL_ColorModel::Run(){ double current_saturation = volB/(volA+volB); double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); - if (rank == 0) printf(" Measured capillary number %f \n ",Ca); + FILE * kr_log_file = fopen("relperm.csv","a"); + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g .5g %.5g %.5g\n",timestep,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); + fclose(kr_log_file); + if (rank == 0) printf(" Measured capillary number %f \n ",Ca); if (SET_CAPILLARY_NUMBER ){ Fx *= capillary_number / Ca; Fy *= capillary_number / Ca; From 77e56572587c15faa317d7c0ee7b84aaef04eafa Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 9 Nov 2018 10:57:31 -0500 Subject: [PATCH 54/79] working on steady state --- models/ColorModel.cpp | 110 +++++++++++++++++++++++++----------------- 1 file changed, 65 insertions(+), 45 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 0ad3edc8..d58164e6 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -399,10 +399,16 @@ void ScaLBL_ColorModel::Run(){ bool SET_CAPILLARY_NUMBER = false; bool MORPH_ADAPT = false; + bool USE_MORPH = false; + int morph_interval; + double morph_delta; int morph_timesteps = 0; int ramp_timesteps = 50000; + double capillary_number; double tolerance = 1.f; - + double krA_previous = 0.f; + double krB_previous = 0.f; + int target_saturation_index=0; std::vector target_saturation; double TARGET_SATURATION = 0.f; @@ -410,8 +416,6 @@ void ScaLBL_ColorModel::Run(){ target_saturation = color_db->getVector( "target_saturation" ); TARGET_SATURATION = target_saturation[0]; } - - double capillary_number; if (color_db->keyExists( "capillary_number" )){ capillary_number = color_db->getScalar( "capillary_number" ); SET_CAPILLARY_NUMBER=true; @@ -423,11 +427,9 @@ void ScaLBL_ColorModel::Run(){ if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 \n"); SET_CAPILLARY_NUMBER=false; } - - int morph_interval; - double morph_delta; if (analysis_db->keyExists( "morph_delta" )){ morph_delta = analysis_db->getScalar( "morph_delta" ); + USE_MORPH = true; } else{ morph_delta=0.5; @@ -436,11 +438,16 @@ void ScaLBL_ColorModel::Run(){ morph_interval = analysis_db->getScalar( "morph_interval" ); } else{ - morph_interval=100000; + morph_interval=10000; + } + if (analysis_db->keyExists( "tolerance" )){ + tolerance = analysis_db->getScalar( "tolerance" ); + } + else{ + tolerance = 0.02; } int analysis_interval = analysis_db->getScalar( "analysis_interval" ); - if (rank==0){ printf("********************************************************\n"); printf("No. of timesteps: %i \n", timestepMax); @@ -542,11 +549,9 @@ void ScaLBL_ColorModel::Run(){ analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); // allow initial ramp-up to get closer to steady state - if (timestep > ramp_timesteps && timestep%analysis_interval == analysis_interval-20 ){ + if (timestep > ramp_timesteps && timestep%analysis_interval == analysis_interval-20 && USE_MORPH){ if ( morph_timesteps > morph_interval ){ - tolerance = 1.f; - MORPH_ADAPT = true; - + double volB = Averages->Volume_w(); double volA = Averages->Volume_n(); double vA_x = Averages->van_global(0); @@ -563,48 +568,63 @@ void ScaLBL_ColorModel::Run(){ double current_saturation = volB/(volA+volB); double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); - FILE * kr_log_file = fopen("relperm.csv","a"); - fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g .5g %.5g %.5g\n",timestep,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); - fclose(kr_log_file); - - if (rank == 0) printf(" Measured capillary number %f \n ",Ca); - if (SET_CAPILLARY_NUMBER ){ - Fx *= capillary_number / Ca; - Fy *= capillary_number / Ca; - Fz *= capillary_number / Ca; + double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); + double krA = muA*volA*flow_rate_A/force_magnitude; + double krB = muB*volB*flow_rate_B/force_magnitude; - double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); + if (fabs(krA - krA_previous) < tolerance && fabs(krB - krB_previous) < tolerance ){ + if (rank==0) printf("** WRITE STEADY POINT *** "); + tolerance = 1.f; + MORPH_ADAPT = true; - if (force_magnitude > 1e-3){ - Fx *= 1e-3/force_magnitude; // impose ceiling for stability - Fy *= 1e-3/force_magnitude; - Fz *= 1e-3/force_magnitude; + FILE * kr_log_file = fopen("relperm.csv","a"); + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g .5g %.5g %.5g\n",timestep,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); + fclose(kr_log_file); + + if (rank == 0) printf(" Measured capillary number %f \n ",Ca); + if (SET_CAPILLARY_NUMBER ){ + Fx *= capillary_number / Ca; + Fy *= capillary_number / Ca; + Fz *= capillary_number / Ca; + + if (force_magnitude > 1e-3){ + Fx *= 1e-3/force_magnitude; // impose ceiling for stability + Fy *= 1e-3/force_magnitude; + Fz *= 1e-3/force_magnitude; + } + if (force_magnitude < 1e-6){ + Fx *= 1e-6/force_magnitude; // impose floor + Fy *= 1e-6/force_magnitude; + Fz *= 1e-6/force_magnitude; + } + tolerance = fabs(Ca-capillary_number)/ Ca ; + if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); + + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); } - if (force_magnitude < 1e-6){ - Fx *= 1e-6/force_magnitude; // impose floor - Fy *= 1e-6/force_magnitude; - Fz *= 1e-6/force_magnitude; + + if (morph_delta > 0.f){ + // wetting phase saturation will decrease + while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ + TARGET_SATURATION = target_saturation[target_saturation_index++]; + } } - tolerance = fabs(Ca-capillary_number)/ Ca ; - if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); - - Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); - } - - if (morph_delta > 0.f){ - // wetting phase saturation will decrease - while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ - TARGET_SATURATION = target_saturation[target_saturation_index++]; + else{ + // wetting phase saturation will increase + while (current_saturation > TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ + TARGET_SATURATION = target_saturation[target_saturation_index++]; + if (rank==0) printf(" Set target saturation as %f (currently %f)\n",TARGET_SATURATION,current_saturation); + } } } else{ - // wetting phase saturation will increase - while (current_saturation > TARGET_SATURATION && target_saturation_index < target_saturation.size() ){ - TARGET_SATURATION = target_saturation[target_saturation_index++]; - if (rank==0) printf(" Set target saturation as %f (currently %f)\n",TARGET_SATURATION,current_saturation); + if (rank==0){ + printf("** Continue to simulate steady *** \n "); + printf("krA = %f, krB = %f, (previous = %f, %f) \n",krA,krB,krA_previous,krB_previous); } + krA_previous = krA; + krB_previous = krB; } - } if (MORPH_ADAPT ){ if (rank==0) printf("***Morphological step with target saturation %f ***\n",TARGET_SATURATION); From c2735c58fc59f81fd30eee9464088e1111ae7d6c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 9 Nov 2018 15:25:31 -0500 Subject: [PATCH 55/79] normalize by volume in kr loop --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index d58164e6..10113fe8 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -569,8 +569,8 @@ void ScaLBL_ColorModel::Run(){ double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); - double krA = muA*volA*flow_rate_A/force_magnitude; - double krB = muB*volB*flow_rate_B/force_magnitude; + double krA = muA*volA*flow_rate_A/force_magnitude/double(Nx*Ny*Nz*nprocs); + double krB = muB*volB*flow_rate_B/force_magnitude/double(Nx*Ny*Nz*nprocs); if (fabs(krA - krA_previous) < tolerance && fabs(krB - krB_previous) < tolerance ){ if (rank==0) printf("** WRITE STEADY POINT *** "); From 05b21b65f38352f7c99664999878f64ab26cdcfe Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 9 Nov 2018 16:16:36 -0500 Subject: [PATCH 56/79] termination criterion for Ca --- models/ColorModel.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 10113fe8..4804f2e0 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -406,8 +406,7 @@ void ScaLBL_ColorModel::Run(){ int ramp_timesteps = 50000; double capillary_number; double tolerance = 1.f; - double krA_previous = 0.f; - double krB_previous = 0.f; + double Ca_previous = 0.f; int target_saturation_index=0; std::vector target_saturation; @@ -569,10 +568,10 @@ void ScaLBL_ColorModel::Run(){ double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs)); double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz); - double krA = muA*volA*flow_rate_A/force_magnitude/double(Nx*Ny*Nz*nprocs); - double krB = muB*volB*flow_rate_B/force_magnitude/double(Nx*Ny*Nz*nprocs); + //double krA = muA*volA*flow_rate_A/force_magnitude/double(Nx*Ny*Nz*nprocs); + //double krB = muB*volB*flow_rate_B/force_magnitude/double(Nx*Ny*Nz*nprocs); - if (fabs(krA - krA_previous) < tolerance && fabs(krB - krB_previous) < tolerance ){ + if (fabs((Ca - Ca_previous)/Ca) < tolerance ){ if (rank==0) printf("** WRITE STEADY POINT *** "); tolerance = 1.f; MORPH_ADAPT = true; @@ -620,10 +619,9 @@ void ScaLBL_ColorModel::Run(){ else{ if (rank==0){ printf("** Continue to simulate steady *** \n "); - printf("krA = %f, krB = %f, (previous = %f, %f) \n",krA,krB,krA_previous,krB_previous); + printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); } - krA_previous = krA; - krB_previous = krB; + Ca_previous = Ca; } } if (MORPH_ADAPT ){ From 342a02235025cfa189cd9be003cba5cd3b60443a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 9 Nov 2018 16:54:32 -0500 Subject: [PATCH 57/79] termination criterion for Ca --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 4804f2e0..715cd8f2 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -621,8 +621,8 @@ void ScaLBL_ColorModel::Run(){ printf("** Continue to simulate steady *** \n "); printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); } - Ca_previous = Ca; } + Ca_previous = Ca; } if (MORPH_ADAPT ){ if (rank==0) printf("***Morphological step with target saturation %f ***\n",TARGET_SATURATION); From d92a6c9dd9aa4d1094f1cc50b65b7d71228d54d7 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 9 Nov 2018 18:04:54 -0500 Subject: [PATCH 58/79] fixe dtolerance --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 715cd8f2..4c98a98b 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -573,7 +573,7 @@ void ScaLBL_ColorModel::Run(){ if (fabs((Ca - Ca_previous)/Ca) < tolerance ){ if (rank==0) printf("** WRITE STEADY POINT *** "); - tolerance = 1.f; + printf("Ca = %f, (previous = %f), tolerance = %f \n",Ca,Ca_previous); MORPH_ADAPT = true; FILE * kr_log_file = fopen("relperm.csv","a"); From 73c9f5fa3760ada80bd95465022f243dc677d65f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 9 Nov 2018 22:53:17 -0500 Subject: [PATCH 59/79] fix morph_interval --- models/ColorModel.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 4c98a98b..d3c6f2e0 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -573,7 +573,7 @@ void ScaLBL_ColorModel::Run(){ if (fabs((Ca - Ca_previous)/Ca) < tolerance ){ if (rank==0) printf("** WRITE STEADY POINT *** "); - printf("Ca = %f, (previous = %f), tolerance = %f \n",Ca,Ca_previous); + printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); MORPH_ADAPT = true; FILE * kr_log_file = fopen("relperm.csv","a"); @@ -621,6 +621,7 @@ void ScaLBL_ColorModel::Run(){ printf("** Continue to simulate steady *** \n "); printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); } + morph_timesteps=0; } Ca_previous = Ca; } From ba7213080c57a65b1a917a42b0cd0d17ae0c4cc8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 10 Nov 2018 07:05:13 -0500 Subject: [PATCH 60/79] fixed parlalel print --- models/ColorModel.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index d3c6f2e0..465e69fc 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -572,8 +572,10 @@ void ScaLBL_ColorModel::Run(){ //double krB = muB*volB*flow_rate_B/force_magnitude/double(Nx*Ny*Nz*nprocs); if (fabs((Ca - Ca_previous)/Ca) < tolerance ){ - if (rank==0) printf("** WRITE STEADY POINT *** "); - printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); + if (rank==0){ + printf("** WRITE STEADY POINT *** "); + printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); + } MORPH_ADAPT = true; FILE * kr_log_file = fopen("relperm.csv","a"); From a390bda2f88404ef115f9a2e44cbdf6c895e9a19 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 10 Nov 2018 13:21:26 -0500 Subject: [PATCH 61/79] fixed tolerance in color model --- models/ColorModel.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 465e69fc..b43b44d7 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -598,9 +598,7 @@ void ScaLBL_ColorModel::Run(){ Fy *= 1e-6/force_magnitude; Fz *= 1e-6/force_magnitude; } - tolerance = fabs(Ca-capillary_number)/ Ca ; - if (rank == 0) printf(" -- adjust force by %f \n ",tolerance); - + if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca); Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha); } @@ -641,7 +639,6 @@ void ScaLBL_ColorModel::Run(){ } MPI_Barrier(comm); morph_timesteps = 0; - tolerance = 1.f; } morph_timesteps += analysis_interval; } From 24dc725d041ead86ac4cc18b247833dda24b63c3 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 10 Nov 2018 15:18:09 -0500 Subject: [PATCH 62/79] fixed tolerance in color model --- models/ColorModel.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index b43b44d7..cfa16bb0 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -572,17 +572,18 @@ void ScaLBL_ColorModel::Run(){ //double krB = muB*volB*flow_rate_B/force_magnitude/double(Nx*Ny*Nz*nprocs); if (fabs((Ca - Ca_previous)/Ca) < tolerance ){ + MORPH_ADAPT = true; if (rank==0){ printf("** WRITE STEADY POINT *** "); printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); + + FILE * kr_log_file = fopen("relperm.csv","a"); + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g .5g %.5g %.5g\n",timestep,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); + fclose(kr_log_file); + + printf(" Measured capillary number %f \n ",Ca); } - MORPH_ADAPT = true; - FILE * kr_log_file = fopen("relperm.csv","a"); - fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g .5g %.5g %.5g\n",timestep,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); - fclose(kr_log_file); - - if (rank == 0) printf(" Measured capillary number %f \n ",Ca); if (SET_CAPILLARY_NUMBER ){ Fx *= capillary_number / Ca; Fy *= capillary_number / Ca; From 67ca2d51610e8039a9504bf535083d95e495c73f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 11 Nov 2018 08:51:01 -0500 Subject: [PATCH 63/79] fixed write out (steady) --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index cfa16bb0..3f94f098 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -578,7 +578,7 @@ void ScaLBL_ColorModel::Run(){ printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); FILE * kr_log_file = fopen("relperm.csv","a"); - fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g .5g %.5g %.5g\n",timestep,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",timestep-analysis_interval+20,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); fclose(kr_log_file); printf(" Measured capillary number %f \n ",Ca); From f1239ada974729fc7dd093628959ed5ca71e4a76 Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Wed, 14 Nov 2018 13:43:59 -0500 Subject: [PATCH 64/79] Updating lbpm_uCT_pp.cpp --- IO/Writer.h | 2 +- analysis/runAnalysis.cpp | 856 +++++++++++++++++++------------------- common/Array.hpp | 2 +- tests/lbpm_uCT_pp.cpp | 862 ++++++++++++++++++--------------------- 4 files changed, 840 insertions(+), 882 deletions(-) diff --git a/IO/Writer.h b/IO/Writer.h index 08ac0775..710fa0d8 100644 --- a/IO/Writer.h +++ b/IO/Writer.h @@ -24,7 +24,7 @@ namespace IO { * silo - Silo * @param[in] append Append any existing data (default is false) */ -void initialize( const std::string& path="", const std::string& format="new", bool append=false ); +void initialize( const std::string& path="", const std::string& format="silo", bool append=false ); /*! diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index af8a3915..14a4fddb 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -15,23 +15,23 @@ AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs) { - lhs = static_cast( - static_cast::type>(lhs) | - static_cast::type>(rhs) - ); - return lhs; + lhs = static_cast( + static_cast::type>(lhs) | + static_cast::type>(rhs) + ); + return lhs; } bool matches( AnalysisType x, AnalysisType y ) { - return static_cast::type>(x) & - static_cast::type>(y) != 0; + return ( static_cast::type>(x) & + static_cast::type>(y) ) != 0; } template void DeleteArray( const TYPE *p ) { - delete [] p; + delete [] p; } @@ -39,36 +39,36 @@ void DeleteArray( const TYPE *p ) class WriteRestartWorkItem: public ThreadPool::WorkItemRet { public: - WriteRestartWorkItem( const char* filename_, std::shared_ptr cDen_, std::shared_ptr cfq_, int N_ ): - filename(filename_), cDen(cDen_), cfq(cfq_), N(N_) {} - virtual void run() { - PROFILE_START("Save Checkpoint",1); - double value; - ofstream File(filename,ios::binary); - for (int n=0; n cDen_, std::shared_ptr cfq_, int N_ ): + filename(filename_), cfq(cfq_), cDen(cDen_), N(N_) {} + virtual void run() { + PROFILE_START("Save Checkpoint",1); + double value; + ofstream File(filename,ios::binary); + for (int n=0; n cfq,cDen; - // const DoubleArray& phase; - //const DoubleArray& dist; + // const DoubleArray& phase; + //const DoubleArray& dist; const int N; }; @@ -78,78 +78,78 @@ static const std::string id_map_filename = "lbpm_id_map.txt"; class BlobIdentificationWorkItem1: public ThreadPool::WorkItemRet { public: - BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, - std::shared_ptr phase_, const DoubleArray& dist_, - BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ): - timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), - phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) + BlobIdentificationWorkItem1( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, + std::shared_ptr phase_, const DoubleArray& dist_, + BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_, runAnalysis::commWrapper&& comm_ ): + timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), + phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) { } - ~BlobIdentificationWorkItem1() { } - virtual void run() { - // Compute the global blob id and compare to the previous version - PROFILE_START("Identify blobs",1); - double vF = 0.0; - double vS = -1.0; // one voxel buffer region around solid - IntArray& ids = new_index->second; - new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm); - PROFILE_STOP("Identify blobs",1); - } + ~BlobIdentificationWorkItem1() { } + virtual void run() { + // Compute the global blob id and compare to the previous version + PROFILE_START("Identify blobs",1); + double vF = 0.0; + double vS = -1.0; // one voxel buffer region around solid + IntArray& ids = new_index->second; + new_index->first = ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,*phase,dist,vF,vS,ids,comm.comm); + PROFILE_STOP("Identify blobs",1); + } private: - BlobIdentificationWorkItem1(); - int timestep; - int Nx, Ny, Nz; - const RankInfoStruct& rank_info; - std::shared_ptr phase; - const DoubleArray& dist; - BlobIDstruct last_id, new_index, new_id; - BlobIDList new_list; - runAnalysis::commWrapper comm; + BlobIdentificationWorkItem1(); + int timestep; + int Nx, Ny, Nz; + const RankInfoStruct& rank_info; + std::shared_ptr phase; + const DoubleArray& dist; + BlobIDstruct last_id, new_index, new_id; + BlobIDList new_list; + runAnalysis::commWrapper comm; }; class BlobIdentificationWorkItem2: public ThreadPool::WorkItemRet { public: - BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, - std::shared_ptr phase_, const DoubleArray& dist_, - BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ): - timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), - phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) + BlobIdentificationWorkItem2( int timestep_, int Nx_, int Ny_, int Nz_, const RankInfoStruct& rank_info_, + std::shared_ptr phase_, const DoubleArray& dist_, + BlobIDstruct last_id_, BlobIDstruct new_index_, BlobIDstruct new_id_, BlobIDList new_list_ , runAnalysis::commWrapper&& comm_ ): + timestep(timestep_), Nx(Nx_), Ny(Ny_), Nz(Nz_), rank_info(rank_info_), + phase(phase_), dist(dist_), last_id(last_id_), new_index(new_index_), new_id(new_id_), new_list(new_list_), comm(std::move(comm_)) { } - ~BlobIdentificationWorkItem2() { } - virtual void run() { - // Compute the global blob id and compare to the previous version - PROFILE_START("Identify blobs maps",1); - const IntArray& ids = new_index->second; - static int max_id = -1; - new_id->first = new_index->first; - new_id->second = new_index->second; - if ( last_id.get()!=NULL ) { - // Compute the timestep-timestep map - const IntArray& old_ids = last_id->second; - ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm); - // Renumber the current timestep's ids - getNewIDs(map,max_id,*new_list); - renumberIDs(*new_list,new_id->second); - writeIDMap(map,timestep,id_map_filename); - } else { - max_id = -1; - ID_map_struct map(new_id->first); - getNewIDs(map,max_id,*new_list); - writeIDMap(map,timestep,id_map_filename); - } - PROFILE_STOP("Identify blobs maps",1); - } + ~BlobIdentificationWorkItem2() { } + virtual void run() { + // Compute the global blob id and compare to the previous version + PROFILE_START("Identify blobs maps",1); + const IntArray& ids = new_index->second; + static int max_id = -1; + new_id->first = new_index->first; + new_id->second = new_index->second; + if ( last_id.get()!=NULL ) { + // Compute the timestep-timestep map + const IntArray& old_ids = last_id->second; + ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,comm.comm); + // Renumber the current timestep's ids + getNewIDs(map,max_id,*new_list); + renumberIDs(*new_list,new_id->second); + writeIDMap(map,timestep,id_map_filename); + } else { + max_id = -1; + ID_map_struct map(new_id->first); + getNewIDs(map,max_id,*new_list); + writeIDMap(map,timestep,id_map_filename); + } + PROFILE_STOP("Identify blobs maps",1); + } private: - BlobIdentificationWorkItem2(); - int timestep; - int Nx, Ny, Nz; - const RankInfoStruct& rank_info; - std::shared_ptr phase; - const DoubleArray& dist; - BlobIDstruct last_id, new_index, new_id; - BlobIDList new_list; - runAnalysis::commWrapper comm; + BlobIdentificationWorkItem2(); + int timestep; + int Nx, Ny, Nz; + const RankInfoStruct& rank_info; + std::shared_ptr phase; + const DoubleArray& dist; + BlobIDstruct last_id, new_index, new_id; + BlobIDList new_list; + runAnalysis::commWrapper comm; }; @@ -190,12 +190,12 @@ public: PROFILE_STOP("Save Vis",1); }; private: - WriteVisWorkItem(); - int timestep; - std::vector& visData; - TwoPhase& Averages; - fillHalo& fillData; - runAnalysis::commWrapper comm; + WriteVisWorkItem(); + int timestep; + std::vector& visData; + TwoPhase& Averages; + fillHalo& fillData; + runAnalysis::commWrapper comm; }; @@ -204,46 +204,46 @@ private: class AnalysisWorkItem: public ThreadPool::WorkItemRet { public: - AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, - BlobIDstruct ids, BlobIDList id_list_, double beta_ ): - type(type_), timestep(timestep_), Averages(Averages_), - blob_ids(ids), id_list(id_list_), beta(beta_) { } - ~AnalysisWorkItem() { } - virtual void run() { - Averages.NumberComponents_NWP = blob_ids->first; - Averages.Label_NWP = blob_ids->second; - Averages.Label_NWP_map = *id_list; - Averages.NumberComponents_WP = 1; - Averages.Label_WP.fill(0.0); - if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { - // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); - } - if ( matches(type,AnalysisType::ComputeAverages) ) { - PROFILE_START("Compute dist",1); - Averages.Initialize(); - Averages.ComputeDelPhi(); - Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn); - Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus); - Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus); - Averages.UpdateMeshValues(); - Averages.ComputeLocal(); - Averages.Reduce(); - Averages.PrintAll(timestep); - Averages.Initialize(); - Averages.ComponentAverages(); - Averages.SortBlobs(); - Averages.PrintComponents(timestep); - PROFILE_STOP("Compute dist",1); - } - } + AnalysisWorkItem( AnalysisType type_, int timestep_, TwoPhase& Averages_, + BlobIDstruct ids, BlobIDList id_list_, double beta_ ): + type(type_), timestep(timestep_), Averages(Averages_), + blob_ids(ids), id_list(id_list_), beta(beta_) { } + ~AnalysisWorkItem() { } + virtual void run() { + Averages.NumberComponents_NWP = blob_ids->first; + Averages.Label_NWP = blob_ids->second; + Averages.Label_NWP_map = *id_list; + Averages.NumberComponents_WP = 1; + Averages.Label_WP.fill(0.0); + if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { + // Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus); + } + if ( matches(type,AnalysisType::ComputeAverages) ) { + PROFILE_START("Compute dist",1); + Averages.Initialize(); + Averages.ComputeDelPhi(); + Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.SDn); + Averages.ColorToSignedDistance(beta,Averages.Phase_tminus,Averages.Phase_tminus); + Averages.ColorToSignedDistance(beta,Averages.Phase_tplus,Averages.Phase_tplus); + Averages.UpdateMeshValues(); + Averages.ComputeLocal(); + Averages.Reduce(); + Averages.PrintAll(timestep); + Averages.Initialize(); + Averages.ComponentAverages(); + Averages.SortBlobs(); + Averages.PrintComponents(timestep); + PROFILE_STOP("Compute dist",1); + } + } private: - AnalysisWorkItem(); - AnalysisType type; - int timestep; - TwoPhase& Averages; - BlobIDstruct blob_ids; - BlobIDList id_list; - double beta; + AnalysisWorkItem(); + AnalysisType type; + int timestep; + TwoPhase& Averages; + BlobIDstruct blob_ids; + BlobIDList id_list; + double beta; }; @@ -251,44 +251,44 @@ private: * MPI comm wrapper for use with analysis * ******************************************************************/ runAnalysis::commWrapper::commWrapper( int tag_, MPI_Comm comm_, runAnalysis* analysis_ ): - comm(comm_), - tag(tag_), - analysis(analysis_) + comm(comm_), + tag(tag_), + analysis(analysis_) { } runAnalysis::commWrapper::commWrapper( commWrapper &&rhs ): - comm(rhs.comm), - tag(rhs.tag), - analysis(rhs.analysis) + comm(rhs.comm), + tag(rhs.tag), + analysis(rhs.analysis) { - rhs.tag = -1; + rhs.tag = -1; } runAnalysis::commWrapper::~commWrapper() { - if ( tag == -1 ) - return; - MPI_Barrier( comm ); - analysis->d_comm_used[tag] = false; + if ( tag == -1 ) + return; + MPI_Barrier( comm ); + analysis->d_comm_used[tag] = false; } runAnalysis::commWrapper runAnalysis::getComm( ) { - // Get a tag from root - int tag = -1; - if ( d_rank == 0 ) { - for (int i=0; i<1024; i++) { - if ( !d_comm_used[i] ) { - tag = i; - break; - } - } - if ( tag == -1 ) - ERROR("Unable to get comm"); - } - MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm ); - d_comm_used[tag] = true; - if ( d_comms[tag] == MPI_COMM_NULL ) - MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] ); - return commWrapper(tag,d_comms[tag],this); + // Get a tag from root + int tag = -1; + if ( d_rank == 0 ) { + for (int i=0; i<1024; i++) { + if ( !d_comm_used[i] ) { + tag = i; + break; + } + } + if ( tag == -1 ) + ERROR("Unable to get comm"); + } + MPI_Bcast( &tag, 1, MPI_INT, 0, d_comm ); + d_comm_used[tag] = true; + if ( d_comms[tag] == MPI_COMM_NULL ) + MPI_Comm_dup( MPI_COMM_WORLD, &d_comms[tag] ); + return commWrapper(tag,d_comms[tag],this); } @@ -296,116 +296,122 @@ runAnalysis::commWrapper runAnalysis::getComm( ) * Constructor/Destructors * ******************************************************************/ runAnalysis::runAnalysis( std::shared_ptr db, - const RankInfoStruct& rank_info, std::shared_ptr ScaLBL_Comm, std::shared_ptr Dm, - int Np, bool Regular, double beta, IntArray Map ): - d_Np( Np ), - d_beta( beta ), - d_regular ( Regular), - d_rank_info( rank_info ), - d_Map( Map ), - d_ScaLBL_Comm( ScaLBL_Comm), - d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1) + const RankInfoStruct& rank_info, std::shared_ptr ScaLBL_Comm, std::shared_ptr Dm, + int Np, bool Regular, double beta, IntArray Map ): + d_Np( Np ), + d_beta( beta ), + d_regular ( Regular), + d_rank_info( rank_info ), + d_Map( Map ), + d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1), + d_ScaLBL_Comm( ScaLBL_Comm) { - char rankString[20]; - sprintf(rankString,"%05d",Dm->rank()); - d_N[0] = Dm->Nx; - d_N[1] = Dm->Ny; - d_N[2] = Dm->Nz; - d_restart_interval = db->getScalar( "restart_interval" ); - d_analysis_interval = db->getScalar( "analysis_interval" ); - d_blobid_interval = db->getScalar( "blobid_interval" ); - d_visualization_interval = db->getScalar( "visualization_interval" ); - auto restart_file = db->getScalar( "restart_file" ); - d_restartFile = restart_file + "." + rankString; - d_rank = MPI_WORLD_RANK(); - writeIDMap(ID_map_struct(),0,id_map_filename); - // Initialize IO for silo - IO::initialize("","silo","false"); - // Create the MeshDataStruct - d_meshData.resize(1); - d_meshData[0].meshName = "domain"; - d_meshData[0].mesh = std::make_shared( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); - auto PhaseVar = std::make_shared(); - auto PressVar = std::make_shared(); - auto VxVar = std::make_shared(); - auto VyVar = std::make_shared(); - auto VzVar = std::make_shared(); - auto SignDistVar = std::make_shared(); - auto BlobIDVar = std::make_shared(); - - PhaseVar->name = "phase"; - PhaseVar->type = IO::VariableType::VolumeVariable; - PhaseVar->dim = 1; - PhaseVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); - d_meshData[0].vars.push_back(PhaseVar); - PressVar->name = "Pressure"; - PressVar->type = IO::VariableType::VolumeVariable; - 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; - SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); - d_meshData[0].vars.push_back(SignDistVar); - BlobIDVar->name = "BlobID"; - BlobIDVar->type = IO::VariableType::VolumeVariable; - BlobIDVar->dim = 1; - BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); - d_meshData[0].vars.push_back(BlobIDVar); - // Initialize the comms - MPI_Comm_dup(MPI_COMM_WORLD,&d_comm); - for (int i=0; i<1024; i++) { - d_comms[i] = MPI_COMM_NULL; - d_comm_used[i] = false; - } - // Initialize the threads - int N_threads = db->getWithDefault( "N_threads", 4 ); - auto method = db->getWithDefault( "load_balance", "default" ); - createThreads( method, N_threads ); + // Ids of work items to use for dependencies + ThreadPool::thread_id_t d_wait_blobID; + ThreadPool::thread_id_t d_wait_analysis; + ThreadPool::thread_id_t d_wait_vis; + ThreadPool::thread_id_t d_wait_restart; + + char rankString[20]; + sprintf(rankString,"%05d",Dm->rank()); + d_N[0] = Dm->Nx; + d_N[1] = Dm->Ny; + d_N[2] = Dm->Nz; + d_restart_interval = db->getScalar( "restart_interval" ); + d_analysis_interval = db->getScalar( "analysis_interval" ); + d_blobid_interval = db->getScalar( "blobid_interval" ); + d_visualization_interval = db->getScalar( "visualization_interval" ); + auto restart_file = db->getScalar( "restart_file" ); + d_restartFile = restart_file + "." + rankString; + d_rank = MPI_WORLD_RANK(); + writeIDMap(ID_map_struct(),0,id_map_filename); + // Initialize IO for silo + IO::initialize("","silo","false"); + // Create the MeshDataStruct + d_meshData.resize(1); + d_meshData[0].meshName = "domain"; + d_meshData[0].mesh = std::make_shared( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); + auto PhaseVar = std::make_shared(); + auto PressVar = std::make_shared(); + auto VxVar = std::make_shared(); + auto VyVar = std::make_shared(); + auto VzVar = std::make_shared(); + auto SignDistVar = std::make_shared(); + auto BlobIDVar = std::make_shared(); + + PhaseVar->name = "phase"; + PhaseVar->type = IO::VariableType::VolumeVariable; + PhaseVar->dim = 1; + PhaseVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + d_meshData[0].vars.push_back(PhaseVar); + PressVar->name = "Pressure"; + PressVar->type = IO::VariableType::VolumeVariable; + 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; + SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + d_meshData[0].vars.push_back(SignDistVar); + BlobIDVar->name = "BlobID"; + BlobIDVar->type = IO::VariableType::VolumeVariable; + BlobIDVar->dim = 1; + BlobIDVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2); + d_meshData[0].vars.push_back(BlobIDVar); + // Initialize the comms + MPI_Comm_dup(MPI_COMM_WORLD,&d_comm); + for (int i=0; i<1024; i++) { + d_comms[i] = MPI_COMM_NULL; + d_comm_used[i] = false; + } + // Initialize the threads + int N_threads = db->getWithDefault( "N_threads", 4 ); + auto method = db->getWithDefault( "load_balance", "default" ); + createThreads( method, N_threads ); } runAnalysis::~runAnalysis( ) { - // Finish processing analysis - finish(); - // Clear internal data - MPI_Comm_free( &d_comm ); - for (int i=0; i<1024; i++) { - if ( d_comms[i] != MPI_COMM_NULL ) - MPI_Comm_free(&d_comms[i]); - } + // Finish processing analysis + finish(); + // Clear internal data + MPI_Comm_free( &d_comm ); + for (int i=0; i<1024; i++) { + if ( d_comms[i] != MPI_COMM_NULL ) + MPI_Comm_free(&d_comms[i]); + } } void runAnalysis::finish( ) { - PROFILE_START("finish"); - // Wait for the work items to finish - d_tpool.wait_pool_finished(); - // Clear the wait ids - d_wait_blobID.reset(); - d_wait_analysis.reset(); - d_wait_vis.reset(); - d_wait_restart.reset(); - // Syncronize - MPI_Barrier( d_comm ); - PROFILE_STOP("finish"); + PROFILE_START("finish"); + // Wait for the work items to finish + d_tpool.wait_pool_finished(); + // Clear the wait ids + d_wait_blobID.reset(); + d_wait_analysis.reset(); + d_wait_vis.reset(); + d_wait_restart.reset(); + // Syncronize + MPI_Barrier( d_comm ); + PROFILE_STOP("finish"); } @@ -414,50 +420,50 @@ void runAnalysis::finish( ) ******************************************************************/ void print( const std::vector& ids ) { - if ( ids.empty() ) - return; - printf("%i",ids[0]); - for (size_t i=1; i0 && timestep%50==0 ) { // Keep a few blob identifications queued up to keep the processors busy, // allowing us to track the blobs as fast as possible @@ -484,28 +490,28 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep ) type |= AnalysisType::IdentifyBlobs; } #endif */ - if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) { - // Copy the averages to the CPU (and identify blobs) - //printf("Copy sim state, timestep=%i \n",timestep); - type |= AnalysisType::CopySimState; - type |= AnalysisType::IdentifyBlobs; - } - if ( timestep%d_analysis_interval == 0 ) { - // Run the analysis - //printf("Compute averages, timestep=%i \n",timestep); - type |= AnalysisType::ComputeAverages; - } - if (timestep%d_restart_interval == 0) { - // Write the restart file - type |= AnalysisType::CreateRestart; - } - if (timestep%d_visualization_interval == 0) { - // Write the visualization data - type |= AnalysisType::WriteVis; - type |= AnalysisType::CopySimState; - type |= AnalysisType::IdentifyBlobs; - } - return type; + if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) { + // Copy the averages to the CPU (and identify blobs) + //printf("Copy sim state, timestep=%i \n",timestep); + type |= AnalysisType::CopySimState; + type |= AnalysisType::IdentifyBlobs; + } + if ( timestep%d_analysis_interval == 0 ) { + // Run the analysis + //printf("Compute averages, timestep=%i \n",timestep); + type |= AnalysisType::ComputeAverages; + } + if (timestep%d_restart_interval == 0) { + // Write the restart file + type |= AnalysisType::CreateRestart; + } + if (timestep%d_visualization_interval == 0) { + // Write the visualization data + type |= AnalysisType::WriteVis; + type |= AnalysisType::CopySimState; + type |= AnalysisType::IdentifyBlobs; + } + return type; } @@ -514,56 +520,56 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep ) * Run the analysis * ******************************************************************/ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi, - double *Pressure, double *Velocity, double *fq, double *Den) + double *Pressure, double *Velocity, double *fq, double *Den) { - int N = d_N[0]*d_N[1]*d_N[2]; + int N = d_N[0]*d_N[1]*d_N[2]; - // Check which analysis steps we need to perform - auto type = computeAnalysisType( timestep ); - if ( type == AnalysisType::AnalyzeNone ) - return; + // Check which analysis steps we need to perform + auto type = computeAnalysisType( timestep ); + if ( type == AnalysisType::AnalyzeNone ) + return; - // Check how may queued items we have - if ( d_tpool.N_queued() > 20 ) { - std::cerr << "Analysis queue is getting behind, waiting ...\n"; - finish(); - } + // Check how may queued items we have + if ( d_tpool.N_queued() > 20 ) { + std::cerr << "Analysis queue is getting behind, waiting ...\n"; + finish(); + } - PROFILE_START("run"); + PROFILE_START("run"); - // Copy the appropriate variables to the host (so we can spawn new threads) - ScaLBL_DeviceBarrier(); - PROFILE_START("Copy data to host",1); - std::shared_ptr phase; - /* if ( matches(type,AnalysisType::CopyPhaseIndicator) || - matches(type,AnalysisType::ComputeAverages) || - matches(type,AnalysisType::CopySimState) || - matches(type,AnalysisType::IdentifyBlobs) ) + // Copy the appropriate variables to the host (so we can spawn new threads) + ScaLBL_DeviceBarrier(); + PROFILE_START("Copy data to host",1); + std::shared_ptr phase; + /* if ( matches(type,AnalysisType::CopyPhaseIndicator) || + matches(type,AnalysisType::ComputeAverages) || + matches(type,AnalysisType::CopySimState) || + matches(type,AnalysisType::IdentifyBlobs) ) { - phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); - //ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double)); - // try 2 d_ScaLBL_Comm.RegulLayout(d_Map,Phi,Averages.Phase); - // memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double)); - int Nx = d_N[0]; - int Ny = d_N[1]; - int Nz = d_N[2]; - double *TmpDat; - TmpDat = new double [d_Np]; - ScaLBL_CopyToHost(&TmpDat[0],&Phi[0], d_Np*sizeof(double)); - for (int k=0; kdata()[n]=value; - } - } - } - } - delete [] TmpDat; + phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); + //ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double)); + // try 2 d_ScaLBL_Comm.RegulLayout(d_Map,Phi,Averages.Phase); + // memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double)); + int Nx = d_N[0]; + int Ny = d_N[1]; + int Nz = d_N[2]; + double *TmpDat; + TmpDat = new double [d_Np]; + ScaLBL_CopyToHost(&TmpDat[0],&Phi[0], d_Np*sizeof(double)); + for (int k=0; kdata()[n]=value; + } + } + } + } + delete [] TmpDat; } */ //if ( matches(type,AnalysisType::CopyPhaseIndicator) ) { @@ -571,57 +577,57 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi, if (d_regular) d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus); else - ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double)); + ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double)); //memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double)); } if ( timestep%d_analysis_interval == 0 ) { if (d_regular) d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus); else - ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double)); + ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double)); //memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double)); } //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 cfq,cDen; //if ( matches(type,AnalysisType::CreateRestart) ) { if (timestep%d_restart_interval==0){ - // Copy restart data to the CPU - cDen = std::shared_ptr(new double[2*d_Np],DeleteArray); - cfq = std::shared_ptr(new double[19*d_Np],DeleteArray); - ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double)); - ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double)); + // Copy restart data to the CPU + cDen = std::shared_ptr(new double[2*d_Np],DeleteArray); + cfq = std::shared_ptr(new double[19*d_Np],DeleteArray); + ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double)); + ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double)); } PROFILE_STOP("Copy data to host",1); // Spawn threads to do blob identification work if ( matches(type,AnalysisType::IdentifyBlobs) ) { - phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); - if (d_regular) + phase = std::shared_ptr(new DoubleArray(d_N[0],d_N[1],d_N[2])); + if (d_regular) d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase); - else - ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double)); + else + ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double)); BlobIDstruct new_index(new std::pair(0,IntArray())); BlobIDstruct new_ids(new std::pair(0,IntArray())); @@ -642,11 +648,11 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi, //if (timestep%d_restart_interval==0){ // if ( matches(type,AnalysisType::ComputeAverages) ) { if ( timestep%d_analysis_interval == 0 ) { - auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta); - work->add_dependency(d_wait_blobID); - work->add_dependency(d_wait_analysis); - work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying - d_wait_analysis = d_tpool.add_work(work); + auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta); + work->add_dependency(d_wait_blobID); + work->add_dependency(d_wait_analysis); + work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying + d_wait_analysis = d_tpool.add_work(work); } // Spawn a thread to write the restart file diff --git a/common/Array.hpp b/common/Array.hpp index f7d4398c..cdf58125 100644 --- a/common/Array.hpp +++ b/common/Array.hpp @@ -1002,7 +1002,7 @@ Array Array::coarsen( throw std::invalid_argument( "Array must be multiple of filter size" ); } Array y( S2 ); - if ( d_size.ndim() <= 3 ) + if ( d_size.ndim() > 3 ) throw std::logic_error( "Function programmed for more than 3 dimensions" ); const auto &Nh = filter.d_size; for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) { diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp index ccb0f5af..80dbb52b 100644 --- a/tests/lbpm_uCT_pp.cpp +++ b/tests/lbpm_uCT_pp.cpp @@ -30,463 +30,415 @@ int main(int argc, char **argv) { - // Initialize MPI - int rank, nprocs; - MPI_Init(&argc,&argv); - MPI_Comm comm = MPI_COMM_WORLD; - MPI_Comm_rank(comm,&rank); - MPI_Comm_size(comm,&nprocs); - { - Utilities::setErrorHandlers(); - PROFILE_START("Main"); - - //std::vector filenames; - if ( argc<2 ) { - if ( rank == 0 ){ - printf("At least one filename must be specified\n"); - } - return 1; - } - std::string filename = std::string(argv[1]); - if ( rank == 0 ){ - printf("Input data file: %s\n",filename.c_str()); - } - - auto db = std::make_shared( filename ); - auto domain_db = db->getDatabase( "Domain" ); - auto uct_db = db->getDatabase( "uCT" ); - auto analysis_db = db->getDatabase( "Analysis" ); - - // Read domain values - auto L = domain_db->getVector( "L" ); - auto size = domain_db->getVector( "n" ); - auto nproc = domain_db->getVector( "nproc" ); - int BoundaryCondition = domain_db->getScalar( "BC" ); - int nx = size[0]; - int ny = size[1]; - int nz = size[2]; - double Lx = L[0]; - double Ly = L[1]; - double Lz = L[2]; - int nprocx = nproc[0]; - int nprocy = nproc[1]; - int nprocz = nproc[2]; - - auto InputFile=uct_db->getScalar( "InputFile" ); - auto target=uct_db->getScalar("target"); - auto background=uct_db->getScalar("background"); - auto rough_cutoff=uct_db->getScalar( "rough_cutoff" ); - auto lamda=uct_db->getScalar( "lamda" ); - auto nlm_sigsq=uct_db->getScalar( "nlm_sigsq" ); - auto nlm_depth=uct_db->getScalar( "nlm_depth" ); - auto cx=uct_db->getScalar( "center_x" ); - auto cy=uct_db->getScalar( "center_y" ); - auto cz=uct_db->getScalar( "center_z" ); - auto CylRad=uct_db->getScalar( "cylinder_radius" ); - - //....................................................................... - // Reading the domain information file - //....................................................................... - // std::shared_ptr Dm (); - //for (int i=0; iNx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = 1; - //Dm->CommInit(); - - // Check that the number of processors >= the number of ranks - if ( rank==0 ) { - printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); - printf("Number of MPI ranks used: %i \n", nprocs); - printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz); - } - if ( nprocs < nprocx*nprocy*nprocz ){ - ERROR("Insufficient number of processors"); - } - - // Determine the maximum number of levels for the desired coarsen ratio - int ratio[3] = {2,2,2}; - //std::vector ratio = {4,4,4}; - // need to set up databases for each level of the mesh - std:vector multidomain_db; - - std::vector Nx(1,nx), Ny(1,ny), Nz(1,nz); - while ( Nx.back()%ratio[0]==0 && Nx.back()>8 && - Ny.back()%ratio[1]==0 && Ny.back()>8 && - Nz.back()%ratio[2]==0 && Nz.back()>8 ) - { - Nx.push_back( Nx.back()/ratio[0] ); - Ny.push_back( Ny.back()/ratio[1] ); - Nz.push_back( Nz.back()/ratio[2] ); - // clone the domain and create coarse version based on Nx,Ny,Nz - //multidomain_db.push_back(); - } - int N_levels = Nx.size(); - - // Initialize the domain - std::vector> Dm(N_levels); - for (int i=0; iid[n] = 1; - } - Dm[i]->CommInit(); - } - - // array containing a distance mask - Array MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2); - - // Create the level data - std::vector> ID(N_levels); - std::vector> LOCVOL(N_levels); - std::vector> Dist(N_levels); - std::vector> MultiScaleSmooth(N_levels); - std::vector> Mean(N_levels); - std::vector> NonLocalMean(N_levels); - std::vector>> fillDouble(N_levels); - std::vector>> fillFloat(N_levels); - std::vector>> fillChar(N_levels); - for (int i=0; i(Nx[i]+2,Ny[i]+2,Nz[i]+2); - LOCVOL[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); - Dist[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); - MultiScaleSmooth[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); - Mean[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); - NonLocalMean[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); - ID[i].fill(0); - LOCVOL[i].fill(0); - Dist[i].fill(0); - MultiScaleSmooth[i].fill(0); - Mean[i].fill(0); - NonLocalMean[i].fill(0); - fillDouble[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) ); - fillFloat[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) ); - fillChar[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) ); - } - - // Read the subvolume of interest on each processor - PROFILE_START("ReadVolume"); - int fid = netcdf::open(InputFile,netcdf::READ); - std::string varname("VOLUME"); - netcdf::VariableType type = netcdf::getVarType( fid, varname ); - std::vector dim = netcdf::getVarDim( fid, varname ); - if ( rank == 0 ) { - printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str()); - printf(" dims = %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2])); - } - { - RankInfoStruct info( rank, nprocx, nprocy, nprocz ); - int x = info.ix*nx; - int y = info.jy*ny; - int z = info.kz*nz; - // Read the local data - Array VOLUME = netcdf::getVar( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} ); - // Copy the data and fill the halos - LOCVOL[0].fill(0); - fillFloat[0]->copy( VOLUME, LOCVOL[0] ); - fillFloat[0]->fill( LOCVOL[0] ); - } - netcdf::close( fid ); - MPI_Barrier(comm); - PROFILE_STOP("ReadVolume"); - if (rank==0) printf("Read complete\n"); - - - // Filter the original data - filter_src( *Dm[0], LOCVOL[0] ); - - // Set up the mask to be distance to cylinder (crop outside cylinder) - // float CylRad=900; - for (int k=0;kiproc(); - int jproc = Dm[0]->jproc(); - int kproc = Dm[0]->kproc(); - - int x=iproc*Nx[0]+i-1; - int y=jproc*Ny[0]+j-1; - int z=kproc*Nz[0]+k-1; - - //int cx = 0.5*nprocx*Nx[0]; - //int cy = 0.5*nprocy*Ny[0]; - //int cz = 0.5*nprocz*Nz[0]; - - // distance from the center line - MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) ); - - } - } - } - - // Compute the means for the high/low regions - // (should use automated mixture model to approximate histograms) - //float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) ); - float THRESHOLD=0.5*float(target+background); - float mean_plus=0; - float mean_minus=0; - int count_plus=0; - int count_minus=0; - for (int k=1;k 0.f ){ - auto tmp = LOCVOL[0](i,j,k); - /* if ((tmp-background)*(tmp-target) > 0){ - // direction to background / target is the same - if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background - else tmp=target; // tmp closer to target - } - */ - if ( tmp > THRESHOLD ) { - mean_plus += tmp; - count_plus++; - } else if ( tmp < -THRESHOLD ) { - mean_minus += tmp; - count_minus++; - } - } - } - } - } - mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus ); - mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus ); - if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus); - - - MPI_Barrier(comm); - // Scale the source data to +-1.0 - for (size_t i=0; i= 0 ) { - LOCVOL[0](i) /= mean_plus; - LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f ); - } else { - LOCVOL[0](i) /= -mean_minus; - LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f ); - } - } - - // Fill the source data for the coarse meshes - PROFILE_START("CoarsenMesh"); - for (int i=1; i filter(ratio[0],ratio[1],ratio[2]); - filter.fill(1.0f/filter.length()); - Array tmp(Nx[i-1],Ny[i-1],Nz[i-1]); - fillFloat[i-1]->copy( LOCVOL[i-1], tmp ); - Array coarse = tmp.coarsen( filter ); - fillFloat[i]->copy( coarse, LOCVOL[i] ); - fillFloat[i]->fill( LOCVOL[i] ); - if (rank==0){ - printf("Coarsen level %i \n",i); - printf(" Nx=%i, Ny=%i, Nz=%i \n",int(tmp.size(0)),int(tmp.size(1)),int(tmp.size(2)) ); - printf(" filter_x=%i, filter_y=%i, filter_z=%i \n",int(filter.size(0)),int(filter.size(1)),int(filter.size(2)) ); - printf(" ratio= %i,%i,%i \n",int(ratio[0]),int(ratio[1]),int(ratio[2]) ); - } - MPI_Barrier(comm); - } - PROFILE_STOP("CoarsenMesh"); - - // Initialize the coarse level - PROFILE_START("Solve full mesh"); - if (rank==0) - printf("Initialize full mesh\n"); - solve( LOCVOL[0], Mean[0], ID[0], Dist[0], MultiScaleSmooth[0], - NonLocalMean[0], *fillFloat[0], *Dm[0], nprocx, - rough_cutoff, lamda, nlm_sigsq, nlm_depth); - PROFILE_STOP("Solve full mesh"); - MPI_Barrier(comm); - - /* // Initialize the coarse level - PROFILE_START("Solve coarse mesh"); - if (rank==0) - printf("Initialize coarse mesh\n"); - solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(), - NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx ); - PROFILE_STOP("Solve coarse mesh"); - MPI_Barrier(comm); - - // Refine the solution - PROFILE_START("Refine distance"); - if (rank==0) - printf("Refine mesh\n"); - for (int i=int(Nx.size())-2; i>=0; i--) { - if (rank==0) - printf(" Refining to level %i\n",int(i)); - refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i], - NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i ); - } - PROFILE_STOP("Refine distance"); - MPI_Barrier(comm); - - // Perform a final filter - PROFILE_START("Filtering final domains"); - if (rank==0) - printf("Filtering final domains\n"); - Array filter_Mean, filter_Dist1, filter_Dist2; - filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 ); - PROFILE_STOP("Filtering final domains"); - */ - //removeDisconnected( ID[0], *Dm[0] ); - - /* - // Write the distance function to a netcdf file - const char* netcdf_filename = "Distance.nc"; + // Initialize MPI + int rank, nprocs; + MPI_Init(&argc,&argv); + MPI_Comm comm = MPI_COMM_WORLD; + MPI_Comm_rank(comm,&rank); + MPI_Comm_size(comm,&nprocs); { - RankInfoStruct info( rank, nprocx, nprocy, nprocz ); - std::vector dim = { Nx[0]*nprocx, Ny[0]*nprocy, Nz[0]*nprocz }; - int fid = netcdf::open( netcdf_filename, netcdf::CREATE, MPI_COMM_WORLD ); - auto dims = netcdf::defDim( fid, {"X", "Y", "Z"}, dim ); - Array data(Nx[0],Ny[0],Nz[0]); - fillFloat[0]->copy( Dist[0], data ); - netcdf::write( fid, "Distance", dims, data, info ); + Utilities::setErrorHandlers(); + PROFILE_START("Main"); + + //std::vector filenames; + if ( argc<2 ) { + if ( rank == 0 ){ + printf("At least one filename must be specified\n"); + } + return 1; + } + std::string filename = std::string(argv[1]); + if ( rank == 0 ){ + printf("Input data file: %s\n",filename.c_str()); + } + + auto db = std::make_shared( filename ); + auto domain_db = db->getDatabase( "Domain" ); + auto uct_db = db->getDatabase( "uCT" ); + auto analysis_db = db->getDatabase( "Analysis" ); + + // Read domain values + auto L = domain_db->getVector( "L" ); + auto size = domain_db->getVector( "n" ); + auto nproc = domain_db->getVector( "nproc" ); + int BoundaryCondition = domain_db->getScalar( "BC" ); + int nx = size[0]; + int ny = size[1]; + int nz = size[2]; + double Lx = L[0]; + double Ly = L[1]; + double Lz = L[2]; + int nprocx = nproc[0]; + int nprocy = nproc[1]; + int nprocz = nproc[2]; + + auto InputFile = uct_db->getScalar( "InputFile" ); + auto target = uct_db->getScalar("target"); + auto background = uct_db->getScalar("background"); + auto rough_cutoff = uct_db->getScalar( "rough_cutoff" ); + auto lamda = uct_db->getScalar( "lamda" ); + auto nlm_sigsq = uct_db->getScalar( "nlm_sigsq" ); + auto nlm_depth = uct_db->getScalar( "nlm_depth" ); + auto center = uct_db->getVector( "center" ); + auto CylRad = uct_db->getScalar( "cylinder_radius" ); + auto maxLevels = uct_db->getScalar( "max_levels" ); + std::vector offset( 3, 0 ); + if ( uct_db->keyExists( "offset" ) ) + offset = uct_db->getVector( "offset" ); + + // Check that the number of processors >= the number of ranks + if ( rank==0 ) { + printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); + printf("Number of MPI ranks used: %i \n", nprocs); + printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz); + } + if ( nprocs < nprocx*nprocy*nprocz ){ + ERROR("Insufficient number of processors"); + } + + // Determine the maximum number of levels for the desired coarsen ratio + int ratio[3] = {2,2,2}; + //std::vector ratio = {4,4,4}; + // need to set up databases for each level of the mesh + std::vector> multidomain_db(1,domain_db); + std::vector Nx(1,nx), Ny(1,ny), Nz(1,nz); + while ( Nx.back()%ratio[0]==0 && Nx.back()>8 && + Ny.back()%ratio[1]==0 && Ny.back()>8 && + Nz.back()%ratio[2]==0 && Nz.back()>8 && + (int) Nx.size() < maxLevels ) + { + Nx.push_back( Nx.back()/ratio[0] ); + Ny.push_back( Ny.back()/ratio[1] ); + Nz.push_back( Nz.back()/ratio[2] ); + // clone the domain and create coarse version based on Nx,Ny,Nz + auto db2 = domain_db->cloneDatabase(); + db2->putVector( "n", { Nx.back(), Ny.back(), Nz.back() } ); + multidomain_db.push_back(db2); + } + int N_levels = Nx.size(); + + // Initialize the domain + std::vector> Dm(N_levels); + for (int i=0; iid[n] = 1; + } + Dm[i]->CommInit(); + } + + // array containing a distance mask + Array MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2); + MASK.fill(0); + + // Create the level data + std::vector> ID(N_levels); + std::vector> LOCVOL(N_levels); + std::vector> Dist(N_levels); + std::vector> MultiScaleSmooth(N_levels); + std::vector> Mean(N_levels); + std::vector> NonLocalMean(N_levels); + std::vector>> fillDouble(N_levels); + std::vector>> fillFloat(N_levels); + std::vector>> fillChar(N_levels); + for (int i=0; i(Nx[i]+2,Ny[i]+2,Nz[i]+2); + LOCVOL[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + Dist[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + MultiScaleSmooth[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + Mean[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + NonLocalMean[i] = Array(Nx[i]+2,Ny[i]+2,Nz[i]+2); + ID[i].fill(0); + LOCVOL[i].fill(0); + Dist[i].fill(0); + MultiScaleSmooth[i].fill(0); + Mean[i].fill(0); + NonLocalMean[i].fill(0); + fillDouble[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) ); + fillFloat[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) ); + fillChar[i].reset(new fillHalo(Dm[i]->Comm,Dm[i]->rank_info,{Nx[i],Ny[i],Nz[i]},{1,1,1},0,1) ); + } + + // Read the subvolume of interest on each processor + PROFILE_START("ReadVolume"); + int fid = netcdf::open(InputFile,netcdf::READ); + std::string varname("VOLUME"); + auto type = netcdf::getVarType( fid, varname ); + auto dim = netcdf::getVarDim( fid, varname ); + if ( rank == 0 ) { + printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str()); + printf(" dims = %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2])); + } + { + // Read the local data + int x = Dm[0]->iproc()*nx + offset[0]; + int y = Dm[0]->jproc()*ny + offset[1]; + int z = Dm[0]->kproc()*nz + offset[2]; + Array VOLUME = netcdf::getVar( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} ); + // Copy the data and fill the halos + LOCVOL[0].fill(0); + fillFloat[0]->copy( VOLUME, LOCVOL[0] ); + fillFloat[0]->fill( LOCVOL[0] ); + } netcdf::close( fid ); + MPI_Barrier(comm); + PROFILE_STOP("ReadVolume"); + if (rank==0) printf("Read complete\n"); + + + // Filter the original data + filter_src( *Dm[0], LOCVOL[0] ); + + // Set up the mask to be distance to cylinder (crop outside cylinder) + for (int k=0;kiproc()*Nx[0]+i-1; + int y=Dm[0]->jproc()*Ny[0]+j-1; + int z=Dm[0]->kproc()*Nz[0]+k-1; + int cx = center[0] - offset[0]; + int cy = center[1] - offset[1]; + int cz = center[2] - offset[2]; + // distance from the center line + MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) ); + } + } + } + + // Compute the means for the high/low regions + // (should use automated mixture model to approximate histograms) + //float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) ); + double THRESHOLD=0.5*(target+background); + double mean_plus=0; + double mean_minus=0; + int count_plus=0; + int count_minus=0; + for (int k=1;k 0.f ){ + auto tmp = LOCVOL[0](i,j,k); + /* if ((tmp-background)*(tmp-target) > 0){ + // direction to background / target is the same + if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background + else tmp=target; // tmp closer to target + } + */ + if ( tmp > THRESHOLD ) { + mean_plus += tmp; + count_plus++; + } else if ( tmp < -THRESHOLD ) { + mean_minus += tmp; + count_minus++; + } + } + } + } + } + ASSERT( count_plus > 0 && count_minus > 0 ); + mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus ); + mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus ); + if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus); + + + MPI_Barrier(comm); + // Scale the source data to +-1.0 + for (size_t i=0; i= 0 ) { + LOCVOL[0](i) /= mean_plus; + LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f ); + } else { + LOCVOL[0](i) /= -mean_minus; + LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f ); + } + } + + + // Fill the source data for the coarse meshes + PROFILE_START("CoarsenMesh"); + for (int i=1; i filter(ratio[0],ratio[1],ratio[2]); + filter.fill(1.0f/filter.length()); + Array tmp(Nx[i-1],Ny[i-1],Nz[i-1]); + fillFloat[i-1]->copy( LOCVOL[i-1], tmp ); + Array coarse = tmp.coarsen( filter ); + fillFloat[i]->copy( coarse, LOCVOL[i] ); + fillFloat[i]->fill( LOCVOL[i] ); + if (rank==0){ + printf("Coarsen level %i \n",i); + printf(" Nx=%i, Ny=%i, Nz=%i \n",int(tmp.size(0)),int(tmp.size(1)),int(tmp.size(2)) ); + printf(" filter_x=%i, filter_y=%i, filter_z=%i \n",int(filter.size(0)),int(filter.size(1)),int(filter.size(2)) ); + printf(" ratio= %i,%i,%i \n",int(ratio[0]),int(ratio[1]),int(ratio[2]) ); + } + MPI_Barrier(comm); + } + PROFILE_STOP("CoarsenMesh"); + + // Initialize the coarse level + PROFILE_START("Solve coarse mesh"); + if (rank==0) + printf("Initialize full mesh\n"); + solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(), + NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx, + rough_cutoff, lamda, nlm_sigsq, nlm_depth); + PROFILE_STOP("Solve coarse mesh"); + MPI_Barrier(comm); + + // Refine the solution + PROFILE_START("Refine distance"); + if (rank==0) + printf("Refine mesh\n"); + for (int i=N_levels-2; i>=0; i--) { + if (rank==0) + printf(" Refining to level %i\n",i); + refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i], + NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i, + rough_cutoff, lamda, nlm_sigsq, nlm_depth); + } + PROFILE_STOP("Refine distance"); + MPI_Barrier(comm); + + // Perform a final filter + PROFILE_START("Filtering final domains"); + if (rank==0) + printf("Filtering final domains\n"); + Array filter_Mean, filter_Dist1, filter_Dist2; + filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 ); + PROFILE_STOP("Filtering final domains"); + //removeDisconnected( ID[0], *Dm[0] ); + + // Write the distance function to a netcdf file + /* const char* netcdf_filename = "Distance.nc"; + { + RankInfoStruct info( rank, nprocx, nprocy, nprocz ); + std::vector dim = { Nx[0]*nprocx, Ny[0]*nprocy, Nz[0]*nprocz }; + int fid = netcdf::open( netcdf_filename, netcdf::CREATE, MPI_COMM_WORLD ); + auto dims = netcdf::defDim( fid, {"X", "Y", "Z"}, dim ); + Array data(Nx[0],Ny[0],Nz[0]); + fillFloat[0]->copy( Dist[0], data ); + netcdf::write( fid, "Distance", dims, data, info ); + netcdf::close( fid ); + } */ + + // Write the results + if (rank==0) printf("Setting up visualization structure \n"); + std::vector meshData(N_levels); + for (size_t i=0; i(Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],Lx,Ly,Lz); + // Source data + auto OrigData = std::make_shared(); + OrigData->name = "Source_Data_" + std::to_string( i ); + OrigData->type = IO::VariableType::VolumeVariable; + OrigData->dim = 1; + OrigData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(OrigData); + fillDouble[i]->copy( LOCVOL[i], OrigData->data ); + // Non-Local Mean + auto NonLocMean = std::make_shared(); + NonLocMean->name = "NonLocal_Mean_" + std::to_string( i ); + NonLocMean->type = IO::VariableType::VolumeVariable; + NonLocMean->dim = 1; + NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(NonLocMean); + fillDouble[i]->copy( NonLocalMean[i], NonLocMean->data ); + // Segmented Data + auto SegData = std::make_shared(); + SegData->name = "Segmented_Data_" + std::to_string( i ); + SegData->type = IO::VariableType::VolumeVariable; + SegData->dim = 1; + SegData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(SegData); + fillDouble[i]->copy( ID[i], SegData->data ); + // Signed Distance + auto DistData = std::make_shared(); + DistData->name = "Signed_Distance_" + std::to_string( i ); + DistData->type = IO::VariableType::VolumeVariable; + DistData->dim = 1; + DistData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(DistData); + fillDouble[i]->copy( Dist[i], DistData->data ); + // Smoothed Data + auto SmoothData = std::make_shared(); + SmoothData->name = "Smoothed_Data_" + std::to_string( i ); + SmoothData->type = IO::VariableType::VolumeVariable; + SmoothData->dim = 1; + SmoothData->data.resize(Nx[i],Ny[i],Nz[i]); + meshData[i].vars.push_back(SmoothData); + fillDouble[i]->copy( MultiScaleSmooth[i], SmoothData->data ); + + } + #if 0 + std::shared_ptr filter_Mean_var( new IO::Variable() ); + filter_Mean_var->name = "Mean"; + filter_Mean_var->type = IO::VariableType::VolumeVariable; + filter_Mean_var->dim = 1; + filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]); + meshData[0].vars.push_back(filter_Mean_var); + fillDouble[0]->copy( filter_Mean, filter_Mean_var->data ); + std::shared_ptr filter_Dist1_var( new IO::Variable() ); + filter_Dist1_var->name = "Dist1"; + filter_Dist1_var->type = IO::VariableType::VolumeVariable; + filter_Dist1_var->dim = 1; + filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]); + meshData[0].vars.push_back(filter_Dist1_var); + fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data ); + std::shared_ptr filter_Dist2_var( new IO::Variable() ); + filter_Dist2_var->name = "Dist2"; + filter_Dist2_var->type = IO::VariableType::VolumeVariable; + filter_Dist2_var->dim = 1; + filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]); + meshData[0].vars.push_back(filter_Dist2_var); + fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data ); + #endif + MPI_Barrier(comm); + if (rank==0) printf("Writing output \n"); + // Write visulization data + IO::writeData( 0, meshData, comm ); + if (rank==0) printf("Finished. \n"); + + // Compute the Minkowski functionals + MPI_Barrier(comm); + auto Averages = std::make_shared(Dm[0]); + + Array phase_label(Nx[0]+2,Ny[0]+2,Nz[0]+2); + Array phase_distance(Nx[0]+2,Ny[0]+2,Nz[0]+2); + // Analyze the wetting fluid + for (int k=1;kid[n] > 0)){ + // Solid phase + phase_label(i,j,k) = 0; + } + else if (Dist[0](i,j,k) < 0.0){ + // wetting phase + phase_label(i,j,k) = 1; + } + else { + // non-wetting phase + phase_label(i,j,k) = 0; + } + phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0; + } + } + } + CalcDist(phase_distance,phase_label,*Dm[0]); + Averages->ComputeScalar(phase_distance,0.f); + Averages->PrintAll(); } - */ - - { - // Write the results - if (rank==0) printf("Setting up visualization structure \n"); - // std::vector meshData(N_levels); - std::vector meshData(1); - // for (size_t i=0; i( new IO::DomainMesh(Dm[0]->rank_info,Nx[0],Ny[0],Nz[0],Lx,Ly,Lz) ); - // Source data - std::shared_ptr OrigData( new IO::Variable() ); - OrigData->name = "Source Data"; - OrigData->type = IO::VariableType::VolumeVariable; - OrigData->dim = 1; - OrigData->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(OrigData); - fillDouble[0]->copy( LOCVOL[0], OrigData->data ); - // Non-Local Mean - std::shared_ptr NonLocMean( new IO::Variable() ); - NonLocMean->name = "Non-Local Mean"; - NonLocMean->type = IO::VariableType::VolumeVariable; - NonLocMean->dim = 1; - NonLocMean->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(NonLocMean); - fillDouble[0]->copy( NonLocalMean[0], NonLocMean->data ); - std::shared_ptr SegData( new IO::Variable() ); - SegData->name = "Segmented Data"; - SegData->type = IO::VariableType::VolumeVariable; - SegData->dim = 1; - SegData->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(SegData); - fillDouble[0]->copy( ID[0], SegData->data ); - // Signed Distance - std::shared_ptr DistData( new IO::Variable() ); - DistData->name = "Signed Distance"; - DistData->type = IO::VariableType::VolumeVariable; - DistData->dim = 1; - DistData->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(DistData); - fillDouble[0]->copy( Dist[0], DistData->data ); - // Smoothed Data - std::shared_ptr SmoothData( new IO::Variable() ); - SmoothData->name = "Smoothed Data"; - SmoothData->type = IO::VariableType::VolumeVariable; - SmoothData->dim = 1; - SmoothData->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(SmoothData); - fillDouble[0]->copy( MultiScaleSmooth[0], SmoothData->data ); - - /*// Segmented Data - std::shared_ptr SegData( new IO::Variable() ); - SegData->name = "Segmented Data"; - SegData->type = IO::VariableType::VolumeVariable; - SegData->dim = 1; - SegData->data.resize(Nx[i],Ny[i],Nz[i]); - meshData[i].vars.push_back(SegData); - fillDouble[i]->copy( ID[i], SegData->data ); - // Signed Distance - std::shared_ptr DistData( new IO::Variable() ); - DistData->name = "Signed Distance"; - DistData->type = IO::VariableType::VolumeVariable; - DistData->dim = 1; - DistData->data.resize(Nx[i],Ny[i],Nz[i]); - meshData[i].vars.push_back(DistData); - fillDouble[i]->copy( Dist[i], DistData->data ); - // Smoothed Data - std::shared_ptr SmoothData( new IO::Variable() ); - SmoothData->name = "Smoothed Data"; - SmoothData->type = IO::VariableType::VolumeVariable; - SmoothData->dim = 1; - SmoothData->data.resize(Nx[i],Ny[i],Nz[i]); - meshData[i].vars.push_back(SmoothData); - fillDouble[i]->copy( MultiScaleSmooth[i], SmoothData->data ); - - } - #if 0 - std::shared_ptr filter_Mean_var( new IO::Variable() ); - filter_Mean_var->name = "Mean"; - filter_Mean_var->type = IO::VariableType::VolumeVariable; - filter_Mean_var->dim = 1; - filter_Mean_var->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(filter_Mean_var); - fillDouble[0]->copy( filter_Mean, filter_Mean_var->data ); - std::shared_ptr filter_Dist1_var( new IO::Variable() ); - filter_Dist1_var->name = "Dist1"; - filter_Dist1_var->type = IO::VariableType::VolumeVariable; - filter_Dist1_var->dim = 1; - filter_Dist1_var->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(filter_Dist1_var); - fillDouble[0]->copy( filter_Dist1, filter_Dist1_var->data ); - std::shared_ptr filter_Dist2_var( new IO::Variable() ); - filter_Dist2_var->name = "Dist2"; - filter_Dist2_var->type = IO::VariableType::VolumeVariable; - filter_Dist2_var->dim = 1; - filter_Dist2_var->data.resize(Nx[0],Ny[0],Nz[0]); - meshData[0].vars.push_back(filter_Dist2_var); - fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data ); - #endif - */ - MPI_Barrier(comm); - if (rank==0) printf("Writing output \n"); - // Write visulization data - IO::writeData( 0, meshData, comm ); - if (rank==0) printf("Finished. \n"); - } - - // Compute the Minkowski functionals - MPI_Barrier(comm); - std::shared_ptr Averages(new Minkowski(Dm[0])); - - Array phase_label(Nx[0],Ny[0],Nz[0]); - Array phase_distance(Nx[0],Ny[0],Nz[0]); - // Analyze the wetting fluid - for (int k=1;kid[n] > 0)){ - // Solid phase - phase_label(i,j,k) = 0; - } - else if (Dist[0](i,j,k) < 0.0){ - // wetting phase - phase_label(i,j,k) = 1; - } - else { - // non-wetting phase - phase_label(i,j,k) = 0; - } - phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0; - } - } - } - CalcDist(phase_distance,phase_label,*Dm[0]); - Averages->ComputeScalar(phase_distance,0.f); - Averages->PrintAll(); - } - PROFILE_STOP("Main"); - PROFILE_SAVE("lbpm_uCT_pp",true); - MPI_Barrier(comm); - MPI_Finalize(); - return 0; + PROFILE_STOP("Main"); + PROFILE_SAVE("lbpm_uCT_pp",true); + MPI_Barrier(comm); + MPI_Finalize(); + return 0; } From 8c8504b496f52ef5aa2b60d0b9d2ec1f54df3f15 Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Thu, 15 Nov 2018 14:50:07 -0500 Subject: [PATCH 65/79] Fixing bug in Array class --- common/Array.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/common/Array.hpp b/common/Array.hpp index cdf58125..6b2dd16f 100644 --- a/common/Array.hpp +++ b/common/Array.hpp @@ -1003,7 +1003,7 @@ Array Array::coarsen( } Array y( S2 ); if ( d_size.ndim() > 3 ) - throw std::logic_error( "Function programmed for more than 3 dimensions" ); + throw std::logic_error( "Function not programmed for more than 3 dimensions" ); const auto &Nh = filter.d_size; for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) { for ( size_t j1 = 0; j1 < y.d_size[1]; j1++ ) { @@ -1037,8 +1037,8 @@ Array Array::coarsen( const std::vec } Array tmp( ratio ); Array y( S2 ); - if ( d_size.ndim() <= 3 ) - throw std::logic_error( "Function programmed for more than 3 dimensions" ); + if ( d_size.ndim() > 3 ) + throw std::logic_error( "Function not programmed for more than 3 dimensions" ); for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) { for ( size_t j1 = 0; j1 < y.d_size[1]; j1++ ) { for ( size_t i1 = 0; i1 < y.d_size[0]; i1++ ) { From e0b00873e9c752ab695492829ee130058591a2b8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sat, 17 Nov 2018 19:53:09 -0500 Subject: [PATCH 66/79] Update to uCT --- analysis/uCT.cpp | 2 +- sample_scripts/config_build_eos | 17 ++++--- tests/lbpm_uCT_pp.cpp | 82 ++++++++++++++++++++------------- 3 files changed, 63 insertions(+), 38 deletions(-) diff --git a/analysis/uCT.cpp b/analysis/uCT.cpp index 315f441f..6a327432 100644 --- a/analysis/uCT.cpp +++ b/analysis/uCT.cpp @@ -148,7 +148,7 @@ void solve( const Array& VOL, Array& Mean, Array& ID, // Compute the median filter on the sparse array Med3D( VOL, Mean ); fillFloat.fill( Mean ); - segment( Mean, ID, 0.01 ); + segment( Mean, ID, threshold ); // Compute the distance using the segmented volume CalcDist( Dist, ID, Dm ); fillFloat.fill(Dist); diff --git a/sample_scripts/config_build_eos b/sample_scripts/config_build_eos index 4b382a63..f4d69f26 100755 --- a/sample_scripts/config_build_eos +++ b/sample_scripts/config_build_eos @@ -7,11 +7,13 @@ module unload PrgEnv-intel module load PrgEnv-gnu/6.0.4 module unload gcc cmake module load gcc/6.3.0 -module load cray-hdf5-parallel/1.10.0.1 +module load cray-hdf5-parallel/1.10.2.0 module load cray-netcdf-hdf5parallel module load mercurial git module load cmake3/3.6.1 +export LD_LIBRARY_PATH=/usr/lib64:/lustre/atlas/proj-shared/geo106/eos/netcdf/lib:/lustre/atlas/proj-shared/geo106/eos/zlib/lib:/lustre/atlas/proj-shared/geo106/eos/hdf5/lib:$LD_LIBRARY_PATH + export MPICH_RDMA_ENABLED_CUDA=0 echo $GNU_VERSION @@ -27,7 +29,8 @@ cmake \ -D CMAKE_C_COMPILER:PATH=cc \ -D CMAKE_CXX_COMPILER:PATH=CC \ -D CMAKE_CXX_COMPILER:PATH=CC \ - -D CXX_STD=11 \ + -D CMAKE_CXX_FLAGS="-fPIC" \ + -D CMAKE_CXX_STANDARD=14 \ -D USE_TIMER=false \ -D TIMER_DIRECTORY=${HOME}/timerutility/build/opt \ -D MPI_COMPILER:BOOL=TRUE \ @@ -35,10 +38,12 @@ cmake \ -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \ -D USE_CUDA=0 \ -D CUDA_FLAGS="-arch sm_35" \ - -DBUILD_SHARED_LIBS=OFF \ - -D USE_NETCDF=1 \ - -D NETCDF_DIRECTORY=$NETCDF_DIR \ - -D HDF5_DIRECTORY=$HDF5_DIR \ + -D USE_SILO=1 \ + -D SILO_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/silo \ + -D HDF5_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/hdf5 \ + -D HDF5_LIB=/lustre/atlas/proj-shared/geo106/eos/hdf5/lib/libhdf5.a \ + -D USE_NETCDF=1 \ + -D NETCDF_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/netcdf \ -D CMAKE_SKIP_RPATH=true \ ~/LBPM-WIA diff --git a/tests/lbpm_uCT_pp.cpp b/tests/lbpm_uCT_pp.cpp index 80dbb52b..c31768ac 100644 --- a/tests/lbpm_uCT_pp.cpp +++ b/tests/lbpm_uCT_pp.cpp @@ -73,8 +73,8 @@ int main(int argc, char **argv) int nprocz = nproc[2]; auto InputFile = uct_db->getScalar( "InputFile" ); - auto target = uct_db->getScalar("target"); - auto background = uct_db->getScalar("background"); + auto target = uct_db->getScalar("target"); + auto background = uct_db->getScalar("background"); auto rough_cutoff = uct_db->getScalar( "rough_cutoff" ); auto lamda = uct_db->getScalar( "lamda" ); auto nlm_sigsq = uct_db->getScalar( "nlm_sigsq" ); @@ -91,6 +91,10 @@ int main(int argc, char **argv) printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz); printf("Number of MPI ranks used: %i \n", nprocs); printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz); + printf("target value = %f \n",target); + printf("background value = %f \n",background); + printf("cylinder center = %i, %i, %i \n",center[0],center[1],center[2]); + printf("cylinder radius = %f \n",CylRad); } if ( nprocs < nprocx*nprocy*nprocz ){ ERROR("Insufficient number of processors"); @@ -188,22 +192,23 @@ int main(int argc, char **argv) PROFILE_STOP("ReadVolume"); if (rank==0) printf("Read complete\n"); - // Filter the original data filter_src( *Dm[0], LOCVOL[0] ); // Set up the mask to be distance to cylinder (crop outside cylinder) + if (rank==0) printf("Cropping with cylinder: %i, %i, %i, radius=%f \n",Dm[0]->nprocx()*Nx[0],Dm[0]->nprocy()*Ny[0],Dm[0]->nprocz()*Nz[0],CylRad); for (int k=0;kiproc()*Nx[0]+i-1; - int y=Dm[0]->jproc()*Ny[0]+j-1; - int z=Dm[0]->kproc()*Nz[0]+k-1; - int cx = center[0] - offset[0]; - int cy = center[1] - offset[1]; - int cz = center[2] - offset[2]; + float x= float(Dm[0]->iproc()*Nx[0]+i-1); + float y= float (Dm[0]->jproc()*Ny[0]+j-1); + float z= float(Dm[0]->kproc()*Nz[0]+k-1); + float cx = float(center[0] - offset[0]); + float cy = float(center[1] - offset[1]); + float cz = float(center[2] - offset[2]); // distance from the center line - MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) ); + MASK(i,j,k) = sqrt((z-cz)*(z-cz) + (y-cy)*(y-cy)); + //if (sqrt(((z-cz)*(z-cz) + (y-cy)*(y-cy)) ) > CylRad) LOCVOL[0](i,j,k)=background; } } } @@ -211,16 +216,21 @@ int main(int argc, char **argv) // Compute the means for the high/low regions // (should use automated mixture model to approximate histograms) //float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) ); - double THRESHOLD=0.5*(target+background); - double mean_plus=0; - double mean_minus=0; + float THRESHOLD=0.5*(target+background); + float mean_plus=0; + float mean_minus=0; + float min_value = LOCVOL[0](0); + float max_value = LOCVOL[0](0); int count_plus=0; int count_minus=0; for (int k=1;k 0.f ){ - auto tmp = LOCVOL[0](i,j,k); + + + //LOCVOL[0](i,j,k) = MASK(i,j,k); + if (MASK(i,j,k) < CylRad ){ + auto tmp = LOCVOL[0](i,j,k); /* if ((tmp-background)*(tmp-target) > 0){ // direction to background / target is the same if (fabs(tmp-target) > fabs(tmp-background)) tmp=background; // tmp closer to background @@ -230,37 +240,47 @@ int main(int argc, char **argv) if ( tmp > THRESHOLD ) { mean_plus += tmp; count_plus++; - } else if ( tmp < -THRESHOLD ) { + } + else { mean_minus += tmp; count_minus++; } - } + if (tmp < min_value) min_value = tmp; + if (tmp > max_value) max_value = tmp; + } } } } + count_plus=sumReduce( Dm[0]->Comm, count_plus); + count_minus=sumReduce( Dm[0]->Comm, count_minus); + if (rank==0) printf("minimum value=%f, max value=%f \n",min_value,max_value); + if (rank==0) printf("plus=%i, minus=%i \n",count_plus,count_minus); ASSERT( count_plus > 0 && count_minus > 0 ); - mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus ); - mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus ); + MPI_Barrier(comm); + mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / count_plus; + mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / count_minus; + MPI_Barrier(comm); if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus); - - MPI_Barrier(comm); - // Scale the source data to +-1.0 + //if (rank==0) printf("Scale the input data (size = %i) \n",LOCVOL[0].length()); for (size_t i=0; i CylRad ){ + LOCVOL[0](i)=background; } - else if ( LOCVOL[0](i) >= 0 ) { - LOCVOL[0](i) /= mean_plus; - LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f ); - } else { - LOCVOL[0](i) /= -mean_minus; - LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f ); + if ( LOCVOL[0](i) >= THRESHOLD ) { + auto tmp = LOCVOL[0](i)/ mean_plus; + LOCVOL[0](i) = std::min( tmp, 1.0f ); + } + else { + auto tmp = -LOCVOL[0](i)/mean_minus; + LOCVOL[0](i) = std::max( tmp, -1.0f ); } + //LOCVOL[0](i) = MASK(i); } - // Fill the source data for the coarse meshes + if (rank==0) printf("Coarsen the mesh for N_levels=%i \n",N_levels); + MPI_Barrier(comm); PROFILE_START("CoarsenMesh"); for (int i=1; i filter(ratio[0],ratio[1],ratio[2]); From 2f6b81ffeebea06ae8c3a558edeb42b41de5c29c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 20 Nov 2018 11:24:53 -0500 Subject: [PATCH 67/79] Added adaptive loop for morph_delta --- models/ColorModel.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 715cd8f2..1542a0d9 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -629,8 +629,14 @@ void ScaLBL_ColorModel::Run(){ double volB = Averages->Volume_w(); double volA = Averages->Volume_n(); double delta_volume = MorphInit(beta,morph_delta); + double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A + // update the volume volA += delta_volume; volB -= delta_volume; + //update size of morphological operation + morph_delta *= delta_volume_target / delta_volume; + if (morph_delta > 1.f) morph_delta = 1.f; + if (morph_delta < -1.f) morph_delta = -1.f; //MORPH_ADAPT = false; if (volB/(volA + volB) > TARGET_SATURATION){ MORPH_ADAPT = false; From 64edb9b18851899598653912c38a6a7d17d26498 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 20 Nov 2018 14:56:16 -0500 Subject: [PATCH 68/79] iterate on morphological change --- models/ColorModel.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 5da1a3b3..9255d69c 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -635,10 +635,12 @@ void ScaLBL_ColorModel::Run(){ // update the volume volA += delta_volume; volB -= delta_volume; - //update size of morphological operation - morph_delta *= delta_volume_target / delta_volume; + //update size of morphological operation -- don't change sign of morph_delta + if (delta_volume_target / delta_volume > 0.f) + morph_delta *= 0.5*delta_volume_target / delta_volume; if (morph_delta > 1.f) morph_delta = 1.f; if (morph_delta < -1.f) morph_delta = -1.f; + if (rank==0) printf(" Adjust morph delta: %f \n", morph_delta); //MORPH_ADAPT = false; if (volB/(volA + volB) > TARGET_SATURATION){ MORPH_ADAPT = false; From 1488a12b892e6897d2db4f29241ab19a8c6d86e8 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Tue, 20 Nov 2018 16:13:17 -0500 Subject: [PATCH 69/79] iterate on morphological change --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 9255d69c..e0a2c994 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -637,7 +637,7 @@ void ScaLBL_ColorModel::Run(){ volB -= delta_volume; //update size of morphological operation -- don't change sign of morph_delta if (delta_volume_target / delta_volume > 0.f) - morph_delta *= 0.5*delta_volume_target / delta_volume; + morph_delta *= 0.2*delta_volume_target / delta_volume; if (morph_delta > 1.f) morph_delta = 1.f; if (morph_delta < -1.f) morph_delta = -1.f; if (rank==0) printf(" Adjust morph delta: %f \n", morph_delta); From 4dedfca4b04483f5fa6ea2263593860693774ce0 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 21 Nov 2018 13:54:30 -0500 Subject: [PATCH 70/79] added minimum for morph_delta --- models/ColorModel.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index e0a2c994..8fdf8582 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -640,6 +640,7 @@ void ScaLBL_ColorModel::Run(){ morph_delta *= 0.2*delta_volume_target / delta_volume; if (morph_delta > 1.f) morph_delta = 1.f; if (morph_delta < -1.f) morph_delta = -1.f; + if (fabs(morph_delta) < 0.05 ) morph_delta = 0.05*(morph_delta)/fabs(morph_delta); // set minimum if (rank==0) printf(" Adjust morph delta: %f \n", morph_delta); //MORPH_ADAPT = false; if (volB/(volA + volB) > TARGET_SATURATION){ From 5f1eb0ec295d5068ddfcf447d593269fde387c47 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 21 Nov 2018 14:25:06 -0500 Subject: [PATCH 71/79] change morph_delta adapt --- models/ColorModel.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8fdf8582..017b7928 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -635,13 +635,13 @@ void ScaLBL_ColorModel::Run(){ // update the volume volA += delta_volume; volB -= delta_volume; - //update size of morphological operation -- don't change sign of morph_delta - if (delta_volume_target / delta_volume > 0.f) - morph_delta *= 0.2*delta_volume_target / delta_volume; - if (morph_delta > 1.f) morph_delta = 1.f; - if (morph_delta < -1.f) morph_delta = -1.f; - if (fabs(morph_delta) < 0.05 ) morph_delta = 0.05*(morph_delta)/fabs(morph_delta); // set minimum - if (rank==0) printf(" Adjust morph delta: %f \n", morph_delta); + if (delta_volume_target / delta_volume > 0.f){ + morph_delta *= 1.01*min((delta_volume_target - delta_volume) / delta_volume, 2.f); + if (morph_delta > 1.f) morph_delta = 1.f; + if (morph_delta < -1.f) morph_delta = -1.f; + if (fabs(morph_delta) < 0.05 ) morph_delta = 0.05*(morph_delta)/fabs(morph_delta); // set minimum + if (rank==0) printf(" Adjust morph delta: %f \n", morph_delta); + } //MORPH_ADAPT = false; if (volB/(volA + volB) > TARGET_SATURATION){ MORPH_ADAPT = false; From 4dd60d9d1e6912769294cb2ad77d9826aa8356e9 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 21 Nov 2018 14:31:43 -0500 Subject: [PATCH 72/79] change morph_delta adapt - fix bug --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 017b7928..8410b0d0 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -636,7 +636,7 @@ void ScaLBL_ColorModel::Run(){ volA += delta_volume; volB -= delta_volume; if (delta_volume_target / delta_volume > 0.f){ - morph_delta *= 1.01*min((delta_volume_target - delta_volume) / delta_volume, 2.f); + morph_delta *= 1.01*min((delta_volume_target - delta_volume) / delta_volume, 2.0); if (morph_delta > 1.f) morph_delta = 1.f; if (morph_delta < -1.f) morph_delta = -1.f; if (fabs(morph_delta) < 0.05 ) morph_delta = 0.05*(morph_delta)/fabs(morph_delta); // set minimum From ae292dc47676b0cd4ff3e9ddf3f46486428d0c67 Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 21 Nov 2018 15:28:53 -0500 Subject: [PATCH 73/79] change morph_delta adapt - adjust condition --- models/ColorModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8410b0d0..6ffebbbc 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -635,7 +635,7 @@ void ScaLBL_ColorModel::Run(){ // update the volume volA += delta_volume; volB -= delta_volume; - if (delta_volume_target / delta_volume > 0.f){ + if ((delta_volume_target - delta_volume) / delta_volume > 0.f){ morph_delta *= 1.01*min((delta_volume_target - delta_volume) / delta_volume, 2.0); if (morph_delta > 1.f) morph_delta = 1.f; if (morph_delta < -1.f) morph_delta = -1.f; From c5c99b52c185684ac35d3b9f29647470f21fcdea Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 21 Nov 2018 22:06:20 -0500 Subject: [PATCH 74/79] minor change in morphology switch --- models/ColorModel.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 6ffebbbc..431fdf8c 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -428,16 +428,17 @@ void ScaLBL_ColorModel::Run(){ } if (analysis_db->keyExists( "morph_delta" )){ morph_delta = analysis_db->getScalar( "morph_delta" ); - USE_MORPH = true; } else{ morph_delta=0.5; } if (analysis_db->keyExists( "morph_interval" )){ morph_interval = analysis_db->getScalar( "morph_interval" ); + USE_MORPH = true; } else{ - morph_interval=10000; + morph_interval=1000000; + USE_MORPH = false; } if (analysis_db->keyExists( "tolerance" )){ tolerance = analysis_db->getScalar( "tolerance" ); From 7feb729f1cf7ddf2c6d4e3749a98ab87deb4368e Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 24 Nov 2018 08:00:57 -0500 Subject: [PATCH 75/79] tweak relperm.csv --- models/ColorModel.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 431fdf8c..85d60891 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -577,9 +577,11 @@ void ScaLBL_ColorModel::Run(){ if (rank==0){ printf("** WRITE STEADY POINT *** "); printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); - + volA /= double(Nx*Ny*Nz*nprocs); + volB /= double(Nx*Ny*Nz*nprocs); FILE * kr_log_file = fopen("relperm.csv","a"); - fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",timestep-analysis_interval+20,muA,muB,5.796*alpha,Fx,Fy,Fz,volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g ",timestep-analysis_interval+20,muA,muB,5.796*alpha,Fx,Fy,Fz); + fprintf(kr_log_file,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z); fclose(kr_log_file); printf(" Measured capillary number %f \n ",Ca); From 3b6a7632ff132559265700550715b5e7e0f68428 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 24 Nov 2018 13:56:03 -0500 Subject: [PATCH 76/79] remove print statement from perm --- models/MRTModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 48060d96..6dd50167 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -204,7 +204,7 @@ void ScaLBL_MRTModel::Run(){ vay /= count; vaz /= count; - if (rank==0) printf("Computing Minkowski functionals \n"); + //if (rank==0) printf("Computing Minkowski functionals \n"); Morphology.ComputeScalar(Distance,0.f); Morphology.PrintAll(); From 61a580f44ddfbfc66a12cccd8dbfb6cf442eb310 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 24 Nov 2018 14:02:31 -0500 Subject: [PATCH 77/79] added permeability log file --- models/MRTModel.cpp | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 6dd50167..d0d61957 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -146,8 +146,12 @@ void ScaLBL_MRTModel::Run(){ Minkowski Morphology(Mask); int SIZE=Np*sizeof(double); - if (rank==0) printf("time Fx Fy Fz mu Vs As Js Xs vx vy vz\n"); - + if (rank==0){ + FILE * log_file = fopen("Permeability.csv","a"); + fprintf(log_file,"time Fx Fy Fz mu Vs As Js Xs vx vy vz\n"); + fclose(log_file); + } + //.......create and start timer............ double starttime,stoptime,cputime; ScaLBL_DeviceBarrier(); MPI_Barrier(comm); @@ -206,11 +210,14 @@ void ScaLBL_MRTModel::Run(){ //if (rank==0) printf("Computing Minkowski functionals \n"); Morphology.ComputeScalar(Distance,0.f); - Morphology.PrintAll(); - + //Morphology.PrintAll(); double mu = (tau-0.5)/3.f; - if (rank==0) printf("%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, - Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz); + if (rank==0) { + FILE * log_file = fopen("Permeability.csv","a"); + fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, + Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz); + fclose(log_file); + } } } //************************************************************************/ From 23d0e09956e2032d54be92a327c231f17814b055 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 24 Nov 2018 15:40:44 -0500 Subject: [PATCH 78/79] added absperm --- models/MRTModel.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index d0d61957..6bab022b 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -213,6 +213,9 @@ void ScaLBL_MRTModel::Run(){ //Morphology.PrintAll(); double mu = (tau-0.5)/3.f; if (rank==0) { + double h = Lz/double(Nz); + double absperm = h*h*mu*sqrt(vax*vax+vay*vay+vaz*vaz)/sqrt(Fx*Fx+Fy*Fy+Fz*Fz); + 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\n",timestep, Fx, Fy, Fz, mu, Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz); From 414abf0a55e303b7fac201fe13c80f566e4c6924 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sat, 24 Nov 2018 15:49:32 -0500 Subject: [PATCH 79/79] added absperm to CSV --- models/MRTModel.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 6bab022b..9e9afcca 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -148,7 +148,7 @@ void ScaLBL_MRTModel::Run(){ if (rank==0){ FILE * log_file = fopen("Permeability.csv","a"); - fprintf(log_file,"time Fx Fy Fz mu Vs As Js Xs vx vy vz\n"); + fprintf(log_file,"time Fx Fy Fz mu Vs As Js Xs vx vy vz k\n"); fclose(log_file); } @@ -217,8 +217,8 @@ void ScaLBL_MRTModel::Run(){ double absperm = h*h*mu*sqrt(vax*vax+vay*vay+vaz*vaz)/sqrt(Fx*Fx+Fy*Fy+Fz*Fz); 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\n",timestep, Fx, Fy, Fz, mu, - Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz); + fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, + Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz, absperm); fclose(log_file); } }