add genericio output plugin, fix memory leak in fftw plans, modernize CMakeLists

CMakeLists.txt

@@ -196,66 +196,72 @@ set_target_properties(${PRGNAME} PROPERTIES CXX_STANDARD 14)
 if(MPI_CXX_FOUND)
   if(CODE_PRECISION STREQUAL "FLOAT")
     if(FFTW3_SINGLE_MPI_FOUND)
-      target_link_libraries(${PRGNAME} ${FFTW3_SINGLE_MPI_LIBRARY})
+      target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_SINGLE_MPI)
-      target_include_directories(${PRGNAME} PRIVATE ${FFTW3_INCLUDE_DIR_PARALLEL})
+      target_compile_definitions(${PRGNAME} PRIVATE "USE_FFTW_MPI")
-      target_compile_options(${PRGNAME} PRIVATE "-DUSE_FFTW_MPI")
     else()
       message(SEND_ERROR "MPI enabled but FFTW3 library not found with MPI support for single precision!")
     endif()
   elseif(CODE_PRECISION STREQUAL "DOUBLE")
     if(FFTW3_DOUBLE_MPI_FOUND)
-      target_link_libraries(${PRGNAME} ${FFTW3_DOUBLE_MPI_LIBRARY})
+      target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_DOUBLE_MPI)
-      target_include_directories(${PRGNAME} PRIVATE ${FFTW3_INCLUDE_DIR_PARALLEL})
+      target_compile_definitions(${PRGNAME} PRIVATE "USE_FFTW_MPI")
-      target_compile_options(${PRGNAME} PRIVATE "-DUSE_FFTW_MPI")
     else()
       message(SEND_ERROR "MPI enabled but FFTW3 library not found with MPI support for double precision!")
     endif()
   elseif(CODE_PRECISION STREQUAL "LONGDOUBLE")
     if(FFTW3_LONGDOUBLE_MPI_FOUND)
-      target_link_libraries(${PRGNAME} ${FFTW3_LONGDOUBLE_MPI_LIBRARY})
+      target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_LONGDOUBLE_MPI)
-      target_include_directories(${PRGNAME} PRIVATE ${FFTW3_INCLUDE_DIR_PARALLEL})
+      target_compile_definitions(${PRGNAME} PRIVATE "USE_FFTW_MPI")
-      target_compile_options(${PRGNAME} PRIVATE "-DUSE_FFTW_MPI")
     else()
       message(SEND_ERROR "MPI enabled but FFTW3 library not found with MPI support for long double precision!")
     endif()
   endif()
 
-  target_include_directories(${PRGNAME} PRIVATE ${MPI_CXX_INCLUDE_PATH})
+  target_link_libraries(${PRGNAME} PRIVATE MPI::MPI_CXX)
-  target_compile_options(${PRGNAME} PRIVATE "-DUSE_MPI")
+  target_compile_definitions(${PRGNAME} PRIVATE "USE_MPI")
-  target_link_libraries(${PRGNAME} ${MPI_LIBRARIES})
 endif(MPI_CXX_FOUND)
-
+# else()
-if(CODE_PRECISION STREQUAL "FLOAT" AND FFTW3_SINGLE_THREADS_FOUND) 
+if(CODE_PRECISION STREQUAL "FLOAT")
-  target_link_libraries(${PRGNAME} ${FFTW3_SINGLE_THREADS_LIBRARY})
+  if(FFTW3_SINGLE_THREADS_FOUND)
-  target_compile_options(${PRGNAME} PRIVATE "-DUSE_FFTW_THREADS")
+    target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_SINGLE_THREADS)
-elseif(CODE_PRECISION STREQUAL "DOUBLE" AND FFTW3_DOUBLE_THREADS_FOUND) 
+  endif()
-  target_link_libraries(${PRGNAME} ${FFTW3_DOUBLE_THREADS_LIBRARY})
+  target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_SINGLE_SERIAL) 
-  target_compile_options(${PRGNAME} PRIVATE "-DUSE_FFTW_THREADS")
+elseif(CODE_PRECISION STREQUAL "DOUBLE")
-elseif(CODE_PRECISION STREQUAL "LONGDOUBLE" AND FFTW3_LONGDOUBLE_THREADS_FOUND) 
+  if(FFTW3_DOUBLE_THREADS_FOUND) 
-  target_link_libraries(${PRGNAME} ${FFTW3_LONGDOUBLE_THREADS_LIBRARY})
+    target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_DOUBLE_THREADS)
-  target_compile_options(${PRGNAME} PRIVATE "-DUSE_FFTW_THREADS")
+  endif()
+  target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_DOUBLE_SERIAL) 
+elseif(CODE_PRECISION STREQUAL "LONGDOUBLE")
+  if(FFTW3_LONGDOUBLE_THREADS_FOUND) 
+    target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_LONGDOUBLE_THREADS)
+  endif()
+  target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_LONGDOUBLE_SERIAL)
 endif()
 
