adapted Swift output to include masses, smoothing lengths, and internal energies as gas particle properties

Makefile

@@ -16,10 +16,10 @@ FF	= gfortran-11# needed only for panphasia, or ## ifort
 LINKER	= g++-11 ## ifort # need to link with ifort if using intel
 OPT     = -Wall -Wno-unknown-pragmas -O3 -mtune=native #-fsanitize=address -fno-omit-frame-pointer
 CFLAGS  =  
-LFLAGS  = -lgsl -lgslcblas #-fsanitize=thread  #-fsanitize=address -fno-omit-frame-pointer
+LFLAGS  = -lgsl -lgslcblas -ltirpc #-fsanitize=thread  #-fsanitize=address -fno-omit-frame-pointer
 FFLAGS  = -ffixed-line-length-132 -O3 -fimplicit-none -g #-fsanitize=address -fno-omit-frame-pointer## use for gfortran
 #FFLAGS = -extend_source -O3 -fimplicit-none -g ## use for ifort
-CPATHS  = -I. -I$(HOME)/local/include -I/opt/local/include -I/usr/local/include -I/usr/include/hdf5/serial
+CPATHS  = -I. -I$(HOME)/local/include -I/opt/local/include -I/usr/local/include -I/usr/include/hdf5/serial -I/usr/include/tirpc
 LPATHS  = -L$(HOME)/local/lib -L/opt/local/lib -L/usr/local/lib -L/usr/lib/x86_64-linux-gnu/hdf5/serial 
 
 

main.cc

@@ -631,6 +631,10 @@ nthreads = omp_get_max_threads();
 				std::cout << "=============================================================\n";
 				std::cout << "   COMPUTING BARYON DENSITY\n";
 				std::cout << "-------------------------------------------------------------\n";
+				if( outformat == "swift"){
+					LOGUSER("Writing baryon data as particle attributes...");
+					the_output_plugin->write_gas_properties(f);
+				}
 				LOGUSER("Computing baryon density...");
 				GenerateDensityHierarchy(	cf, the_transfer_function_plugin, baryon , rh_TF, rand, f, false, bbshift );
 				coarsen_density(rh_Poisson, f, bspectral_sampling);
@@ -900,6 +904,12 @@ nthreads = omp_get_max_threads();
 				the_output_plugin->write_dm_mass(f);
 			}
 
+			// //For our current use with Richings, this is unnecessary here. Only leaving it b/c we might need it everywhere we also write_dm_mass.
+			// if( do_baryons && outformat == "swift"){
+			// 		LOGUSER("Writing baryon data as particle attributes...");
+			// 		the_output_plugin->write_gas_properties(f);
+			// 	}
+
 			u1 = f;	u1.zero();
 
 			//... compute 1LPT term
@@ -1077,6 +1087,11 @@ nthreads = omp_get_max_threads();
 				the_output_plugin->write_dm_mass(f);
 				u1 = f;	u1.zero();
 
+				if( do_baryons && outformat == "swift"){
+					LOGUSER("Writing baryon data as particle attributes...");
+					the_output_plugin->write_gas_properties(f);
+				}
+
 				if(bdefd)
 					f2LPT=f;
 

output.hh

@@ -116,6 +116,9 @@ public:
 	
 	//! purely virtual prototype to write the baryon gravitational potential (from which displacements are computed in 1LPT)
 	virtual void write_gas_potential( const grid_hierarchy& gh )  = 0;
+
+	//! purely virtual prototype to write the baryon masses, smoothing lengths and internal energy
+	virtual void write_gas_properties( const grid_hierarchy& gh )  = 0;
 	
 	//! purely virtual prototype for all things to be done at the very end
 	virtual void finalize( void ) = 0;

plugins/output_arepo.cc

@@ -654,6 +654,11 @@ public:
   void write_gas_potential(const grid_hierarchy &gh) { /* skip */
   }
 
