diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 4c82b377..c7bad7f3 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -607,6 +607,9 @@ runAnalysis::runAnalysis( std::shared_ptr input_db, const RankInfoStru auto db = input_db->getDatabase( "Analysis" ); auto vis_db = input_db->getDatabase( "Visualization" ); + + /* set the I/O format */ + format = vis_db->getWithDefault( "format", "silo" ); // Ids of work items to use for dependencies ThreadPool::thread_id_t d_wait_blobID; @@ -645,8 +648,10 @@ runAnalysis::runAnalysis( std::shared_ptr input_db, const RankInfoStru d_rank = d_comm.getRank(); writeIDMap( ID_map_struct(), 0, id_map_filename ); + // Initialize IO for silo - IO::initialize( "", "silo", "false" ); + //std::string format = "silo"; + IO::initialize( "", format, "false" ); // Create the MeshDataStruct d_meshData.resize( 1 ); @@ -786,7 +791,11 @@ runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel) d_rank = d_comm.getRank(); writeIDMap( ID_map_struct(), 0, id_map_filename ); // Initialize IO for silo - IO::initialize( "", "silo", "false" ); + //std::string format = "silo"; + + format = vis_db->getWithDefault( "format", "silo" ); + + IO::initialize( "", format, "false" ); // Create the MeshDataStruct d_meshData.resize( 1 ); diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h index c7c4ce71..9600a297 100644 --- a/analysis/runAnalysis.h +++ b/analysis/runAnalysis.h @@ -100,6 +100,8 @@ private: int d_subphase_analysis_interval; double d_beta; bool d_regular; + std::string format; // IO format string "silo" or "hdf5" + ThreadPool d_tpool; RankInfoStruct d_rank_info; IntArray d_Map; diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index b94b4b1f..8f924aa2 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -754,7 +754,6 @@ double ScaLBL_ColorModel::Run(int returntime){ timestep = INITIAL_TIMESTEP; TRIGGER_FORCE_RESCALE = true; if (rank == 0) printf(" Capillary number missed target value = %f (measured value was Ca = %f) \n ",capillary_number, Ca); - } if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_TIMESTEP > RESCALE_FORCE_AFTER_TIMESTEP){ @@ -925,198 +924,7 @@ void ScaLBL_ColorModel::Run(){ if (analysis_db->keyExists( "analysis_interval" )){ analysis_interval = analysis_db->getScalar( "analysis_interval" ); } - /* - int IMAGE_INDEX = 0; - int IMAGE_COUNT = 0; - std::vector ImageList; - bool SET_CAPILLARY_NUMBER = false; - bool RESCALE_FORCE = false; - bool MORPH_ADAPT = false; - bool USE_MORPH = false; - bool USE_SEED = false; - bool USE_DIRECT = false; - bool USE_MORPHOPEN_OIL = false; - int MAX_MORPH_TIMESTEPS = 50000; // maximum number of LBM timesteps to spend in morphological adaptation routine - int MIN_STEADY_TIMESTEPS = 100000; - int MAX_STEADY_TIMESTEPS = 200000; - int RESCALE_FORCE_AFTER_TIMESTEP = 0; - int RAMP_TIMESTEPS = 0;//50000; // number of timesteps to run initially (to get a reasonable velocity field before other pieces kick in) - int CURRENT_MORPH_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time) - int CURRENT_STEADY_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time) - int morph_interval = 100000; - int morph_timesteps = 0; - double morph_delta = 0.0; - double seed_water = 0.0; - double capillary_number = 0.0; - double tolerance = 0.01; - double Ca_previous = 0.f; - double initial_volume = 0.0; - double delta_volume = 0.0; - double delta_volume_target = 0.0; - - double KRA_MORPH_FACTOR=0.5; - double volA_prev = 0.0; - double log_krA_prev = 1.0; - double log_krA_target = 1.0; - double log_krA = 1.0; - double slope_krA_volume = 0.0; - if (color_db->keyExists( "vol_A_previous" )){ - volA_prev = color_db->getScalar( "vol_A_previous" ); - } - if (color_db->keyExists( "log_krA_previous" )){ - log_krA_prev = color_db->getScalar( "log_krA_previous" ); - } - if (color_db->keyExists( "krA_morph_factor" )){ - KRA_MORPH_FACTOR = color_db->getScalar( "krA_morph_factor" ); - } - if (color_db->keyExists( "capillary_number" )){ - capillary_number = color_db->getScalar( "capillary_number" ); - SET_CAPILLARY_NUMBER=true; - } - if (color_db->keyExists( "rescale_force_after_timestep" )){ - RESCALE_FORCE_AFTER_TIMESTEP = color_db->getScalar( "rescale_force_after_timestep" ); - RESCALE_FORCE = true; - } - if (color_db->keyExists( "timestep" )){ - timestep = color_db->getScalar( "timestep" ); - } - if (BoundaryCondition != 0 && BoundaryCondition != 5 && SET_CAPILLARY_NUMBER==true){ - if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 or 5 \n"); - SET_CAPILLARY_NUMBER=false; - } - if (analysis_db->keyExists( "seed_water" )){ - seed_water = analysis_db->getScalar( "seed_water" ); - if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water); - } - if (analysis_db->keyExists( "morph_delta" )){ - morph_delta = analysis_db->getScalar( "morph_delta" ); - if (rank == 0) printf("Target volume change %f (morph_delta) \n",morph_delta); - } - if (analysis_db->keyExists( "morph_interval" )){ - 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" ); - } - - if (analysis_db->keyExists( "min_steady_timesteps" )){ - MIN_STEADY_TIMESTEPS = analysis_db->getScalar( "min_steady_timesteps" ); - } - if (analysis_db->keyExists( "max_steady_timesteps" )){ - MAX_STEADY_TIMESTEPS = analysis_db->getScalar( "max_steady_timesteps" ); - } - if (analysis_db->keyExists( "max_morph_timesteps" )){ - MAX_MORPH_TIMESTEPS = analysis_db->getScalar( "max_morph_timesteps" ); - } - - 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; - USE_SEED = false; - } - else if (protocol == "seed water"){ - morph_delta = -0.05; - seed_water = 0.01; - USE_SEED = true; - USE_MORPH = true; - } - else if (protocol == "open connected oil"){ - morph_delta = -0.05; - USE_SEED = false; - USE_MORPH = true; - USE_MORPHOPEN_OIL = true; - } - else if (protocol == "shell aggregation"){ - morph_delta = -0.05; - USE_MORPH = true; - USE_SEED = false; - } - else if (protocol == "fractional flow"){ - USE_MORPH = false; - USE_SEED = false; - } - else if (protocol == "centrifuge"){ - USE_MORPH = false; - USE_SEED = false; - } - else if (protocol == "core flooding"){ - USE_MORPH = false; - USE_SEED = false; - if (SET_CAPILLARY_NUMBER){ - double MuB = rhoB*(tauB - 0.5)/3.0; - double IFT = 6.0*alpha; - double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2)); - flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB; - } - } if (rank==0){ - printf("********************************************************\n"); - if (protocol == "image sequence"){ - printf(" using protocol = image sequence \n"); - printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); - printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); - printf(" tolerance = %f \n",tolerance); - std::string first_image = ImageList[IMAGE_INDEX]; - printf(" first image in sequence: %s ***\n", first_image.c_str()); - } - else if (protocol == "seed water"){ - printf(" using protocol = seed water \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(" 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); - } - else if (protocol == "fractional flow"){ - printf(" using protocol = fractional flow \n"); - printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); - printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); - printf(" tolerance = %f \n",tolerance); - } - else if (protocol == "centrifuge"){ - printf(" using protocol = centrifuge \n"); - printf(" driving force = %f \n",Fz); - if (Fz < 0){ - printf(" Component B displacing component A \n"); - } - else if (Fz > 0){ - printf(" Component A displacing component B \n"); - } - } - else if (protocol == "core flooding"){ - printf(" using protocol = core flooding \n"); - printf(" capillary number = %f \n", capillary_number); - } - printf("No. of timesteps: %i \n", timestepMax); - fflush(stdout); - } - */ //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); PROFILE_START("Loop"); @@ -1293,26 +1101,5 @@ void ScaLBL_ColorModel::WriteDebug(){ fwrite(PhaseField.data(),8,N,VELZ_FILE); fclose(VELZ_FILE); -/* ScaLBL_Comm->RegularLayout(Map,&ColorGrad[0],PhaseField); - FILE *CGX_FILE; - sprintf(LocalRankFilename,"Gradient_X.%05i.raw",rank); - CGX_FILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,CGX_FILE); - fclose(CGX_FILE); - - ScaLBL_Comm->RegularLayout(Map,&ColorGrad[Np],PhaseField); - FILE *CGY_FILE; - sprintf(LocalRankFilename,"Gradient_Y.%05i.raw",rank); - CGY_FILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,CGY_FILE); - fclose(CGY_FILE); - - ScaLBL_Comm->RegularLayout(Map,&ColorGrad[2*Np],PhaseField); - FILE *CGZ_FILE; - sprintf(LocalRankFilename,"Gradient_Z.%05i.raw",rank); - CGZ_FILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,CGZ_FILE); - fclose(CGZ_FILE); -*/ } diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index e1a451e2..2cbcd224 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -20,6 +20,7 @@ void ScaLBL_MRTModel::ReadParams(string filename){ db = std::make_shared( filename ); domain_db = db->getDatabase( "Domain" ); mrt_db = db->getDatabase( "MRT" ); + vis_db = db->getDatabase( "Visualization" ); tau = 1.0; timestepMax = 100000; @@ -371,51 +372,8 @@ void ScaLBL_MRTModel::Run(){ void ScaLBL_MRTModel::VelocityField(){ -/* Minkowski Morphology(Mask); - int SIZE=Np*sizeof(double); - ScaLBL_D3Q19_Momentum(fq,Velocity, Np); - ScaLBL_DeviceBarrier(); comm.barrier(); - ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE); + auto format = vis_db->getWithDefault( "format", "silo" ); - memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double)); - Morphology.Initialize(); - Morphology.UpdateMeshValues(); - Morphology.ComputeLocal(); - Morphology.Reduce(); - - double count_loc=0; - double count; - double vax,vay,vaz; - double vax_loc,vay_loc,vaz_loc; - vax_loc = vay_loc = vaz_loc = 0.f; - for (int n=0; nLastExterior(); n++){ - vax_loc += VELOCITY[n]; - vay_loc += VELOCITY[Np+n]; - vaz_loc += VELOCITY[2*Np+n]; - count_loc+=1.0; - } - - for (int n=ScaLBL_Comm->FirstInterior(); nLastInterior(); n++){ - vax_loc += VELOCITY[n]; - vay_loc += VELOCITY[Np+n]; - vaz_loc += VELOCITY[2*Np+n]; - count_loc+=1.0; - } - MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm); - MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm); - MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm); - MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm); - - vax /= count; - vay /= count; - vaz /= count; - - double mu = (tau-0.5)/3.f; - if (rank==0) printf("Fx Fy Fz mu Vs As Js Xs vx vy vz\n"); - if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Fx, Fy, Fz, mu, - Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz); - */ - std::vector visData; fillHalo fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1); @@ -424,7 +382,7 @@ void ScaLBL_MRTModel::VelocityField(){ auto VzVar = std::make_shared(); auto SignDistVar = std::make_shared(); - IO::initialize("","silo","false"); + IO::initialize("",format,"false"); // Create the MeshDataStruct visData.resize(1); visData[0].meshName = "domain"; diff --git a/models/MRTModel.h b/models/MRTModel.h index 4c41b746..f44ccf8a 100644 --- a/models/MRTModel.h +++ b/models/MRTModel.h @@ -49,6 +49,7 @@ public: std::shared_ptr db; std::shared_ptr domain_db; std::shared_ptr mrt_db; + std::shared_ptr vis_db; IntArray Map; DoubleArray Distance;