diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 13840d1b..282f6c69 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -491,7 +491,7 @@ 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, +runAnalysis::runAnalysis(std::shared_ptr input_db, const RankInfoStruct& rank_info, std::shared_ptr ScaLBL_Comm, std::shared_ptr Dm, int Np, bool Regular, IntArray Map ): d_Np( Np ), d_regular ( Regular), @@ -501,6 +501,7 @@ runAnalysis::runAnalysis(std::shared_ptr db, const RankInfoStruct& ran d_ScaLBL_Comm( ScaLBL_Comm) { + auto db = input_db->getDatabase( "Analysis" ); // Ids of work items to use for dependencies ThreadPool::thread_id_t d_wait_blobID; ThreadPool::thread_id_t d_wait_analysis; @@ -729,11 +730,12 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep ) /****************************************************************** * Run the analysis * ******************************************************************/ -void runAnalysis::run( std::shared_ptr db, TwoPhase& Averages, const double *Phi, +void runAnalysis::run( std::shared_ptr input_db, TwoPhase& Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den) { int N = d_N[0]*d_N[1]*d_N[2]; + auto db = input_db->getDatabase( "Analysis" ); int timestep = db->getWithDefault( "timestep", 0 ); // Check which analysis steps we need to perform @@ -871,12 +873,13 @@ void runAnalysis::run( std::shared_ptr db, TwoPhase& Averages, const d // if ( matches(type,AnalysisType::CreateRestart) ) { if (timestep%d_restart_interval==0){ - if (d_rank==0) { - FILE *Rst = fopen("Restart.txt","w"); - fprintf(Rst,"%i\n",timestep+4); - fclose(Rst); - } - // Write the restart file (using a seperate thread) + if (d_rank==0) { + input_db->putScalar( "Restart", true ); + std::ofstream OutStream("Restart.db"); + input_db->print(OutStream, ""); + OutStream.close(); + } + // Write the restart file (using a seperate thread) auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cDen,cfq,d_Np); work->add_dependency(d_wait_restart); d_wait_restart = d_tpool.add_work(work); @@ -899,12 +902,14 @@ void runAnalysis::run( std::shared_ptr db, TwoPhase& Averages, const d /****************************************************************** * Run the analysis * ******************************************************************/ -void runAnalysis::basic( std::shared_ptr db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den) +void runAnalysis::basic( std::shared_ptr input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den) { int N = d_N[0]*d_N[1]*d_N[2]; // Check which analysis steps we need to perform - int timestep = db->getWithDefault( "timestep", 0 ); + auto color_db = input_db->getDatabase( "Color" ); + + int timestep = color_db->getWithDefault( "timestep", 0 ); auto type = computeAnalysisType( timestep ); if ( type == AnalysisType::AnalyzeNone ) return; @@ -973,9 +978,12 @@ void runAnalysis::basic( std::shared_ptr db, SubPhase &Averages, const ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double)); if (d_rank==0) { - FILE *Rst = fopen("Restart.txt","w"); - fprintf(Rst,"%i\n",timestep+4); - fclose(Rst); + color_db->putScalar( "Restart", true ); + input_db->putDatabase("Color", color_db); + std::ofstream OutStream("Restart.db"); + input_db->print(OutStream, ""); + OutStream.close(); + } // Write the restart file (using a seperate thread) auto work1 = new WriteRestartWorkItem(d_restartFile.c_str(),cDen,cfq,d_Np); diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index f8f316df..f3f6ad48 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -132,14 +132,38 @@ void ScaLBL_ColorModel::ReadParams(string filename){ inletB=0.f; outletA=0.f; outletB=1.f; + if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details) BoundaryCondition = 0; if (domain_db->keyExists( "BC" )){ BoundaryCondition = domain_db->getScalar( "BC" ); } - if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details) - + + // Override user-specified boundary condition for specific protocols + auto protocol = color_db->getWithDefault( "protocol", "none" ); + if (protocol == "seed water"){ + if (BoundaryCondition != 0 ){ + BoundaryCondition = 0; + if (rank==0) printf("WARNING: protocol (seed water) supports only full periodic boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + } + else if (protocol == "open connected oil"){ + if (BoundaryCondition != 0 ){ + BoundaryCondition = 0; + if (rank==0) printf("WARNING: protocol (open connected oil) supports only full periodic boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + } + else if (protocol == "shell aggregation"){ + if (BoundaryCondition != 0 ){ + BoundaryCondition = 0; + if (rank==0) printf("WARNING: protocol (shell aggregation) supports only full periodic boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + } } + void ScaLBL_ColorModel::SetDomain(){ Dm = std::shared_ptr(new Domain(domain_db,comm)); // full domain for analysis Mask = std::shared_ptr(new Domain(domain_db,comm)); // mask domain removes immobile phases @@ -255,7 +279,7 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) AFFINITY=AffinityList[idx]; label_count[idx] += 1.0; idx = NLABELS; - Mask->id[n] = 0; // set mask to zero since this is an immobile component + //Mask->id[n] = 0; // set mask to zero since this is an immobile component } } // fluid labels are reserved @@ -265,9 +289,10 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) } } } + // Set Dm to match Mask for (int i=0; iid[i] = Mask->id[i]; - + for (int idx=0; idxComm, label_count[idx]); if (rank==0){ @@ -387,17 +412,8 @@ void ScaLBL_ColorModel::Initialize(){ if (Restart == true){ if (rank==0){ printf("Reading restart file! \n"); - ifstream restart("Restart.txt"); - if (restart.is_open()){ - restart >> timestep; - printf("Restarting from timestep =%i \n",timestep); - } - else{ - printf("WARNING:No Restart.txt file, setting timestep=0 \n"); - timestep=0; - } } - MPI_Bcast(×tep,1,MPI_INT,0,comm); + // Read in the restart file to CPU buffers int *TmpMap; TmpMap = new int[Np]; @@ -504,13 +520,15 @@ void ScaLBL_ColorModel::Run(){ double delta_volume_target = 0.0; double RESIDUAL_ENDPOINT_THRESHOLD = 0.04; - auto protocol = color_db->getWithDefault( "protocol", "default" ); + auto protocol = color_db->getWithDefault( "protocol", "none" ); if (protocol == "image sequence"){ // Get the list of images USE_DIRECT = true; ImageList = color_db->getVector( "image_sequence"); IMAGE_INDEX = color_db->getWithDefault( "image_index", 0 ); IMAGE_COUNT = ImageList.size(); + morph_interval = 10000; + USE_MORPH = true; } else if (protocol == "seed water"){ morph_delta = 0.05; @@ -535,6 +553,9 @@ void ScaLBL_ColorModel::Run(){ capillary_number = color_db->getScalar( "capillary_number" ); SET_CAPILLARY_NUMBER=true; } + if (color_db->keyExists( "timestep" )){ + timestep = color_db->getScalar( "timestep" ); + } 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; @@ -621,7 +642,8 @@ void ScaLBL_ColorModel::Run(){ PROFILE_START("Loop"); //std::shared_ptr analysis_db; bool Regular = false; - runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); + auto current_db = db->cloneDatabase(); + 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)); } @@ -704,8 +726,9 @@ void ScaLBL_ColorModel::Run(){ // Run the analysis //analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den ); - analysis_db->putScalar("timestep",timestep); - analysis.basic( analysis_db, *Averages, Phi, Pressure, Velocity, fq, Den ); + color_db->putScalar("timestep",timestep); + current_db->putDatabase("Color", color_db); + analysis.basic( current_db, *Averages, Phi, Pressure, Velocity, fq, Den ); if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){ printf("....inlet pressure=%f \n",din); @@ -765,13 +788,45 @@ void ScaLBL_ColorModel::Run(){ 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 pB = Averages->gwb.p; - - double h = Dm->voxel_length; - double kAeff = h*h*muA*flow_rate_A/(rhoA*force_mag); - double kBeff = h*h*muB*flow_rate_B/(rhoB*force_mag); + 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); + // connected contribution + double Vol_nc = Averages->gnc.V/Dm->Volume; + double Vol_wc = Averages->gwc.V/Dm->Volume; + double Vol_nd = Averages->gnd.V/Dm->Volume; + double Vol_wd = Averages->gwd.V/Dm->Volume; + double Mass_n = Averages->gnc.M + Averages->gnd.M; + double Mass_w = Averages->gwc.M + Averages->gwd.M; + double vAc_x = Averages->gnc.Px/Mass_n; + double vAc_y = Averages->gnc.Py/Mass_n; + double vAc_z = Averages->gnc.Pz/Mass_n; + double vBc_x = Averages->gwc.Px/Mass_w; + double vBc_y = Averages->gwc.Py/Mass_w; + double vBc_z = Averages->gwc.Pz/Mass_w; + // disconnected contribution + double vAd_x = Averages->gnd.Px/Mass_n; + double vAd_y = Averages->gnd.Py/Mass_n; + double vAd_z = Averages->gnd.Pz/Mass_n; + double vBd_x = Averages->gwd.Px/Mass_w; + double vBd_y = Averages->gwd.Py/Mass_w; + double vBd_z = Averages->gwd.Pz/Mass_w; + + double flow_rate_A_connected = Vol_nc*(vAc_x*dir_x + vAc_y*dir_y + vAc_z*dir_z); + double flow_rate_B_connected = Vol_wc*(vBc_x*dir_x + vBc_y*dir_y + vBc_z*dir_z); + double flow_rate_A_disconnected = Vol_nd*(vAd_x*dir_x + vAd_y*dir_y + vAd_z*dir_z); + double flow_rate_B_disconnected = Vol_wd*(vBd_x*dir_x + vBd_y*dir_y + vBd_z*dir_z); + + double kAeff_connected = h*h*muA*flow_rate_A_connected/(rhoA*force_mag); + double kBeff_connected = h*h*muB*flow_rate_B_connected/(rhoB*force_mag); + + double kAeff = h*h*muA*(flow_rate_A_connected+flow_rate_A_disconnected)/(rhoA*force_mag); + double kBeff = h*h*muB*(flow_rate_B_connected+flow_rate_B_disconnected)/(rhoB*force_mag); double viscous_pressure_drop = (rhoA*volA + rhoB*volB)*force_mag; double Mobility = muA/muB; @@ -783,9 +838,9 @@ void ScaLBL_ColorModel::Run(){ WriteHeader=true; kr_log_file = fopen("relperm.csv","a"); if (WriteHeader) - fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water cap.pressure pressure.drop Ca M\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility); + fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected cap.pressure cap.pressure.connected pressure.drop Ca M\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility); - fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility); + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility); fclose(kr_log_file); printf(" Measured capillary number %f \n ",Ca); @@ -807,6 +862,7 @@ void ScaLBL_ColorModel::Run(){ } if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca); Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta); + color_db->putVector("F",{Fx,Fy,Fz}); } CURRENT_STEADY_TIMESTEPS = 0; } @@ -829,6 +885,7 @@ void ScaLBL_ColorModel::Run(){ if (IMAGE_INDEX < IMAGE_COUNT){ std::string next_image = ImageList[IMAGE_INDEX]; if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX); + color_db->putScalar("image_index",IMAGE_INDEX); ImageInit(next_image); } else{ @@ -916,46 +973,42 @@ double ScaLBL_ColorModel::ImageInit(std::string Filename){ if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str()); Mask->Decomp(Filename); for (int i=0; iid[i]; // save what was read - + for (int i=0; iid[i] = Mask->id[i]; // save what was read + double *PhaseLabel; PhaseLabel = new double[Nx*Ny*Nz]; AssignComponentLabels(PhaseLabel); - // consistency check double Count = 0.0; double PoreCount = 0.0; - for (int k=0; kSDs(i,j,k); - if (distance > 0.0){ - if (id[Nx*Ny*k+Nx*j+i] == 2){ - PoreCount++; - Count++; - } - else if (id[Nx*Ny*k+Nx*j+i] == 1){ - PoreCount++; - } -/* else if (suppress == false){ - printf("WARNING (ScaLBLColorModel::ImageInit) image input file sequence may not be labeled correctly (rank=%i) \n",rank); - suppress = true; - } - */ + for (int k=1; kComm, Count); PoreCount=sumReduce( Dm->Comm, PoreCount); - if (rank==0) printf(" new saturation: %f \n", Count / PoreCount); + if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount); ScaLBL_CopyToDevice(Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double)); MPI_Barrier(comm); + ScaLBL_D3Q19_Init(fq, Np); 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); MPI_Barrier(comm); + ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double)); + double saturation = Count/PoreCount; return saturation;