Merge pull request #123 from pynbody/gadget-hdf

GadgetHDF support

.github/workflows/build-test.yaml

@@ -33,30 +33,33 @@ jobs:
       with:
         python-version: ${{ matrix.python-version }}
     - uses: actions/checkout@v2
+      with:
+        submodules: recursive
     - name: Install dependencies
       if: matrix.os == 'ubuntu-latest'
       shell: bash
       run: |
         sudo add-apt-repository ppa:ubuntu-toolchain-r/test
         sudo apt-get update -qq
-        sudo apt install libfftw3-dev libgsl0-dev
+        sudo apt install libfftw3-dev libgsl0-dev libhdf5-serial-dev
         sudo apt install gcc-10 g++-10 gcc-11 g++-11 gcc-12 g++-12
     - name: Install dependencies
       if: matrix.os == 'macos-latest'
       shell: bash
       run: |
-        brew install gcc@12 fftw gsl
+        brew install gcc@12 fftw gsl hdf5 
         brew link gsl
     - name: Compile code
       working-directory: genetIC
       run: make
     - name: Install python dependencies
       shell: bash
       run: |
-        export CC=gcc-10
+        export CC=gcc-12
+        export CXX=g++-12
         python -m pip install --upgrade pip setuptools wheel
         python -m pip install numpy scipy cython
-        python -m pip install pynbody
+        python -m pip install --pre pynbody
     - name: Run tests
       shell: bash
       working-directory: genetIC/tests

.gitmodules

@@ -0,0 +1,4 @@
+[submodule "genetIC/HighFive"]
+	path = genetIC/HighFive
+	url = https://github.com/BlueBrain/HighFive
+	branch = v2.9.0

CMakeLists.txt

@@ -2,11 +2,12 @@ cmake_minimum_required(VERSION 3.2)
 
 project(IC)
 
+
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/genetIC)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14  -g ")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++17 -g ")
 