+  void write_gas_properties( const grid_hierarchy& gh )
+  {
+    //... we don't care about gas properties for Arepo  
+  }
+
   void finalize(void) {
     // generate and add contiguous IDs for each particle type we have written
     generateAndWriteIDs();

plugins/output_art.cc

@@ -820,6 +820,7 @@ public:
 
   void write_gas_density(const grid_hierarchy &gh) {}
   void write_gas_potential(const grid_hierarchy &gh) {}
+  void write_gas_properties( const grid_hierarchy& gh ) {}
 
   void finalize(void) {
     this->write_header_file();

plugins/output_cart.cc

@@ -795,6 +795,8 @@ class cart_output_plugin : public output_plugin
 
 		void write_gas_position( int coord, const grid_hierarchy& gh ) {}
 
+		void write_gas_properties( const grid_hierarchy& gh ){} //... we don't care about gas properties for art  
+
 		void write_gas_velocity( int coord, const grid_hierarchy& gh )
 		{
 			size_t nptot = gh.count_leaf_cells(gh.levelmin(), gh.levelmax());

plugins/output_enzo.cc

@@ -592,6 +592,8 @@ public:
 
   void write_gas_potential( const grid_hierarchy& gh )
   { }
+  void write_gas_properties( const grid_hierarchy& gh )
+  { }
 
 
   void write_gas_velocity( int coord, const grid_hierarchy& gh )

plugins/output_gadget2.cc

@@ -1226,6 +1226,11 @@ public:
   {
     //... we don't care about gas potential for Gadget  
   }
+
+  void write_gas_properties( const grid_hierarchy& gh )
+  {
+    //... we don't care about gas properties for Gadget  
+  }
   
   //... write data for gas -- don't do this
   void write_gas_velocity( int coord, const grid_hierarchy& gh )

plugins/output_gadget2_2comp.cc

@@ -1305,7 +1305,9 @@ public:
 	
 	void write_gas_potential( const grid_hierarchy& gh )
 	{ }
-	
+
+	void write_gas_properties( const grid_hierarchy& gh )
+  	{ }
 	
 	
 	//... write data for gas -- don't do this

plugins/output_gadget_tetmesh.cc

@@ -1562,6 +1562,9 @@ public:
 	
 	void write_gas_density( const grid_hierarchy& gh )
 	{ }
+
+    void write_gas_properties( const grid_hierarchy& gh )
+    { }
 	
 	void finalize( void )
 	{	

plugins/output_generic.cc

@@ -228,6 +228,9 @@ public:
 	
 	void write_gas_position( int coord, const grid_hierarchy& gh )
 	{	}
+
+	void write_gas_properties( const grid_hierarchy& gh )
+  	{   }
 	
 	void write_gas_density( const grid_hierarchy& gh )
 	{	

plugins/output_grafic2.cc

@@ -580,7 +580,10 @@ public:
 	
 	void write_gas_position( int coord, const grid_hierarchy& gh )
 	{	/* do nothing, not used... */ }
-	
+
+	void write_gas_properties( const grid_hierarchy& gh )
+  	{ }
+
 	void finalize( void )
 	{	}
 	

plugins/output_swift.cc

@@ -46,6 +46,7 @@ protected:
   // parameter file hints
   int pmgrid, gridboost;
   float softening, Tini;
+  double gamma, YHe, Tcmb0; //for gas properties
 
   using output_plugin::cf_;
 
@@ -497,6 +498,38 @@ protected:
     writeHDF5_b("Coordinates", coord, GAS_PARTTYPE, data); // write highres gas
   }
 
+  template <typename T> void __write_gas_properties(const grid_hierarchy &gh) {
+    countLeafCells(gh);
+
+    std::vector<T> masses(npfine);
+    std::vector<T> smoothing_lengths(npfine);
+    std::vector<T> internal_energies(npfine);
+
+    T gas_mass = omega_b * rhoCrit * pow(boxSize * posFac, 3.0) / pow(2, 3 * levelmax_);
+    T smoothing_length = boxSize / hubbleParam / pow(2, levelmax_);
+
+    // calculate internal energy for gas
+    double npol = (fabs(1.0 - gamma) > 1e-7) ? 1.0 / (gamma - 1.) : 1.0;
+    double astart = 1.0 / (1.0 + redshift);
+    double h2 = hubbleParam * hubbleParam;
+    double adec = 1.0 / (160.0 * pow(omega_b * h2 / 0.022, 2.0 / 5.0));
+
+    Tini = astart < adec ? Tcmb0 / astart : Tcmb0 / astart / astart * adec;
+
+    const double mu = (Tini > 1.e4) ? 4.0 / (8. - 5. * YHe) : 4.0 / (1. + 3. * (1. - YHe));
+    T internal_energy = 1.3806e-16 / 1.6726e-24 * Tini * npol / mu / UnitVelocity_in_cm_per_s / UnitVelocity_in_cm_per_s;
+
+    for (size_t i = 0; i < npfine; i++){
+      masses[i] = gas_mass;
+      smoothing_lengths[i] = smoothing_length;
+      internal_energies[i] = internal_energy;
+    }
+
+    writeHDF5_a("Masses", GAS_PARTTYPE, masses); // write gas masses and other data
+    writeHDF5_a("SmoothingLength", GAS_PARTTYPE, smoothing_lengths);
+    writeHDF5_a("InternalEnergy", GAS_PARTTYPE, internal_energies);
+  }
+
 public:
   swift_output_plugin(config_file &cf) : output_plugin(cf) {
     // ensure that everyone knows we want to do SPH, implies: bsph=1, bbshift=1, decic_baryons=1
@@ -561,11 +594,13 @@ public:
 
     // calculate Tini for gas
     hubbleParam = cf.getValue<double>("cosmology", "H0") / 100.0;
+    Tcmb0 = cf.getValueSafe<double>("cosmology", "Tcmb0", 2.7255); //from monofonIC Planck2018EE+BAO+SN
+    gamma = cf.getValueSafe<double>("cosmology", "gamma", 5.0 / 3.0);
+    YHe = cf.getValueSafe<double>("cosmology", "YHe", 0.245421); //from monofonIC Planck2018EE+BAO+SN
 
     double astart = 1.0 / (1.0 + redshift);
     double h2 = hubbleParam * hubbleParam;
     double adec = 1.0 / (160.0 * pow(omega_b * h2 / 0.022, 2.0 / 5.0));
-    double Tcmb0 = 2.726;
 
     rhoCrit *= h2;
     posFac /= hubbleParam;
@@ -668,6 +703,13 @@ public:
       __write_gas_position<double>(coord, gh);
   }
 
+  void write_gas_properties(const grid_hierarchy &gh) {
+    if (!doublePrec)
+      __write_gas_properties<float>(gh);
+    else
+      __write_gas_properties<double>(gh);
+  }
+
   void write_gas_density(const grid_hierarchy &gh) {
     // if only saving highres gas, then all gas cells have the same initial mass
     // do not write out densities as we write out displacements

plugins/output_tipsy.cc

@@ -947,6 +947,11 @@ public:
     {
         //... we don't care about baryon potential for TIPSY
     }
+
+	void write_gas_properties( const grid_hierarchy& gh )
+  	{
+    	//... we don't care about gas properties for TIPSY  
+  	}
     
     //... write data for gas -- don't do this
     void write_gas_velocity( int coord, const grid_hierarchy& gh )

plugins/output_tipsy_resample.cc

@@ -1157,6 +1157,11 @@ public:
         //... we don't care about baryon potential for TIPSY
     }
 
+    void write_gas_properties( const grid_hierarchy& gh )
+    {
+        //... we don't care about gas properties for TIPSY  
+    }
+
     //... write data for gas
     void write_gas_velocity (int coord, const grid_hierarchy & gh)
     {