+
 if(HDF5_FOUND)
   # target_link_libraries(${PRGNAME} ${HDF5_C_LIBRARY_DIRS})
-  target_link_libraries(${PRGNAME} ${HDF5_LIBRARIES})
+  target_link_libraries(${PRGNAME} PRIVATE ${HDF5_LIBRARIES})
   target_include_directories(${PRGNAME} PRIVATE ${HDF5_INCLUDE_DIRS})
-  target_compile_options(${PRGNAME} PRIVATE "-DUSE_HDF5")
+  target_compile_definitions(${PRGNAME} PRIVATE "USE_HDF5")
 endif(HDF5_FOUND)
 
 if(ENABLE_PANPHASIA)
-target_compile_options(${PRGNAME} PRIVATE "-DUSE_PANPHASIA")
+  target_compile_definitions(${PRGNAME} PRIVATE "USE_PANPHASIA")
 endif(ENABLE_PANPHASIA)
 
 if(ENABLE_PLT)
-target_compile_options(${PRGNAME} PRIVATE "-DENABLE_PLT")
+  target_compile_definitions(${PRGNAME} PRIVATE "ENABLE_PLT")
 endif(ENABLE_PLT)
 
-target_link_libraries(${PRGNAME} ${FFTW3_LIBRARIES})
+target_link_libraries(${PRGNAME} PRIVATE GSL::gsl)
-target_include_directories(${PRGNAME} PRIVATE ${FFTW3_INCLUDE_DIRS})
 
-target_link_libraries(${PRGNAME} ${GSL_LIBRARIES})
+# GenericIO
-target_include_directories(${PRGNAME} PRIVATE ${GSL_INCLUDE_DIR})
+include(CMakeDependentOption)
+cmake_dependent_option(ENABLE_GENERICIO "Enable GenericIO (HACC) output support" off "ENABLE_MPI;MPI_CXX_FOUND" off)
 
-target_link_libraries(${PRGNAME} ${HDF5_LIBRARIES})
+if(ENABLE_GENERICIO)
-target_include_directories(${PRGNAME} PRIVATE ${HDF5_INCLUDE_DIR})
+  include(${CMAKE_CURRENT_SOURCE_DIR}/external/genericio.cmake)
+  target_link_libraries(${PRGNAME} PRIVATE genericio::genericio_mpi)
+  target_compile_definitions(${PRGNAME} PRIVATE "ENABLE_GENERICIO")
+endif()

FindFFTW3.cmake

@@ -245,3 +245,14 @@ find_package_handle_standard_args(FFTW3
     VERSION_VAR FFTW3_VERSION_STRING
     HANDLE_COMPONENTS
 )
+
+if(FFTW3_FOUND)
+  foreach(component ${FFTW3_FIND_COMPONENTS})
+    if(NOT TARGET FFTW3::FFTW3_${component})
+      add_library(FFTW3::FFTW3_${component} UNKNOWN IMPORTED)
+      set_target_properties(FFTW3::FFTW3_${component} PROPERTIES 
+        IMPORTED_LOCATION "${FFTW3_${component}_LIBRARY}"
+        INTERFACE_INCLUDE_DIRECTORIES ${FFTW3_INCLUDE_DIR})
+    endif()
+  endforeach()
+endif()

external/class.cmake

@@ -104,10 +104,10 @@ if(ENABLE_CLASS)
     target_include_directories(class PRIVATE ${CLASS_INCLUDE_DIR})
     target_include_directories(class PRIVATE ${CLASS_INCLUDE_CPP_DIR})
     target_compile_options(class PRIVATE "-ffast-math")
-    target_compile_options(class PRIVATE "-D__CLASSDIR__=\"${CMAKE_CURRENT_LIST_DIR}/class\"")
+    target_compile_definitions(class PRIVATE "__CLASSDIR__=\"${CMAKE_CURRENT_LIST_DIR}/class\"")
     set_property(TARGET class PROPERTY POSITION_INDEPENDENT_CODE ON)
     target_include_directories(class PRIVATE ${CMAKE_CURRENT_LIST_DIR}/class/hyrec)
-    target_compile_options(class PRIVATE "-DHYREC")
+    target_compile_definitions(class PRIVATE "HYREC")
     set_target_properties(class PROPERTIES CXX_STANDARD 14 C_STANDARD 11)
     # target_compile_options(class PRIVATE "-Wall")
     # target_compile_options(class PRIVATE "-Wextra")
@@ -127,9 +127,9 @@ macro(target_setup_class target_name)
       target_link_libraries(${target_name} class_cpp)
       add_dependencies(${target_name} class_objects)
     else(USE_CLASS_MAKEFILE)
