From 7442a317840f43fcf32d063a460fc4e1363965ff Mon Sep 17 00:00:00 2001 From: JamesEMcclure Date: Fri, 12 Jul 2019 09:21:16 -0400 Subject: [PATCH] write id in vis --- analysis/SubPhase.cpp | 95 ++++++++++++++++++++++++++ analysis/SubPhase.h | 3 +- analysis/runAnalysis.cpp | 4 +- models/ColorModel.cpp | 143 +++++++++++---------------------------- 4 files changed, 140 insertions(+), 105 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 9c2f4529..be07002b 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -624,4 +624,99 @@ void SubPhase::Full(){ } +void SubPhase::AggregateLabels(char *FILENAME){ + + int nx = Dm->Nx; + int ny = Dm->Ny; + int nz = Dm->Nz; + + int npx = Dm->nprocx(); + int npy = Dm->nprocy(); + int npz = Dm->nprocz(); + + int ipx = Dm->iproc(); + int ipy = Dm->jproc(); + int ipz = Dm->kproc(); + + int full_nx = npx*(nx-2); + int full_ny = npy*(ny-2); + int full_nz = npz*(nz-2); + int local_size = (nx-2)*(ny-2)*(nz-2); + long int full_size = long(full_nx)*long(full_ny)*long(full_nz); + + signed char *LocalID; + LocalID = new signed char [local_size]; + + //printf("aggregate labels: local size=%i, global size = %i",local_size, full_size); + // assign the ID for the local sub-region + for (int k=1; kid[n]; + if (local_id_val > 0){ + double value = Phi(i,j,k); + if (value > 0.0) local_id_val = 1; + else local_id_val = 2; + } + LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val; + } + } + } + MPI_Barrier(Dm->Comm); + + // populate the FullID + if (Dm->rank() == 0){ + signed char *FullID; + FullID = new signed char [full_size]; + // first handle local ID for rank 0 + for (int k=1; krank(); rnk++){ + ipz = rnk / (npx*npy); + ipy = (rnk - ipz*npx*npy) / npx; + ipx = (rnk - ipz*npx*npy - ipy*npx); + int tag = 15+rnk; + MPI_Recv(FullID,local_size,MPI_CHAR,rnk,tag,Dm->Comm,MPI_STATUS_IGNORE); + for (int k=1; krank(); + int dstrank = 0; + MPI_Send(LocalID,local_size,MPI_CHAR,dstrank,tag,Dm->Comm); + } + MPI_Barrier(Dm->Comm); + +} + + diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h index 738d2188..683fc46a 100644 --- a/analysis/SubPhase.h +++ b/analysis/SubPhase.h @@ -101,7 +101,8 @@ public: void Basic(); void Full(); void Write(int time); - + void AggregateLabels(char *FILENAME); + private: FILE *TIMELOG; FILE *SUBPHASE; diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 319014c5..d765eb25 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -232,6 +232,9 @@ public: fillData.copy(Averages.Vel_z,VelzData); fillData.copy(Averages.morph_n->label,BlobData); IO::writeData( timestep, visData, comm.comm ); + char CurrentIDFilename[40]; + sprintf(CurrentIDFilename,"id_t%d.raw",timestep); + Averages.AggregateLabels(CurrentIDFilename); PROFILE_STOP("Save Vis",1); }; private: @@ -1005,4 +1008,3 @@ void runAnalysis::WriteVisData( int timestep, SubPhase &Averages, const double * PROFILE_STOP("write vis"); } - diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index e2c44bc0..7c774957 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -633,6 +633,7 @@ void ScaLBL_ColorModel::Run(){ double volA = Averages->gnb.V; volA /= Dm->Volume; volB /= Dm->Volume;; + initial_volume = volA*Dm->Volume; double vA_x = Averages->gnb.Px/Averages->gnb.M; double vA_y = Averages->gnb.Py/Averages->gnb.M; double vA_z = Averages->gnb.Pz/Averages->gnb.M; @@ -666,10 +667,9 @@ void ScaLBL_ColorModel::Run(){ isSteady = true; if ( isSteady ){ - initial_volume = volA*Dm->Volume; MORPH_ADAPT = true; CURRENT_MORPH_TIMESTEPS=0; - delta_volume_target = Dm->Volume*morph_delta; //*volA ???? // set target volume change + delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change Averages->Full(); Averages->Write(timestep); analysis.WriteVisData( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); @@ -748,8 +748,7 @@ void ScaLBL_ColorModel::Run(){ else if (USE_SEED){ delta_volume = volA*Dm->Volume - initial_volume; CURRENT_MORPH_TIMESTEPS += analysis_interval; - double effective_seed_water = seed_water*(1.0+volB/volA); - double massChange = SeedPhaseField(effective_seed_water); + double massChange = SeedPhaseField(seed_water); if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", seed_water, delta_volume, delta_volume_target); } else if (USE_MORPHOPEN_OIL){ @@ -968,16 +967,13 @@ double ScaLBL_ColorModel::SeedPhaseField(const double seed_water_in_oil){ srand(time(NULL)); double mass_loss =0.f; double count =0.f; - double *Aq_tmp, *Bq_tmp, *Vel_tmp; - double SEED_THRESHOLD = -0.9; + double *Aq_tmp, *Bq_tmp; Aq_tmp = new double [7*Np]; Bq_tmp = new double [7*Np]; - Vel_tmp = new double [3*Np]; ScaLBL_CopyToHost(Aq_tmp, Aq, 7*Np*sizeof(double)); ScaLBL_CopyToHost(Bq_tmp, Bq, 7*Np*sizeof(double)); - ScaLBL_CopyToHost(Vel_tmp, Velocity, 3*Np*sizeof(double)); /* for (int k=1; kLastExterior(); n++){ - * ux = Vel[n]; - * uy = Vel[n+Np]; - * uz = Vel[n+2*Np]; - * magforce = sqrt(Fx*Fx+ Fy*Fy + Fz*Fz); - * dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce; - * if ( dotprod > flow_rate_A + flow_rate_B ){ - * - * } - * } - * - * for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){ - * } - * - */ - double count_voxels = 0.0; - double average_flow_rate=0.0; - double ux,uy,uz,dotprod; - double magforce = sqrt(Fx*Fx+ Fy*Fy + Fz*Fz); - for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){ - ux = Vel_tmp[n]; - uy = Vel_tmp[n+Np]; - uz = Vel_tmp[n+2*Np]; - dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce; - average_flow_rate += dotprod; - count_voxels += 1.0; - } - for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){ - ux = Vel_tmp[n]; - uy = Vel_tmp[n+Np]; - uz = Vel_tmp[n+2*Np]; - dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce; - average_flow_rate += dotprod; - count_voxels += 1.0; - } - count_voxels=sumReduce( Dm->Comm, count_voxels); - average_flow_rate=sumReduce( Dm->Comm, average_flow_rate); - average_flow_rate /= count_voxels; - - double oil_value = 0.0; - double water_value = 1.0; for (int n=0; n < ScaLBL_Comm->LastExterior(); n++){ + double random_value = seed_water_in_oil*double(rand())/ RAND_MAX; double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; double phase_id = (dA - dB) / (dA + dB); - double random_value = double(rand())/ RAND_MAX; - - // bias toward seed higher flow rate voxels - ux = Vel_tmp[n]; - uy = Vel_tmp[n+Np]; - uz = Vel_tmp[n+2*Np]; - dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce; - - if (dotprod > average_flow_rate && phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){ - Aq_tmp[n] = 0.3333333333333333*oil_value; - Aq_tmp[n+Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value; + if (phase_id > 0.0){ + Aq_tmp[n] -= 0.3333333333333333*random_value; + Aq_tmp[n+Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value; - Bq_tmp[n] = 0.3333333333333333*water_value; - Bq_tmp[n+Np] = 0.1111111111111111*water_value; - Bq_tmp[n+2*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+3*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+4*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+5*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+6*Np] = 0.1111111111111111*water_value; - mass_loss += oil_value - dA; - count++; + Bq_tmp[n] += 0.3333333333333333*random_value; + Bq_tmp[n+Np] += 0.1111111111111111*random_value; + Bq_tmp[n+2*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+3*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+4*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+5*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+6*Np] += 0.1111111111111111*random_value; } + mass_loss += random_value*seed_water_in_oil; } for (int n=ScaLBL_Comm->FirstInterior(); n < ScaLBL_Comm->LastInterior(); n++){ - double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; - double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; + double random_value = seed_water_in_oil*double(rand())/ RAND_MAX; + double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; + double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; double phase_id = (dA - dB) / (dA + dB); - double random_value = double(rand())/ RAND_MAX; - - // bias toward seed higher flow rate voxels - ux = Vel_tmp[n]; - uy = Vel_tmp[n+Np]; - uz = Vel_tmp[n+2*Np]; - dotprod = (Fx*ux + Fy*uy + Fz*uz) / magforce; - if (dotprod > average_flow_rate && phase_id > SEED_THRESHOLD && random_value < seed_water_in_oil){ - Aq_tmp[n] = 0.3333333333333333*oil_value; - Aq_tmp[n+Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+2*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+3*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+4*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+5*Np] = 0.1111111111111111*oil_value; - Aq_tmp[n+6*Np] = 0.1111111111111111*oil_value; + if (phase_id > 0.0){ + Aq_tmp[n] -= 0.3333333333333333*random_value; + Aq_tmp[n+Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value; - Bq_tmp[n] = 0.3333333333333333*water_value; - Bq_tmp[n+Np] = 0.1111111111111111*water_value; - Bq_tmp[n+2*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+3*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+4*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+5*Np] = 0.1111111111111111*water_value; - Bq_tmp[n+6*Np] = 0.1111111111111111*water_value; - mass_loss += oil_value - dA; - count++; + Bq_tmp[n] += 0.3333333333333333*random_value; + Bq_tmp[n+Np] += 0.1111111111111111*random_value; + Bq_tmp[n+2*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+3*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+4*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+5*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+6*Np] += 0.1111111111111111*random_value; } + mass_loss += random_value*seed_water_in_oil; } count= sumReduce( Dm->Comm, count);