Merge pull request #124 from pynbody/gadget-hdf

Provide internal energy and smoothing output for GadgetHDF + baryon runs

genetIC/Dockerfile

@@ -1,8 +1,8 @@
-FROM ubuntu:20.04
+FROM ubuntu:24.04
 COPY ./ /genetIC
 
 RUN apt-get update && DEBIAN_FRONTEND=noninteractive apt-get install -y g++-12 libgsl-dev libfftw3-dev python3-pip pkg-config libhdf5-serial-dev
-RUN pip3 install numpy pynbody
+RUN pip3 install pynbody --break-system-packages # should use a venv but this is simpler for now
 RUN cd /genetIC && make clean && make
 
 ENTRYPOINT ["/genetIC/genetIC"]

genetIC/Makefile

@@ -24,7 +24,7 @@ CODEOPTIONS  = -DDOUBLEPRECISION -DOUTPUT_IN_DOUBLEPRECISION -DCUBIC_INTERPOLATI
 CPATH ?= /opt/local/include/
 LPATH ?= /opt/local/lib
 
-GSLFLAGS = -lgsl -lgslcblas
+GSLFLAGS = $(shell pkg-config --libs-only-L gsl) -lgsl -lgslcblas
 
 # Use -DFFTW3_THREADS=n to use threads. This will use n threads
 # UNLESS you are compiling with openmp, in which case you do not need to set n
@@ -36,13 +36,13 @@ GSLFLAGS = -lgsl -lgslcblas
 # Note that genetIC no longer supports the use of FFTW2
 
 FFTW = -DFFTW3 -DFFTW_THREADS
-FFTWLIB = -lfftw3 -lfftw3f -lfftw3_threads
+FFTWLIB = $(shell pkg-config --libs-only-L fftw3) -lfftw3 -lfftw3f -lfftw3_threads
 
 # if pkg-config is installed, add $(pkg-config --libs-only-L hdf5) to HDFLIB:
 HDFLIB = $(shell pkg-config --libs-only-L hdf5) -lhdf5
 
-# add also $(pkg-config --cflags hdf5) to CFLAGS:
-CFLAGS += $(shell pkg-config --cflags hdf5)
+# add hdf5, fftw3, and gsl include paths to CFLAGS
+CFLAGS += $(shell pkg-config --cflags hdf5 fftw3 gsl)
 
 
 # Provide information on the git repository, to be printed out on each run
@@ -93,7 +93,7 @@ ifeq ($(HOST), clema)
  endif
 
 ifeq ($(HOST), butte)
-	CXX = g++-13
+	CXX = g++-14
 	CPATH = /opt/homebrew/include/
 	LPATH = /opt/homebrew/lib/
 	CFLAGS += -Wl,-ld_classic # new macOS linker doesn't like C++, see XCode bug FB13097713

genetIC/src/cosmology/camb.hpp

@@ -65,6 +65,8 @@ namespace cosmology {
 
   public:
 
+    virtual ~PowerSpectrum() = default;
+
     //! \brief Evaluate power spectrum for a given species at wavenumber k (Mpc/h), including the normalisation
     //! transferType specifies the transfer function to use (currently dm or baryon)
     virtual CoordinateType operator()(CoordinateType k, particle::species transferType) const = 0;

genetIC/src/cosmology/parameters.hpp

@@ -28,8 +28,29 @@ namespace cosmology {
     FloatType sigma8; //!< Sigma8 parameter
     FloatType ns; //!< Scalar spectral index
     FloatType TCMB; //!< CMB temperature today, in K
+
+    FloatType smoothNorm; //!< Smoothing scale for the initial conditions, as multiple of grid spacing
   };
 
+  //! Estimate the temperature of the gas in the initial conditions, in K
+  template<typename FloatType>
+  FloatType getTemperature(const CosmologicalParameters<FloatType> &cosmo) {
+    if(cosmo.scalefactor > cosmo.scalefactorAtDecoupling)
+      return cosmo.TCMB * cosmo.scalefactorAtDecoupling / (cosmo.scalefactor * cosmo.scalefactor);
+    else
+      return cosmo.TCMB / cosmo.scalefactor;
+  }
+
+  //! Estimate the specific internal energy of gas in the initial conditions, in km^2 s^-2 a [gadget units]
+  template<typename FloatType>
+  FloatType getInternalEnergy(const CosmologicalParameters<FloatType> &cosmo) {
+    FloatType H_mass_fraction = 0.76; // Hydrogen mass fraction
+    FloatType mean_molecular_weight = 4.0 / (1.0 + 3.0 * H_mass_fraction);
+    FloatType boltzmann_constant = 1.380649e-29; // Boltzmann constant in km^2 kg s^-2 K^-1
+    FloatType proton_mass = 1.6726219e-27; // Proton mass in kg
+    return 1.5 * getTemperature(cosmo) * boltzmann_constant / ( proton_mass * mean_molecular_weight) / cosmo.scalefactor;
+  }
+
   //! Computes an estimate of the linear growth factor.
   template<typename FloatType>
   FloatType growthFactor(const CosmologicalParameters<FloatType> &cosmology) {

genetIC/src/ic.hpp

@@ -205,6 +205,7 @@ class ICGenerator {
     cosmology.TCMB = 2.725;
     cosmology.scalefactorAtDecoupling = 1./150.; // NB - this is NOT photon decoupling, it's baryon decoupling
     cosmology.OmegaM0 = 0.3;
+    setNumNeighbours(64);
     outputFolder = "./";
 
     // only an approximate value is needed to initialise a gas temperature
@@ -1259,6 +1260,16 @@ class ICGenerator {
     }
   }
 
+  //! \brief Set the number of neighbours to use for smoothing
+  /*!
+    Note that this is only used for output smoothing length estimates, not by any calculations within
+    genetIC itself.
+  */
+  void setNumNeighbours(unsigned int n) {
+    // density of particles is 1, so we want (4 pi smoothMultiple^3/3) = n
+    cosmology.smoothNorm = pow(3.0 * n / (4.0 * M_PI), 1.0 / 3.0);
+  }
+
 
   //! Outputs the ICs in the defined format, creating appropriate particle generators if required
   virtual void write() {

genetIC/src/io/gadgethdf.hpp

@@ -219,12 +219,23 @@ namespace io {
             }, "Masses");
         }
 
+        if (this->cosmology.OmegaBaryons0 > 0) {
+          OutputFloatType internalEnergy = cosmology::getInternalEnergy(this->cosmology);
+          saveBlock<OutputFloatType>(
+            particle::species::baryon,
+            [internalEnergy](auto &localIterator) {
+              return internalEnergy;
+            }, "InternalEnergy");
 
-
+          saveBlock<OutputFloatType>(
+            particle::species::baryon,
+            [](auto &localIterator) {
+              return localIterator.getSmoothingScale();
+            }, "SmoothingLength");
+        }
 
       }
 
