From 2e47e9c306ef361b841c5b31bbfb78083324fc79 Mon Sep 17 00:00:00 2001 From: James McClure Date: Sun, 7 Jun 2020 20:54:09 -0400 Subject: [PATCH 1/8] fix bug in Euler characteristic at domain boundary --- analysis/Minkowski.cpp | 6 ++-- models/ColorModel.cpp | 74 +++++++++++++++++++++--------------------- 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index faac6142..e98cfdc7 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -58,9 +58,9 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue) //int Nx = Field.size(0); //int Ny = Field.size(1); //int Nz = Field.size(2); - for (int k=1; k 0.01*volume_initial) - target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume); - delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor); + if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial); + double target_delta_volume_incremental = target_delta_volume; + if (fabs(target_delta_volume) > 0.01*volume_initial) + target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume); + delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor); - for (int k=0; kSDs(i,j,k) > 0.f){ - if (d < 3.f){ - //phase(i,j,k) = -1.0; - phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f); - } - } - } + for (int k=0; kSDs(i,j,k) > 0.f){ + if (d < 3.f){ + //phase(i,j,k) = -1.0; + phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f); + } + } + } + } + } + fillDouble.fill(phase); //} count = 0.f; From 5fae57157854561753c7290883052d6690e460b8 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 18 Jun 2020 19:56:38 -0400 Subject: [PATCH 2/8] fix connected pc --- models/ColorModel.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 2784e28a..583ca230 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -854,8 +854,8 @@ void ScaLBL_ColorModel::Run(){ double pB = Averages->gwb.p; double pAc = Averages->gnc.p; double pBc = Averages->gwc.p; - double pAB = (pA-pB)/(h*5.796*alpha); - double pAB_connected = (pAc-pBc)/(h*5.796*alpha); + double pAB = (pA-pB)/(h*6.0*alpha); + double pAB_connected = (pAc-pBc)/(h*6.0*alpha); // connected contribution double Vol_nc = Averages->gnc.V/Dm->Volume; double Vol_wc = Averages->gwc.V/Dm->Volume; From cc858ea87b9b777092af5e3e5e31ca9ce16d8846 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 19 Jun 2020 21:36:28 -0400 Subject: [PATCH 3/8] Revert "fix bug in Euler characteristic at domain boundary" This reverts commit 2e47e9c306ef361b841c5b31bbfb78083324fc79. --- analysis/Minkowski.cpp | 6 ++-- models/ColorModel.cpp | 74 +++++++++++++++++++++--------------------- 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index e98cfdc7..faac6142 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -58,9 +58,9 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue) //int Nx = Field.size(0); //int Ny = Field.size(1); //int Nz = Field.size(2); - for (int k=0; k 0.01*volume_initial) - target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume); - delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor); + if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial); + double target_delta_volume_incremental = target_delta_volume; + if (fabs(target_delta_volume) > 0.01*volume_initial) + target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume); + delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental, WallFactor); - for (int k=0; kSDs(i,j,k) > 0.f){ + if (d < 3.f){ + //phase(i,j,k) = -1.0; + phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f); + } + } + } } } - } - - CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance - - // 5. Update phase indicator field based on new distnace - for (int k=0; kSDs(i,j,k) > 0.f){ - if (d < 3.f){ - //phase(i,j,k) = -1.0; - phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f); - } - } - } - } - } - fillDouble.fill(phase); + fillDouble.fill(phase); //} count = 0.f; From 3ab182daf795b8fdcdbcb6088d3b134813c379f2 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 19 Jun 2020 21:39:54 -0400 Subject: [PATCH 4/8] fixed half edge rule for Euler characteristic --- analysis/Minkowski.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index faac6142..43ab6ee8 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -88,7 +88,7 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue) // printf(" (%i,%i,%i) RP(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e3,object.halfedge.twin(e3),P3.x,P3.y,P3.z,P1.x,P1.y,P1.z,a3,s3); //} // Euler characteristic (half edge rule: one face - 0.5*(three edges)) - Xi -= 0.5; + Xi -= 1.5; } // Euler characteristic -- each vertex shared by four cubes Xi += 0.25*double(object.VertexCount); From 435869740b87ef53c2e7e0a76f9a23d8fb44a293 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Mon, 22 Jun 2020 12:44:50 -0400 Subject: [PATCH 5/8] debug scalar MF --- analysis/Minkowski.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index 43ab6ee8..9652db35 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -58,9 +58,9 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue) //int Nx = Field.size(0); //int Ny = Field.size(1); //int Nz = Field.size(2); - for (int k=1; k Date: Wed, 24 Jun 2020 13:48:53 -0400 Subject: [PATCH 6/8] add Minkowski to TestTorus --- tests/TestTorus.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/TestTorus.cpp b/tests/TestTorus.cpp index 2d486774..d3cffa05 100644 --- a/tests/TestTorus.cpp +++ b/tests/TestTorus.cpp @@ -63,6 +63,8 @@ int main(int argc, char **argv) auto Dm = std::make_shared(db,comm); // const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); auto Averages = std::make_shared(Dm); + auto Object = std::make_shared(Dm); + Nx += 2; Ny += 2; Nz += 2; @@ -124,9 +126,11 @@ int main(int argc, char **argv) //Averages->Phase(i,j,k) = - Averages->Phase(i,j,k); if (Averages->Phase(i,j,k) > 0.0){ Dm->id[n] = 2; + Object->id(i,j,k) = 1; } else{ Dm->id[n] = 1; + Object->id(i,j,k) = 0; } Averages->SDn(i,j,k) = Averages->Phase(i,j,k); Averages->Phase(i,j,k) = Averages->SDn(i,j,k); @@ -163,6 +167,17 @@ int main(int argc, char **argv) Averages->Reduce(); Averages->PrintAll(int(5)); // Averages->Reduce(); + + Object->MeasureObject(); + double Vi = Object.V(); + double Ai = Object.A(); + double Hi = Object.H(); + double Xi = Object.X(); + Vi=sumReduce( Dm->Comm, Vi); + Ai=sumReduce( Dm->Comm, Ai); + Hi=sumReduce( Dm->Comm, Hi); + Xi=sumReduce( Dm->Comm, Xi); + printf("Vi=%f, Ai=%f, Hi=%f, Xi=%f \n", Vi,Ai,Hi,Xi); } // Limit scope so variables that contain communicators will free before MPI_Finialize MPI_Barrier(comm); From 1a5b46156dc1569024870eadf916bb783a08cc3e Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 24 Jun 2020 14:01:10 -0400 Subject: [PATCH 7/8] fix Torus --- analysis/Minkowski.cpp | 2 +- tests/TestTorus.cpp | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/analysis/Minkowski.cpp b/analysis/Minkowski.cpp index 9652db35..e98cfdc7 100644 --- a/analysis/Minkowski.cpp +++ b/analysis/Minkowski.cpp @@ -88,7 +88,7 @@ void Minkowski::ComputeScalar(const DoubleArray& Field, const double isovalue) // printf(" (%i,%i,%i) RP(%i,%i)={%f,%f,%f} {%f,%f,%f} a=%f l=%f \n",i,j,k,e3,object.halfedge.twin(e3),P3.x,P3.y,P3.z,P1.x,P1.y,P1.z,a3,s3); //} // Euler characteristic (half edge rule: one face - 0.5*(three edges)) - Xi -= 1.5; + Xi -= 0.5; } // Euler characteristic -- each vertex shared by four cubes Xi += 0.25*double(object.VertexCount); diff --git a/tests/TestTorus.cpp b/tests/TestTorus.cpp index d3cffa05..0cad888a 100644 --- a/tests/TestTorus.cpp +++ b/tests/TestTorus.cpp @@ -169,10 +169,10 @@ int main(int argc, char **argv) // Averages->Reduce(); Object->MeasureObject(); - double Vi = Object.V(); - double Ai = Object.A(); - double Hi = Object.H(); - double Xi = Object.X(); + double Vi = Object->V(); + double Ai = Object->A(); + double Hi = Object->H(); + double Xi = Object->X(); Vi=sumReduce( Dm->Comm, Vi); Ai=sumReduce( Dm->Comm, Ai); Hi=sumReduce( Dm->Comm, Hi); From 0d292a52b2b849afff065b689978aa72d2463cef Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Wed, 24 Jun 2020 23:33:36 -0400 Subject: [PATCH 8/8] fix volume bug --- common/Domain.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/common/Domain.cpp b/common/Domain.cpp index 3dec0128..eadc4592 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -232,7 +232,7 @@ void Domain::initialize( std::shared_ptr db ) if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0; // Fill remaining variables N = Nx*Ny*Nz; - Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0; + Volume = nx*ny*nz*nproc[0]*nproc[1]*nproc[2]*1.0; if (myrank==0) printf("voxel length = %f micron \n", voxel_length); @@ -266,8 +266,8 @@ void Domain::Decomp( const std::string& Filename ) outlet_layers_x = 0; outlet_layers_y = 0; outlet_layers_z = 0; - inlet_layers_phase=1; - outlet_layers_phase=2; + inlet_layers_phase=1; + outlet_layers_phase=2; checkerSize = 32; // Read domain parameters