-      target_link_libraries(${target_name} class)
+      target_link_libraries(${target_name} PRIVATE class)
     endif(USE_CLASS_MAKEFILE)
-    target_compile_options(${target_name} PRIVATE "-DUSE_CLASS")
+    target_compile_definitions(${target_name} PRIVATE "USE_CLASS")
   endif(ENABLE_CLASS)
 endmacro(target_setup_class)
 

external/genericio.cmake

@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 3.11)
+include(FetchContent)
+FetchContent_Declare(
+    genericio
+    GIT_REPOSITORY https://xgitlab.cels.anl.gov/mbuehlmann/genericio.git
+    GIT_TAG master
+)
+
+FetchContent_GetProperties(genericio)
+if(NOT genericio_POPULATED)
+    FetchContent_Populate(genericio)
+    add_subdirectory(${genericio_SOURCE_DIR} ${genericio_BINARY_DIR})
+endif()

include/general.hh

@@ -28,12 +28,12 @@
 #include <fftw3.h>
 #endif
 
-#include <config_file.hh>
+#include "config_file.hh"
 
 #define _unused(x) ((void)(x))
 
 // include CMake controlled configuration settings
-#include <cmake_config.hh>
+#include "cmake_config.hh"
 
 #if defined(USE_PRECISION_FLOAT)
 using real_t = float;

include/grid_fft.hh

@@ -80,12 +80,13 @@ public:
     // avoid implicit copying of data
     Grid_FFT(const grid_fft_t &g) = delete;
 
-    ~Grid_FFT()
+    ~Grid_FFT() { reset(); }
+
+    void reset()
     {
-        if (data_ != nullptr)
+        if (data_ != nullptr) { fftw_free(data_); }
-        {
+        if (plan_ != nullptr) { fftw_destroy_plan(plan_); }
-            fftw_free(data_);
+        if (iplan_ != nullptr) { fftw_destroy_plan(iplan_); }
-        }
     }
 
     const grid_fft_t *get_grid(size_t ilevel) const { return this; }

include/ic_generator.hh