-
     };
 
 

genetIC/src/io/swift.hpp

@@ -39,6 +39,10 @@ namespace io {
       HighFive::SilenceHDF5 silence; // suppresses HDF5 warnings while in scope
       SwiftHDFOutput<GridDataType, OutputFloatType> output(Boxlength, mapper, generators, cosmology, nFiles);
       output(name);
+      logging::entry() << "Swift output saved.";
+      logging::entry() << " Note that Swift output follows the GadgetHDF format and unit conventions.  You should";
+      logging::entry() << " ensure that cleanup_h_factors and cleanup_velocity_factors are both enabled within your";
+      logging::entry() << " Swift parameter file.";
 
     }
 

genetIC/src/io/tipsy.hpp

@@ -79,7 +79,7 @@ namespace io {
       //! \brief Initialise baryonic particles:
       template<typename T>
       void initialise(gas &p, const cosmology::CosmologicalParameters<T> &cosmo) {
-        p.temp = cosmo.TCMB * cosmo.scalefactorAtDecoupling / (cosmo.scalefactor * cosmo.scalefactor);
+        p.temp = cosmology::getTemperature(cosmo);
         p.metals = 0.0;
         p.rho = 0.0;
       }

genetIC/src/main.cpp

@@ -60,6 +60,8 @@ void setup_parser(tools::ClassDispatch<ICType, void> &dispatch) {
   dispatch.add_class_route("supersample_gas", &ICType::setSupersampleGas);
   dispatch.add_class_route("subsample", &ICType::setSubsample);
   dispatch.add_class_route("eps_norm", &ICType::setEpsNorm);
+  dispatch.add_class_route("num_neighbors", &ICType::setNumNeighbours);
+  dispatch.add_class_route("num_neighbours", &ICType::setNumNeighbours);
 
   // Grafic options
   dispatch.add_class_route("pvar", &ICType::setpvarValue);

genetIC/src/simulation/particles/generator.hpp

@@ -34,6 +34,9 @@ namespace particle {
     //! Gets the cell softening scale of the grid for this evaluator
     virtual T getEps() const = 0;
 
+    //! Gets the smoothing scale for this evaluator
+    virtual T getSmoothingScale() const = 0;
+
     //! Returns the particle at index id on the grid, without offset
     virtual Particle <T> getParticleNoOffset(size_t id) const = 0;
 

genetIC/src/simulation/particles/mapper/mapperiterator.hpp

@@ -243,12 +243,23 @@ namespace particle {
       }
 
 
-      //! Gets the mass in the cell currently pointed at.
-      T getMass() const {
+
+      template<typename Function>
+      inline auto getParticleProperty(Function func) const {
         EvaluatorPtrType pEval;
         size_t id;
         std::tie(pEval, id) = getParticleEvaluatorAndIndex();
-        return pEval->getMass();
+        return (pEval.get()->*func)();
+      }
+
+      //! Gets the mass in the cell currently pointed at.
+      T getMass() const {
+        return getParticleProperty(&ParticleEvaluator<GridDataType>::getMass);
+      }
+
+      //! Gets the smoothing scale in the cell currently pointed at.
+      T getSmoothingScale() const {
+        return getParticleProperty(&ParticleEvaluator<GridDataType>::getSmoothingScale);
       }
 
       //! Returns a pointer to the pair obtained by dereferencing the iterator at its current position

genetIC/src/simulation/particles/offsetgenerator.hpp

@@ -53,6 +53,10 @@ namespace particle {
     GridDataType getEps() const override {
       return underlying->getEps();
     }
+
+    GridDataType getSmoothingScale() const override {
+      return underlying->getSmoothingScale();
+    }
   };
 
   /*! \class OffsetMultiLevelParticleGenerator

genetIC/src/simulation/particles/zeldovich.hpp

@@ -118,6 +118,11 @@ namespace particle {
       return onGrid->cellSize * onGrid->cellSofteningScale * epsNorm;
     }
 
+    //! Gets an SPH smoothing scale estimate for a single particle
+    virtual T getSmoothingScale() const override {
+      return onGrid->cellSize * cosmology.smoothNorm;
+    }
+
   };
 
 

genetIC/tests/test_26a_gadgethdf_multifile_baryons/paramfile.txt

@@ -1,4 +1,4 @@
-# Based on test 10c, but use multi-file gadget output
+# Multi-file gadgethdf output with baryons
 
 
 # output parameters