-set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -DDEBUG_INFO")
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -DOPENMP -fopenmp -DDEBUG_INFO")
 set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fopenmp -O3 -DOPENMP")
 
 set(SOURCE_FILES
@@ -20,6 +21,7 @@ set(SOURCE_FILES
         genetIC/src/io.hpp
         genetIC/src/io/grafic.hpp
         genetIC/src/io/gadget.hpp
+        genetIC/src/io/gadgethdf.hpp
         genetIC/src/io/input.hpp
         genetIC/src/io/tipsy.hpp
         genetIC/src/main.cpp
@@ -61,12 +63,17 @@ set(SOURCE_FILES
         genetIC/src/simulation/modifications/quadraticmodification.hpp
         genetIC/src/simulation/multilevelgrid/mask.hpp
         genetIC/src/tools/memmap.hpp
-        genetIC/src/tools/numerics/tricubic.hpp genetIC/src/tools/logging.hpp genetIC/src/tools/logging.cpp genetIC/src/simulation/modifications/splice.hpp genetIC/src/tools/lru_cache.hpp)
+        genetIC/src/tools/numerics/tricubic.hpp genetIC/src/tools/logging.hpp genetIC/src/tools/logging.cpp genetIC/src/simulation/modifications/splice.hpp genetIC/src/tools/lru_cache.hpp
+        genetIC/src/io/swift.hpp)
 
+include_directories(/opt/homebrew/include)
+link_directories(/opt/homebrew/lib)
 include_directories( /opt/local/include )
 link_directories(/opt/local/lib )
+
 include_directories( genetIC/  )
-link_libraries(fftw3 m fftw3f fftw3_threads gsl gslcblas)
+include_directories( genetIC/HighFive/include )
+link_libraries(fftw3 m fftw3f fftw3_threads gsl gslcblas hdf5 )
 
 
 exec_program(
@@ -116,4 +123,4 @@ add_definitions(-DCUBIC_INTERPOLATION
         -DZELDOVICH_GRADIENT_FOURIER_SPACE
         )
 add_compile_options(-Wextra)
-add_executable(genetIC ${SOURCE_FILES})
+add_executable(genetIC ${SOURCE_FILES})

README.md

@@ -6,9 +6,10 @@ GenetIC
 GenetIC is a code to generate initial conditions for cosmological simulations, especially for zoom simulations of galaxies. It provides support for 'genetic modifications' as described by e.g. [Roth et al 2015](https://arxiv.org/abs/1504.07250), [Rey & Pontzen 2018](https://arxiv.org/abs/1706.04615). It also supports 'splicing' as described by [Cadiou et al 2021](https://arxiv.org/abs/2107.03407). 
 * The main code is written in C++ and can be found (with build instructions and a test suite) in the subfolder `genetIC`. 
 * The code is also available as [docker](http://docker.io) image `apontzen/genetic` (see the user manual for further information)
+* A **hands-on tutorial** for genetic modification can be found at https://github.com/pynbody/genetic_example
+* **User documentation** can be found in the form of a PDF manual, also at https://github.com/pynbody/genetIC/releases/latest
 * A set of python tools for development and testing is provided in `tools`.
 * If you wish to use the code for science, we recommend adopting a particular (numbered) version. The most recent release can be downloaded at https://github.com/pynbody/genetIC/releases/latest
-* Documentation can be found in the form of a PDF user manual, also at https://github.com/pynbody/genetIC/releases/latest
 * The technical underpinnings of the code are described in a release paper at https://arxiv.org/abs/2006.01841
 
 License and acknowledgement policy

genetIC/Dockerfile

@@ -1,7 +1,7 @@
 FROM ubuntu:20.04
 COPY ./ /genetIC
 
-RUN apt-get update && apt-get install -y g++-9 libgsl-dev libfftw3-dev python3-pip
+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 cd /genetIC && make clean && make
 

genetIC/Makefile

@@ -1,5 +1,5 @@
 # Reasonably sensible defaults
-CFLAGS ?= -Wall -g -lpthread -O3 -fopenmp -std=c++14 -fdiagnostics-color=auto -I`pwd`
+CFLAGS ?= -Wall -g -lpthread -O3 -fopenmp -std=c++17 -fdiagnostics-color=auto -I`pwd` -I`pwd`/HighFive/include
 
 # Here are code options. Probably don't fiddle with these unless you know what you're doing.
 #
@@ -38,6 +38,13 @@ GSLFLAGS = -lgsl -lgslcblas
 FFTW = -DFFTW3 -DFFTW_THREADS
 FFTWLIB = -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)
+
+
 # Provide information on the git repository, to be printed out on each run
 GIT_MODIFIED = $(shell git ls-files -m 2>/dev/null | tr "\n" " " )
 GIT_VERSION = $(shell git rev-parse --short HEAD 2>/dev/null)
@@ -58,14 +65,14 @@ ifeq ($(USE_CUFFT), 1)
 	FFTW = -DFFTW3
 	FFTWLIB = -lcufft -lcufftw
 	OLD_CFLAGS=$(CFLAGS)
-	CFLAGS=-O3 -I`pwd` -Xcompiler="-O3 -fopenmp -std=c++14 -Wall"
+	CFLAGS=-O3 -I`pwd` -I`pwd`/HighFive/include -Xcompiler="-O3 -fopenmp -std=c++14 -Wall"
 	CXX=nvcc
 endif
 ifeq ($(HOST3), pfe)
 	CXX = /nasa/pkgsrc/2016Q2/gcc5/bin/g++
 	CPATH = /u/apontzen/genetIC/genetIC:/nasa/gsl/1.14/include/:/nasa/intel/Compiler/2016.2.181/compilers_and_libraries_2016.2.181/linux/mkl/include/fftw
 	LPATH =  /nasa/gsl/1.14/lib/ -L/nasa/intel/Compiler/2016.2.181/compilers_and_libraries_2016.2.181/linux/mkl/lib/intel64/
-	CFLAGS = -Wall -g -O3 -fopenmp -std=c++14 -fdiagnostics-color=auto -I`pwd` -DOPENMP -DDOUBLEPRECISION
+	CFLAGS = -Wall -g -O3 -fopenmp -std=c++17 -fdiagnostics-color=auto -I`pwd` -I`pwd`/HighFive/include -DOPENMP -DDOUBLEPRECISION
 	FFTWLIB = -ldl -lpthread -lm -lstdc++ -lgfortran -lmkl_intel_ilp64 -lmkl_gnu_thread -lmkl_core
 	FFTW = -DFFTW3 -DFFTW_THREADS
 	GSLFLAGS = /nasa/gsl/1.14/lib/libgsl.a /nasa/gsl/1.14/lib/libgslcblas.a
@@ -74,20 +81,27 @@ endif
 ifeq ($(HOST), Snowd)
 	CPATH = /opt/local/include/
 	LPATH = /opt/local/lib/
-	CFLAGS = -O3 -fopenmp -std=c++14 -fdiagnostics-color=auto -I`pwd` -DOPENMP -DDOUBLEPRECISION
+	CFLAGS = -O3 -fopenmp -std=c++17 -fdiagnostics-color=auto -I`pwd` -I`pwd`/HighFive/include -DOPENMP -DDOUBLEPRECISION
 	CXX = /opt/local/bin/clang++-mp-6.0
 endif
 
 ifeq ($(HOST), clema)
         CPATH = /opt/local/include/
         LPATH = /opt/local/lib/
-	      CFLAGS = -O2 -g -fopenmp -std=c++14 -fdiagnostics-color=auto -I`pwd` -DOPENMP -DDOUBLEPRECISION
+	      CFLAGS = -O2 -g -fopenmp -std=c++17 -fdiagnostics-color=auto -I`pwd` -I`pwd`/HighFive/include -DOPENMP -DDOUBLEPRECISION
         CXX = /opt/local/bin/clang++-mp-6.0
  endif
 
+ifeq ($(HOST), butte)
+	CXX = g++-13
+	CPATH = /opt/homebrew/include/
+	LPATH = /opt/homebrew/lib/
+	CFLAGS += -Wl,-ld_classic # new macOS linker doesn't like C++, see XCode bug FB13097713
+endif
+
 ifeq ($(HOST), jasmi)
-	CXX = /opt/local/bin/g++-mp-11   -L/opt/local/lib/
-	CFLAGS = -O3 -g -fopenmp -std=c++14 -I`pwd` -DOPENMP -DDOUBLEPRECISION
+	CXX = /opt/local/bin/g++-mp-13 -L/opt/local/lib/
+	CFLAGS = -O3 -g -fopenmp -std=c++17 -I`pwd` -I`pwd`/HighFive/include -DOPENMP -DDOUBLEPRECISION
 endif
 
 ifeq ($(HOST), Rhodo)
@@ -134,7 +148,7 @@ all: genetIC
 %.o: %.cpp ; $(CXX) $(CFLAGS) $(CODEOPTIONS) $(GIT_VARIABLES) -I$(CPATH) $(FFTW) -c $< -o $@
 
 genetIC: src/main.o src/tools/filesystem.o src/tools/progress/progress.o src/tools/logging.o
-		$(CXX) $(CFLAGS) -o genetIC $(GIT_VARIABLES) -I$(CPATH) $(FFTW) src/main.o src/tools/filesystem.o src/tools/progress/progress.o src/tools/logging.o -L$(LPATH) $(GSLFLAGS) -lm $(FFTWLIB)
+		$(CXX) $(CFLAGS) -o genetIC $(GIT_VARIABLES) -I$(CPATH) $(FFTW) src/main.o src/tools/filesystem.o src/tools/progress/progress.o src/tools/logging.o -L$(LPATH) $(GSLFLAGS) -lm $(FFTWLIB) $(HDFLIB)
 
 clean:
 	rm -f genetIC

genetIC/src/ic.hpp

@@ -1285,6 +1285,15 @@ class ICGenerator {
                             cosmology, static_cast<int>(outputFormat),
                             nGadgetFiles);
         break;
+      case OutputFormat::gadgethdf:
+        gadgethdf::save<float>(getOutputPath(), boxlen, *pMapper, pParticleGenerator, cosmology,
+                               nGadgetFiles);
+        break;
+      case OutputFormat::swift:
+        swift::save<float>(getOutputPath(), boxlen, *pMapper, pParticleGenerator, cosmology,
+                               nGadgetFiles);
+        break;
+
       case OutputFormat::tipsy:
         tipsy::save(getOutputPath() + ".tipsy", boxlen, pParticleGenerator,
                     pMapper, cosmology);

genetIC/src/io.hpp

@@ -15,6 +15,8 @@ namespace cosmology {
 
 #include "io/input.hpp"
 #include "io/gadget.hpp"
+#include "io/gadgethdf.hpp"
+#include "io/swift.hpp"
 #include "io/tipsy.hpp"
 #include "io/grafic.hpp"
 
@@ -31,7 +33,7 @@ namespace cosmology {
 namespace io {
 
   enum class OutputFormat {
-    unknown = 1, gadget2 = 2, gadget3 = 3, tipsy = 4, grafic = 5
+    unknown = 1, gadget2 = 2, gadget3 = 3, tipsy = 4, grafic = 5, gadgethdf = 6, swift = 7
   };
 
   std::ostream &operator<<(std::ostream &outputStream, const OutputFormat &format) {
@@ -50,6 +52,13 @@ namespace io {
         break;
       case OutputFormat::grafic:
         outputStream << "grafic";
+        break;
+      case OutputFormat::gadgethdf:
+        outputStream << "gadgethdf";
+        break;
+      case OutputFormat::swift:
+        outputStream << "swift";
+        break;
     }
     return outputStream;
   }
@@ -60,7 +69,7 @@ namespace io {
     try {
       int i = std::stoi(s);
       format = static_cast<OutputFormat>(i);
-    } catch (std::invalid_argument e) {
+    } catch (const std::invalid_argument & e) {
       if (s == "gadget2") {
         format = OutputFormat::gadget2;
       } else if (s == "gadget3") {
@@ -69,8 +78,12 @@ namespace io {
         format = OutputFormat::tipsy;
       } else if (s == "grafic") {
         format = OutputFormat::grafic;
+      } else if (s == "gadgethdf") {
+        format = OutputFormat::gadgethdf;
+      } else if (s == "swift") {
+        format = OutputFormat::swift;
       } else {
-        inputStream.setstate(std::ios::failbit);
+          inputStream.setstate(std::ios::failbit);
       }
     }
     return inputStream;

genetIC/src/io/gadget.hpp

@@ -273,7 +273,7 @@ namespace io {
       void saveGadgetBlock(particle::species forSpecies,
                            std::function<WriteType(const particle::mapper::MapperIterator<GridDataType> &)> getData) {
 
-        size_t current_n = 0;
+
         size_t nTotalForThisBlock = 0;
 
         std::vector<decltype(writers[0].template getMemMapFortran<WriteType>(std::declval<size_t>()))> currentWriteBlocks;
@@ -301,37 +301,48 @@ namespace io {
         for(int i=0; i<nFiles; i++)
           currentWriteBlocks.push_back(writers[i].template getMemMapFortran<WriteType>(nPerFile[i]));
 
+        assert(this->template genericSaveBlock(particleTypes, getData, currentWriteBlocks) == nTotalForThisBlock);
+
+      }
+
+      template<typename TargetType, typename WriteType>
+      size_t genericSaveBlock(const std::vector<unsigned int> & particleTypes,
+                            std::function<WriteType(const particle::mapper::MapperIterator<GridDataType> &)> getData,
+                            std::vector<TargetType> & currentWriteBlocks) {
+        size_t current_n = 0;
+
         for (auto particle_type: particleTypes) {
           auto begin = mapper.beginParticleType(*generators[gadgetTypeToSpecies[particle_type]], particle_type);
           auto end = mapper.endParticleType(*generators[gadgetTypeToSpecies[particle_type]], particle_type);
           size_t nMax = end.getIndex() - begin.getIndex();
 
           current_n += begin.parallelIterate(
-            [&](size_t n_offset, const particle::mapper::MapperIterator<GridDataType> &localIterator) {
-              int fileNum = 0;
-              size_t addr = n_offset + current_n;
-
-              // Now figure out which file to dump this into.
-              // The following approach looks slow (doing it for every particle!),
-              // but code profiling suggests it's not a significant overhead.
-              // Easier to do this than to try and have thread-local variables
-              while(addr >= nPartPerFile[fileNum]) {
-                addr-=nPartPerFile[fileNum];
-                fileNum++;
-              }
-
-              currentWriteBlocks[fileNum][addr] = getData(localIterator);
-            }, nMax);
-
+              [&](size_t n_offset, const particle::mapper::MapperIterator<GridDataType> &localIterator) {
+                int fileNum = 0;
+                size_t addr = n_offset + current_n;
+
+                // Now figure out which file to dump this into.
+                // The following approach looks slow (doing it for every particle!),
+                // but code profiling suggests it's not a significant overhead.
+                // Easier to do this than to try and have thread-local variables
+                while (addr >= nPartPerFile[fileNum]) {
+                  addr -= nPartPerFile[fileNum];
+                  fileNum++;
+                }
+
+                currentWriteBlocks[fileNum][addr] = getData(localIterator);
+              }, nMax);
         }
-
-        assert(current_n == nTotalForThisBlock);
+        return current_n;
 
       }
 
+      //! \brief Returns false in base class. Override in child class to force variable mass output
+      virtual bool forceVariableMass() const { return false; }
+
       //! \brief Scan through data to obtain information on particles and their masses.
-      void preScanForMassesAndParticleNumbers() {
-        variableMass = false;
+      virtual void preScanForMassesAndParticleNumbers() {
+        variableMass = this->forceVariableMass();
         masses = vector<InternalFloatType>(6, 0.0);
         nPartPerType = vector<size_t>(6, 0);
         nTotal = 0;
@@ -401,8 +412,20 @@ namespace io {
 
       }
 
+      virtual void writeHeaderOneFile(size_t fileNumber, std::vector<size_t> nPartPerTypeThisFile) {
+        if (gadgetVersion == 3) {
+          writers[fileNumber].writeFortran(createGadget3Header<OutputFloatType>(masses, nPartPerTypeThisFile, nPartPerType, nFiles,
+                                                                       boxLength, cosmology));
+        } else if (gadgetVersion == 2) {
+          writers[fileNumber].writeFortran(createGadget2Header<OutputFloatType>(masses, nPartPerTypeThisFile, nPartPerType, nFiles,
+                                                                       boxLength, cosmology));
+        } else {
+          throw std::runtime_error("Unknown gadget format");
+        }
+      }
+
       //! \brief Output the gadget3 or gadget2 header:
-      void writeHeader() {
+      virtual void writeHeader() {
         std::vector<size_t> nPartPerTypeThisFile(6, 0);
         std::vector<size_t> nPartRemainingPerType = nPartPerType;
         size_t offset = 0;
@@ -433,15 +456,7 @@ namespace io {
 
           assert(nPartRemainingInFile == 0); // ensure we have assigned all the particles to a type
 
-          if (gadgetVersion == 3) {
-            writers[i].writeFortran(createGadget3Header<OutputFloatType>(masses, nPartPerTypeThisFile, nPartPerType, nFiles,
-                                                                         boxLength, cosmology));
-          } else if (gadgetVersion == 2) {
-            writers[i].writeFortran(createGadget2Header<OutputFloatType>(masses, nPartPerTypeThisFile, nPartPerType, nFiles,
-                                                                         boxLength, cosmology));
-          } else {
-            throw std::runtime_error("Unknown gadget format");
-          }
+          writeHeaderOneFile(i, nPartPerTypeThisFile);
 
           this->nPartPerTypePerFile.push_back(nPartPerTypeThisFile);
         }
@@ -481,7 +496,7 @@ namespace io {
       }
 
       //! \brief Operation to save gadget particles
-      void operator()(const std::string &name) {
+      virtual void operator()(const std::string &name) {
 
         preScanForMassesAndParticleNumbers();