@@ -29,6 +29,8 @@ namespace ic_generator{
     
     int Initialise( config_file& the_config );
 
+    void reset();
+
     extern std::unique_ptr<RNG_plugin> the_random_number_generator;
     extern std::unique_ptr<output_plugin> the_output_plugin;
     extern std::unique_ptr<cosmology::calculator>  the_cosmo_calc;

src/grid_fft.cc

@@ -386,10 +386,7 @@ void Grid_FFT<data_t, bdistributed>::Read_from_HDF5(const std::string Filename,
         this->n_[i] = dimsize[i];
     this->space_ = rspace_id;
 
-    if (data_ != nullptr)
+    this->reset();
-    {
-        fftw_free(data_);
-    }
     this->allocate();
 
     //... copy data to internal array ...

src/ic_generator.cc

@@ -49,6 +49,12 @@ int Initialise( config_file& the_config )
     return 0;
 }
 
+void reset () {
+    the_random_number_generator.reset();
+    the_output_plugin.reset();
+    the_cosmo_calc.reset();
+}
+
 int Run( config_file& the_config )
 {
     //--------------------------------------------------------------------------------------------------------

src/main.cc

@@ -244,6 +244,9 @@ int main( int argc, char** argv )
     // particle::test_plt();
     ///////////////////////////////////////////////////////////////////////
 
+    // call the destructor of plugins before tearing down MPI
+    ic_generator::reset();
+
 #if defined(USE_MPI)
     MPI_Barrier(MPI_COMM_WORLD);
     MPI_Finalize();

src/plugins/output_genericio.cc

@@ -0,0 +1,137 @@
+#ifdef ENABLE_GENERICIO
+
+#include <output_plugin.hh>
+#include <GenericIO.h>
+#include <mpi.h>
+
+class genericio_output_plugin : public output_plugin
+{
+protected:
+    real_t lunit_, vunit_;
+    bool hacc_hydro_;
+    float hacc_etamax_;
+    float hh_value_, rho_value_, mu_value_;
+    std::vector<int64_t> ids;
+    std::vector<float> xx, yy, zz;
+    std::vector<float> vx, vy, vz;
+    std::vector<float> mass, hh, uu;
+    std::vector<float> mu, phi, rho;
+    std::vector<int16_t> mask;
+
+public:
+    //! constructor
+    explicit genericio_output_plugin(config_file &cf)
+            : output_plugin(cf, "GenericIO")
+    {
+        real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
+        lunit_ = cf_.get_value<double>("setup", "BoxLength");
+        vunit_ = lunit_;
+        hacc_hydro_ = cf_.get_value_safe<bool>("output", "GenericIO_HACCHydro", false);
+        hacc_etamax_ = cf_.get_value_safe<float>("output", "GenericIO_ETAMAX", 1.0f);
+        hh_value_ = 4.0f * hacc_etamax_ * lunit_ / cf_.get_value<float>("setup", "GridRes");
+        mu_value_ = 4.0 / (8.0 - 5.0 * (1.0 - 0.75)); // neutral value. FIXME: account for ionization?
+        
+        double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
+        rho_value_ = cf_.get_value<double>("cosmology", "Omega_b") * rhoc;
+    }
+
+    output_type write_species_as(const cosmo_species &) const { return output_type::particles; }
+
+    real_t position_unit() const { return lunit_; }
+
+    real_t velocity_unit() const { return vunit_; }
+
+    bool has_64bit_reals() const { return false; }
+
+    bool has_64bit_ids() const { return true; }
+
+    void write_particle_data(const particle::container &pc, const cosmo_species &s, double Omega_species)
+    {
+        double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
+        double boxmass = Omega_species * rhoc * std::pow(cf_.get_value<double>("setup", "BoxLength"), 3.);
+        double particle_mass = boxmass / pc.get_global_num_particles();
+
+        size_t npart = pc.get_local_num_particles();
+        xx.reserve(xx.size() + npart);
+        yy.reserve(yy.size() + npart);
+        zz.reserve(zz.size() + npart);
+        vx.reserve(vx.size() + npart);
+        vy.reserve(vy.size() + npart);
+        vz.reserve(vz.size() + npart);
+        ids.reserve(ids.size() + npart);
+        
+        
+
+        auto _pos = reinterpret_cast<const float*>(pc.get_pos32_ptr());
+        auto _vel = reinterpret_cast<const float*>(pc.get_vel32_ptr());
+        auto _ids = reinterpret_cast<const uint64_t*>(pc.get_ids64_ptr());
+        
+        for(size_t i=0; i<npart; ++i) {
+            xx.push_back(_pos[3*i + 0]);
+            yy.push_back(_pos[3*i + 1]);
+            zz.push_back(_pos[3*i + 2]);
+            vx.push_back(_vel[3*i + 0]);
+            vy.push_back(_vel[3*i + 1]);
+            vz.push_back(_vel[3*i + 2]);
+        }
+
+        // phi doesn't need to be initialized, just needs to be present in data
+        phi.resize(phi.size() + npart);
+        std::copy(_ids, _ids+npart, std::back_inserter(ids));
+
+        if(hacc_hydro_) {
+            size_t prev_size = mass.size();
+            size_t new_size = prev_size + npart;
+            
+            mass.resize(new_size);
+            std::fill(mass.begin() + prev_size, mass.end(), particle_mass);
+
+            mask.resize(new_size);
+            std::fill(mask.begin() + prev_size, mask.end(), s == cosmo_species::baryon ? 1<<2 : 0);
+
+            hh.resize(new_size);
+            std::fill(hh.begin() + prev_size, hh.end(), s == cosmo_species::baryon ? hh_value_ : 0.0f);
+
+            uu.resize(new_size);
+            std::fill(uu.begin() + prev_size, uu.end(), s == cosmo_species::baryon ? 0.0f : 0.0f);
+
+            rho.resize(new_size);
+            std::fill(rho.begin() + prev_size, rho.end(), s == cosmo_species::baryon ? rho_value_ : 0.0f);
+
+            mu.resize(new_size);
+            std::fill(mu.begin() + prev_size, mu.end(), s == cosmo_species::baryon ? mu_value_ : 0.0f);
+        }
+    }
+
+    ~genericio_output_plugin() override {
+        gio::GenericIO writer(MPI_COMM_WORLD, fname_, gio::GenericIO::FileIO::FileIOMPI);
+        writer.setPhysOrigin(0., -1);
+        writer.setPhysScale(lunit_, -1);
+        writer.setNumElems(xx.size());
+        writer.addVariable("x", xx);
+        writer.addVariable("y", yy);
+        writer.addVariable("z", zz);
+        writer.addVariable("vx", vx);
+        writer.addVariable("vy", vy);
+        writer.addVariable("vz", vz);
+        writer.addVariable("id", ids);
+        writer.addVariable("phi", phi);
+        if(hacc_hydro_) {
+            writer.addVariable("mass", mass);
+            writer.addVariable("mask", mask);
+            writer.addVariable("hh", hh);
+            writer.addVariable("uu", uu);
+            writer.addVariable("rho", rho);
+            writer.addVariable("mu", mu);
+        }
+        writer.write();
+        MPI_Barrier(MPI_COMM_WORLD);
+    }
+};
+
+namespace
+{
+output_plugin_creator_concrete<genericio_output_plugin> creator("genericio");
+} // namespace
+
+#endif // ENABLE_GENERICIO