update I/O format options to include hdf5

analysis/runAnalysis.cpp

@@ -607,6 +607,9 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> 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<bool>( "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<Database> 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<bool>( "format", "silo" );
+ 
+    IO::initialize( "", format, "false" );
     // Create the MeshDataStruct
     d_meshData.resize( 1 );
 

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;

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<int>( "analysis_interval" );
 	}
-	/*
-	int IMAGE_INDEX = 0;
-	int IMAGE_COUNT = 0;
-	std::vector<std::string> 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<double>( "vol_A_previous" );
-	}
-	if (color_db->keyExists( "log_krA_previous" )){
-		log_krA_prev  = color_db->getScalar<double>( "log_krA_previous" );
-	}
-	if (color_db->keyExists( "krA_morph_factor" )){
-		KRA_MORPH_FACTOR  = color_db->getScalar<double>( "krA_morph_factor" );
-	}
 
-	if (color_db->keyExists( "capillary_number" )){
-		capillary_number = color_db->getScalar<double>( "capillary_number" );
-		SET_CAPILLARY_NUMBER=true;
-	}
-	if (color_db->keyExists( "rescale_force_after_timestep" )){
-		RESCALE_FORCE_AFTER_TIMESTEP = color_db->getScalar<int>( "rescale_force_after_timestep" );
-		RESCALE_FORCE = true;
-	}
-	if (color_db->keyExists( "timestep" )){
-		timestep = color_db->getScalar<int>( "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<double>( "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<double>( "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<int>( "morph_interval" );
-		USE_MORPH = true;
-	}
-	if (analysis_db->keyExists( "use_morphopen_oil" )){
-		USE_MORPHOPEN_OIL = analysis_db->getScalar<bool>( "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<double>( "tolerance" );
-	}
-
-	if (analysis_db->keyExists( "min_steady_timesteps" )){
-		MIN_STEADY_TIMESTEPS = analysis_db->getScalar<int>( "min_steady_timesteps" );
-	}
-	if (analysis_db->keyExists( "max_steady_timesteps" )){
-		MAX_STEADY_TIMESTEPS = analysis_db->getScalar<int>( "max_steady_timesteps" );
-	}
-	if (analysis_db->keyExists( "max_morph_timesteps" )){
-		MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
-	}
-	
-	auto protocol = color_db->getWithDefault<std::string>( "protocol", "none" );
-	if (protocol == "image sequence"){
-		// Get the list of images
-		USE_DIRECT = true;
-		ImageList = color_db->getVector<std::string>( "image_sequence");
-		IMAGE_INDEX = color_db->getWithDefault<int>( "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);
-*/
 }
 

models/MRTModel.cpp

@@ -20,6 +20,7 @@ void ScaLBL_MRTModel::ReadParams(string filename){
 	db = std::make_shared<Database>( 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);
+    auto format = vis_db->getWithDefault<string>( "format", "silo" );
-	int SIZE=Np*sizeof(double);
-	ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
-	ScaLBL_DeviceBarrier(); comm.barrier();
-	ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
 
-	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; n<ScaLBL_Comm->LastExterior(); 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(); n<ScaLBL_Comm->LastInterior(); 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<IO::MeshDataStruct> visData;
 	fillHalo<double> 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<IO::Variable>();
 	auto SignDistVar = std::make_shared<IO::Variable>();
 
-	IO::initialize("","silo","false");
+	IO::initialize("",format,"false");
 	// Create the MeshDataStruct	
 	visData.resize(1);
 	visData[0].meshName = "domain";

models/MRTModel.h

@@ -49,6 +49,7 @@ public:
     std::shared_ptr<Database> db;
     std::shared_ptr<Database> domain_db;
     std::shared_ptr<Database> mrt_db;
+    std::shared_ptr<Database> vis_db;
 
     IntArray Map;
     DoubleArray Distance;