diff --git a/analysis/GreyPhase.cpp b/analysis/GreyPhase.cpp index 6e7a958f..ce1e2547 100644 --- a/analysis/GreyPhase.cpp +++ b/analysis/GreyPhase.cpp @@ -1,13 +1,15 @@ -#include "analysis/SubPhase.h" +#include "analysis/GreyPhase.h" // Constructor -GreyPhase::GreyPhase(std::shared_ptr dm): +GreyPhaseAnalysis::GreyPhaseAnalysis(std::shared_ptr dm): Dm(dm) { Nx=dm->Nx; Ny=dm->Ny; Nz=dm->Nz; Volume=(Nx-2)*(Ny-2)*(Nz-2)*Dm->nprocx()*Dm->nprocy()*Dm->nprocz()*1.0; // Global arrays + SDs.resize(Nx,Ny,Nz); SDs.fill(0); + Porosity.resize(Nx,Ny,Nz); Porosity.fill(0); PhaseID.resize(Nx,Ny,Nz); PhaseID.fill(0); Rho_n.resize(Nx,Ny,Nz); Rho_n.fill(0); Rho_w.resize(Nx,Ny,Nz); Rho_w.fill(0); @@ -39,17 +41,17 @@ GreyPhase::GreyPhase(std::shared_ptr dm): // Destructor -GreyPhase::~GreyPhase() +GreyPhaseAnalysis::~GreyPhaseAnalysis() { } -void GreyPhase::Write(int timestep) +void GreyPhaseAnalysis::Write(int timestep) { } -void GreyPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double B) +void GreyPhaseAnalysis::SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double B) { Fx = force_x; Fy = force_y; @@ -62,7 +64,7 @@ void GreyPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, do beta = B; } -void GreyPhase::Basic(){ +void GreyPhaseAnalysis::Basic(){ int i,j,k,n,imin,jmin,kmin,kmax; // If external boundary conditions are set, do not average over the inlet @@ -72,22 +74,8 @@ void GreyPhase::Basic(){ double count_w = 0.0; double count_n = 0.0; - for (k=0; kid[n] > 0 ){ - // compute density - double nA = Rho_n(n); - double nB = Rho_w(n); - double phi = (nA-nB)/(nA+nB); - Phi(n) = phi; - } - } - } - } - + Water_local.reset(); + Oil_local.reset(); for (k=kmin; kid[n] > 0 ){ // compute density - double nA = Rho_n(n); - double nB = Rho_w(n); - double phi = (nA-nB)/(nA+nB); - /* if ( phi > 0.0 ){ - nA = 1.0; - nb.V += 1.0; - nb.M += nA*rho_n; - // velocity - nb.Px += rho_n*nA*Vel_x(n); - nb.Py += rho_n*nA*Vel_y(n); - nb.Pz += rho_n*nA*Vel_z(n); - } - else{ - nB = 1.0; - wb.M += nB*rho_w; - wb.V += 1.0; + double pressure = Pressure(n); + double rho_n = Rho_n(n); + double rho_w = Rho_w(n); + double porosity = Porosity(n); + double vel_x = Vel_x(n); + double vel_y = Vel_y(n); + double vel_z = Vel_z(n); + Water_local.p += pressure*(rho_w)/(rho_n+rho_w); + Water_local.M += rho_w*porosity; + Water_local.Px += rho_w*vel_x; + Water_local.Py += rho_w*vel_y; + Water_local.Pz += rho_w*vel_z; + Oil_local.p += pressure*(rho_n)/(rho_n+rho_w); + Oil_local.M += rho_n*porosity; + Oil_local.Px += rho_n*vel_x; + Oil_local.Py += rho_n*vel_y; + Oil_local.Pz += rho_n*vel_z; - // velocity - wb.Px += rho_w*nB*Vel_x(n); - wb.Py += rho_w*nB*Vel_y(n); - wb.Pz += rho_w*nB*Vel_z(n); - } - if ( phi > 0.99 ){ - nb.p += Pressure(n); - count_n += 1.0; - } - else if ( phi < -0.99 ){ - wb.p += Pressure(n); - count_w += 1.0; - } - */ } } } } - /* gwb.V=sumReduce( Dm->Comm, wb.V); - gnb.V=sumReduce( Dm->Comm, nb.V); - gwb.M=sumReduce( Dm->Comm, wb.M); - gnb.M=sumReduce( Dm->Comm, nb.M); - gwb.Px=sumReduce( Dm->Comm, wb.Px); - gwb.Py=sumReduce( Dm->Comm, wb.Py); - gwb.Pz=sumReduce( Dm->Comm, wb.Pz); - gnb.Px=sumReduce( Dm->Comm, nb.Px); - gnb.Py=sumReduce( Dm->Comm, nb.Py); - gnb.Pz=sumReduce( Dm->Comm, nb.Pz); - - count_w=sumReduce( Dm->Comm, count_w); - count_n=sumReduce( Dm->Comm, count_n); - if (count_w > 0.0) - gwb.p=sumReduce( Dm->Comm, wb.p) / count_w; - else - gwb.p = 0.0; - if (count_n > 0.0) - gnb.p=sumReduce( Dm->Comm, nb.p) / count_n; - else - gnb.p = 0.0; + Oil.M=sumReduce( Dm->Comm, Oil_local.M); + Oil.p=sumReduce( Dm->Comm, Oil_local.p); + Oil.Px=sumReduce( Dm->Comm, Oil_local.Px); + Oil.Py=sumReduce( Dm->Comm, Oil_local.Py); + Oil.Pz=sumReduce( Dm->Comm, Oil_local.Pz); + + Water.M=sumReduce( Dm->Comm, Water_local.M); + Water.p=sumReduce( Dm->Comm, Water_local.p); + Water.Px=sumReduce( Dm->Comm, Water_local.Px); + Water.Py=sumReduce( Dm->Comm, Water_local.Py); + Water.Pz=sumReduce( Dm->Comm, Water_local.Pz); + + Oil.p /= Oil.M; + Water.p /= Water.M; // check for NaN bool err=false; - if (gwb.V != gwb.V) err=true; - if (gnb.V != gnb.V) err=true; - if (gwb.p != gwb.p) err=true; - if (gnb.p != gnb.p) err=true; - if (gwb.Px != gwb.Px) err=true; - if (gwb.Py != gwb.Py) err=true; - if (gwb.Pz != gwb.Pz) err=true; - if (gnb.Px != gnb.Px) err=true; - if (gnb.Py != gnb.Py) err=true; - if (gnb.Pz != gnb.Pz) err=true; + if (Water.M != Water.M) err=true; + if (Water.p != Water.p) err=true; + if (Water.Px != Water.Px) err=true; + if (Water.Py != Water.Py) err=true; + if (Water.Pz != Water.Pz) err=true; + + if (Oil.M != Oil.M) err=true; + if (Oil.p != Oil.p) err=true; + if (Oil.Px != Oil.Px) err=true; + if (Oil.Py != Oil.Py) err=true; + if (Oil.Pz != Oil.Pz) err=true; if (Dm->rank() == 0){ double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz); @@ -194,23 +163,21 @@ void GreyPhase::Basic(){ dir_z = 1.0; force_mag = 1.0; } - double saturation=gwb.V/(gwb.V + gnb.V); - double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume; - double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume; - //double total_flow_rate = water_flow_rate + not_water_flow_rate; - //double fractional_flow = water_flow_rate / total_flow_rate; + saturation=Water.M/(Water.M + Oil.M); // assume constant density + water_flow_rate=saturation*(Water.Px*dir_x + Water.Py*dir_y + Water.Pz*dir_z)/Water.M; + oil_flow_rate=(1.0-saturation)*(Oil.Px*dir_x + Oil.Py*dir_y + Oil.Pz*dir_z)/Oil.M; double h = Dm->voxel_length; - double krn = h*h*nu_n*not_water_flow_rate / force_mag ; + double krn = h*h*nu_n*oil_flow_rate / force_mag ; double krw = h*h*nu_w*water_flow_rate / force_mag; //printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow); - fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, gwb.p, gnb.p); + fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*water_flow_rate,h*oil_flow_rate, Water.p, Oil.p); fflush(TIMELOG); } -*/ + if (err==true){ // exception if simulation produceds NaN - printf("GreyPhase.cpp: NaN encountered, may need to check simulation parameters \n"); + printf("GreyPhaseAnalysis.cpp: NaN encountered, may need to check simulation parameters \n"); } ASSERT(err==false); } diff --git a/analysis/GreyPhase.h b/analysis/GreyPhase.h index ddf43261..3019648e 100644 --- a/analysis/GreyPhase.h +++ b/analysis/GreyPhase.h @@ -6,7 +6,7 @@ #define GreyPhase_INC #include -#include "common/Domain.h" +#include "common/ScaLBL.h" #include "common/Communication.h" #include "analysis/analysis.h" #include "common/Utilities.h" @@ -16,6 +16,17 @@ #include "IO/Writer.h" class GreyPhase{ + public: + double p; + double M,Px,Py,Pz; + void reset(){ + p=M=Px=Py=Pz=0.0; + } + + private: +}; + +class GreyPhaseAnalysis{ public: std::shared_ptr Dm; double Volume; @@ -24,9 +35,16 @@ public: double nu_n, nu_w; double gamma_wn, beta; double Fx, Fy, Fz; + // outputs + double saturation,water_flow_rate, oil_flow_rate; + //simulation outputs (averaged values) + GreyPhase Water, Oil; + GreyPhase Water_local, Oil_local; //........................................................................... - int Nx,Ny,Nz; + int Nx,Ny,Nz; + IntArray SDs; // contains porosity map + IntArray Porosity; // contains porosity map IntArray PhaseID; // Phase ID array DoubleArray Rho_n; // density field DoubleArray Rho_w; // density field @@ -37,8 +55,8 @@ public: DoubleArray Vel_y; DoubleArray Vel_z; - GreyPhase(std::shared_ptr Dm); - ~GreyPhase(); + GreyPhaseAnalysis(std::shared_ptr Dm); + ~GreyPhaseAnalysis(); void SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha, double beta); void Basic(); diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp index e35aff4a..b83c0c53 100644 --- a/models/GreyscaleColorModel.cpp +++ b/models/GreyscaleColorModel.cpp @@ -135,7 +135,7 @@ void ScaLBL_GreyscaleColorModel::SetDomain(){ N = Nx*Ny*Nz; id = new signed char [N]; for (int i=0; iid[i] = 1; // initialize this way - Averages = std::shared_ptr ( new SubPhase(Dm) ); // TwoPhase analysis object + Averages = std::shared_ptr ( new GreyPhaseAnalysis(Dm) ); // TwoPhase analysis object MPI_Barrier(comm); Dm->CommInit(); MPI_Barrier(comm); @@ -831,15 +831,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ USE_SEED = true; USE_MORPH = true; } - else if (protocol == "open connected oil"){ - morph_delta = -0.05; - USE_MORPH = true; - USE_MORPHOPEN_OIL = true; - } - else if (protocol == "shell aggregation"){ - morph_delta = -0.05; - USE_MORPH = true; - } + if (greyscaleColor_db->keyExists( "capillary_number" )){ capillary_number = greyscaleColor_db->getScalar( "capillary_number" ); SET_CAPILLARY_NUMBER=true; @@ -868,11 +860,6 @@ void ScaLBL_GreyscaleColorModel::Run(){ morph_interval = analysis_db->getScalar( "morph_interval" ); USE_MORPH = true; } - if (analysis_db->keyExists( "use_morphopen_oil" )){ - USE_MORPHOPEN_OIL = analysis_db->getScalar( "use_morphopen_oil" ); - if (rank == 0 && USE_MORPHOPEN_OIL) printf("Volume change by morphological opening \n"); - USE_MORPH = true; - } if (analysis_db->keyExists( "tolerance" )){ tolerance = analysis_db->getScalar( "tolerance" ); } @@ -908,20 +895,6 @@ void ScaLBL_GreyscaleColorModel::Run(){ printf(" morph_delta = %f \n",morph_delta); printf(" seed_water = %f \n",seed_water); } - else if (protocol == "open connected oil"){ - printf(" using protocol = open connected oil \n"); - printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); - printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); - printf(" tolerance = %f \n",tolerance); - printf(" morph_delta = %f \n",morph_delta); - } - else if (protocol == "shell aggregation"){ - printf(" using protocol = shell aggregation \n"); - printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); - printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); - printf(" tolerance = %f \n",tolerance); - printf(" morph_delta = %f \n",morph_delta); - } printf("No. of timesteps: %i \n", timestepMax); fflush(stdout); } @@ -938,7 +911,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ //std::shared_ptr analysis_db; bool Regular = false; auto current_db = db->cloneDatabase(); - runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); + //runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); //analysis.createThreads( analysis_method, 4 ); while (timestep < timestepMax ) { //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); } @@ -1076,26 +1049,20 @@ void ScaLBL_GreyscaleColorModel::Run(){ if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition == 4){ printf("%i %f \n",timestep,din); + ScaLBL_Comm->RegularLayout(Map,Pressure,Averages->Pressure); + ScaLBL_Comm->RegularLayout(Map,&Den[0],Averages->Rho_n); + ScaLBL_Comm->RegularLayout(Map,&Den[Np],Averages->Rho_w); + ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Averages->Vel_x); + ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Averages->Vel_y); + ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Averages->Vel_z); + Averages->Basic(); } - // Run the analysis - analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den ); // allow initial ramp-up to get closer to steady state if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){ - analysis.finish(); + //analysis.finish(); CURRENT_STEADY_TIMESTEPS += analysis_interval; - double volB = Averages->gwb.V; - 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; - double vB_x = Averages->gwb.Px/Averages->gwb.M; - double vB_y = Averages->gwb.Py/Averages->gwb.M; - double vB_z = Averages->gwb.Pz/Averages->gwb.M; double muA = rhoA*(tauA-0.5)/3.f; double muB = rhoB*(tauB-0.5)/3.f; double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz); @@ -1109,10 +1076,12 @@ void ScaLBL_GreyscaleColorModel::Run(){ dir_z = 1.0; force_mag = 1.0; } - double current_saturation = volB/(volA+volB); - double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z); - double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z); - double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha); + double current_saturation = Averages->saturation; + double volA = current_saturation*Mask->Porosity()*Mask->Volume; + double volB = (1.0-current_saturation)*Mask->Porosity()*Mask->Volume; + double flow_rate_A = Averages->oil_flow_rate; + double flow_rate_B = Averages->water_flow_rate; + double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(6.0*alpha); if ( morph_timesteps > morph_interval ){ @@ -1165,17 +1134,17 @@ void ScaLBL_GreyscaleColorModel::Run(){ volA_prev = volA; //******************************** **/ /** compute averages & write data **/ - Averages->Full(); + /*Averages->Full(); Averages->Write(timestep); analysis.WriteVisData(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den ); analysis.finish(); - + */ if (rank==0){ printf("** WRITE STEADY POINT *** "); printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); double h = Dm->voxel_length; // pressures - double pA = Averages->gnb.p; + /*double pA = Averages->gnb.p; double pB = Averages->gwb.p; double pAc = Averages->gnc.p; double pBc = Averages->gwc.p; @@ -1231,7 +1200,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility); fclose(kr_log_file); - + */ printf(" Measured capillary number %f \n ",Ca); } if (SET_CAPILLARY_NUMBER ){ @@ -1283,16 +1252,6 @@ void ScaLBL_GreyscaleColorModel::Run(){ double massChange = SeedPhaseField(seed_water); if (rank==0) printf("***Seed water in oil %f, volume change %f / %f ***\n", massChange, delta_volume, delta_volume_target); } - else if (USE_MORPHOPEN_OIL){ - delta_volume = volA*Dm->Volume - initial_volume; - if (rank==0) printf("***Morphological opening of connected oil, target volume change %f ***\n", delta_volume_target); - MorphOpenConnected(delta_volume_target); - } - else { - if (rank==0) printf("***Shell aggregation, target volume change %f ***\n", delta_volume_target); - //double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A - delta_volume += MorphInit(beta,delta_volume_target-delta_volume); - } if ( (delta_volume - delta_volume_target)/delta_volume_target > 0.0 ){ MORPH_ADAPT = false; @@ -1316,7 +1275,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ } MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL); } - analysis.finish(); + //analysis.finish(); PROFILE_STOP("Loop"); PROFILE_SAVE("lbpm_color_simulator",1); //************************************************************************ @@ -1380,149 +1339,6 @@ double ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){ return saturation; } - -double ScaLBL_GreyscaleColorModel::MorphOpenConnected(double target_volume_change){ - - int nx = Nx; - int ny = Ny; - int nz = Nz; - int n; - int N = nx*ny*nz; - double volume_change=0.0; - - if (target_volume_change < 0.0){ - Array id_solid(nx,ny,nz); - Array phase_label(nx,ny,nz); - DoubleArray distance(Nx,Ny,Nz); - DoubleArray phase(nx,ny,nz); - signed char *id_connected; - id_connected = new signed char [nx*ny*nz]; - - ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); - - // Extract only the connected part of NWP - BlobIDstruct new_index; - double vF=0.0; double vS=0.0; - ComputeGlobalBlobIDs(nx-2,ny-2,nz-2,Dm->rank_info,phase,Averages->SDs,vF,vS,phase_label,Dm->Comm); - MPI_Barrier(Dm->Comm); - - long long count_connected=0; - long long count_porespace=0; - long long count_water=0; - for (int k=1; k 0){ - count_porespace++; - } - if (id[n] == 2){ - count_water++; - } - } - } - } - count_connected=sumReduce( Dm->Comm, count_connected); - count_porespace=sumReduce( Dm->Comm, count_porespace); - count_water=sumReduce( Dm->Comm, count_water); - - for (int k=0; kSDs(i,j,k) > 0.f){ - if (d < 3.f){ - phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f); - } - } - } - } - } - - int count_morphopen=0.0; - for (int k=1; kComm, count_morphopen); - volume_change = double(count_morphopen - count_connected); - - if (rank==0) printf(" opening of connected oil %f \n",volume_change/count_connected); - - ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double)); - ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){ - if (Dm->kproc()==0){ - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1); - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2); - } - if (Dm->kproc() == nprocz-1){ - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2); - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); - } - } - } - return(volume_change); -} double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil){ srand(time(NULL)); double mass_loss =0.f; @@ -1598,219 +1414,6 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil return(mass_loss); } -double ScaLBL_GreyscaleColorModel::MorphInit(const double beta, const double target_delta_volume){ - const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); - - double vF = 0.f; - double vS = 0.f; - double delta_volume; - double WallFactor = 0.0; - bool USE_CONNECTED_NWP = false; - - DoubleArray phase(Nx,Ny,Nz); - IntArray phase_label(Nx,Ny,Nz);; - DoubleArray phase_distance(Nx,Ny,Nz); - Array phase_id(Nx,Ny,Nz); - fillHalo fillDouble(Dm->Comm,Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1); - - - // Basic algorithm to - // 1. Copy phase field to CPU - ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double)); - - double count = 0.f; - for (int k=1; k 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f; - } - } - } - double volume_initial = sumReduce( Dm->Comm, count); - /* - sprintf(LocalRankFilename,"phi_initial.%05i.raw",rank); - FILE *INPUT = fopen(LocalRankFilename,"wb"); - fwrite(phase.data(),8,N,INPUT); - fclose(INPUT); - */ - // 2. Identify connected components of phase field -> phase_label - - double volume_connected = 0.0; - double second_biggest = 0.0; - if (USE_CONNECTED_NWP){ - BlobIDstruct new_index; - ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm); - MPI_Barrier(Dm->Comm); - - // only operate on component "0" - count = 0.0; - - for (int k=0; kComm, count); - second_biggest = sumReduce( Dm->Comm, second_biggest); - } - else { - // use the whole NWP - for (int k=0; kSDs(i,j,k) > 0.f){ - if (phase(i,j,k) > 0.f ){ - phase_id(i,j,k) = 0; - } - else { - phase_id(i,j,k) = 1; - } - } - else { - phase_id(i,j,k) = 1; - } - } - } - } - } - - /*int reach_x, reach_y, reach_z; - 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; - } - // erase the original object - phase(i,j,k) = -1.0; - } - } - } - } - - if (USE_CONNECTED_NWP){ - if (volume_connected - second_biggest < 2.0*fabs(target_delta_volume) && target_delta_volume < 0.0){ - // if connected volume is less than 2% just delete the whole thing - if (rank==0) printf("Connected region has shrunk! \n"); - REVERSE_FLOW_DIRECTION = true; - } - -/* else{*/ - if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest ); - } - 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); - } - } - } - } - } - fillDouble.fill(phase); - //} - - count = 0.f; - for (int k=1; k 0.f && Averages->SDs(i,j,k) > 0.f){ - count+=1.f; - } - } - } - } - double volume_final= sumReduce( Dm->Comm, count); - - delta_volume = (volume_final-volume_initial); - if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial); - if (rank == 0) printf(" new saturation = %f \n", volume_final/(0.238323*double((Nx-2)*(Ny-2)*(Nz-2)*nprocs))); - - // 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)); - /* - sprintf(LocalRankFilename,"dist_final.%05i.raw",rank); - FILE *DIST = fopen(LocalRankFilename,"wb"); - fwrite(phase_distance.data(),8,N,DIST); - fclose(DIST); - - sprintf(LocalRankFilename,"phi_final.%05i.raw",rank); - FILE *PHI = fopen(LocalRankFilename,"wb"); - fwrite(phase.data(),8,N,PHI); - fclose(PHI); - */ - // 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 == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){ - if (Dm->kproc()==0){ - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1); - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2); - } - if (Dm->kproc() == nprocz-1){ - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2); - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); - } - } - return delta_volume; -} - void ScaLBL_GreyscaleColorModel::WriteDebug(){ // Copy back final phase indicator field and convert to regular layout DoubleArray PhaseField(Nx,Ny,Nz); diff --git a/models/GreyscaleColorModel.h b/models/GreyscaleColorModel.h index b9e8d11f..d7298c10 100644 --- a/models/GreyscaleColorModel.h +++ b/models/GreyscaleColorModel.h @@ -10,8 +10,7 @@ Implementation of two-fluid greyscale color lattice boltzmann model #include #include "common/Communication.h" -#include "analysis/TwoPhase.h" -#include "analysis/runAnalysis.h" +#include "analysis/GreyPhase.h" #include "common/MPI_Helpers.h" #include "ProfilerApp.h" #include "threadpool/thread_pool.h" @@ -50,7 +49,7 @@ public: std::shared_ptr Mask; // this domain is for lbm std::shared_ptr ScaLBL_Comm; std::shared_ptr ScaLBL_Comm_Regular; - std::shared_ptr Averages; + std::shared_ptr Averages; // input database std::shared_ptr db;