Tabulated optically thin radiative cooling module with Townsend integration

.github/workflows/idefix-ci-jobs.yml

@@ -55,6 +55,8 @@ jobs:
         run: scripts/ci/run-tests $IDEFIX_DIR/test/HD/ViscousDisk -all $TESTME_OPTIONS
       - name: Thermal diffusion
         run: scripts/ci/run-tests $IDEFIX_DIR/test/HD/thermalDiffusion -all $TESTME_OPTIONS
+      - name: CoolBox
+        run: scripts/ci/run-tests $IDEFIX_DIR/test/HD/CoolBox -all $TESTME_OPTIONS
 
   ShocksMHD:
     runs-on: self-hosted

CMakeLists.txt

@@ -48,6 +48,31 @@ if(Kokkos_ENABLE_CUDA)
   set(Kokkos_ENABLE_CUDA_LAMBDA ON CACHE BOOL "Idefix requires lambdas on Cuda" FORCE)
   # CUDA_MALLOC_ASYNC disbaled by default in Kokkos 4.5, so not required here
   #set(Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC OFF CACHE BOOL "Disable Async malloc to avoid bugs on PSM2" FORCE)
+
+  # On recent glibc (>= 2.34) the realtime functions were merged into libc and
+  # /usr/lib64/librt.a is only an empty archive stub. CMake's FindCUDAToolkit still
+  # resolves librt via find_library(CUDAToolkit_rt_LIBRARY rt) and attaches the
+  # resulting /usr/lib64/librt.a to CUDA::cudart_static_deps, which is pulled in by
+  # CUDA::cudart and CUDA::cuda_driver (the targets Kokkos links). The device link
+  # then passes that absolute archive path to nvlink, which cannot open the empty
+  # stub and fails with:
+  #   nvlink fatal: Could not open input file '/usr/lib64/librt.a'
+  # Pre-seed the cache so FindCUDAToolkit links the shared librt (host-only, ignored
+  # by nvlink) instead of the static stub. We point at the shared object when it can
+  # be located, otherwise fall back to the "-lrt" link flag form; both avoid handing
+  # an absolute *.a path to nvlink, and the realtime symbols are provided by libc.
+  if(NOT DEFINED CUDAToolkit_rt_LIBRARY)
+    find_library(IDEFIX_RT_SHARED_LIBRARY
+                 NAMES librt.so.1 librt.so
+                 PATHS /lib64 /usr/lib64 /lib /usr/lib)
+    if(IDEFIX_RT_SHARED_LIBRARY)
+      set(CUDAToolkit_rt_LIBRARY "${IDEFIX_RT_SHARED_LIBRARY}"
+          CACHE FILEPATH "Shared librt used by CUDA::cudart_static_deps" FORCE)
+    else()
+      set(CUDAToolkit_rt_LIBRARY "rt"
+          CACHE FILEPATH "librt link flag used by CUDA::cudart_static_deps" FORCE)
+    endif()
+  endif()
 endif()
 
 # Add kokkos CMAKE files (required early since these set compiler options)
@@ -113,12 +138,64 @@ if(Idefix_HDF5)
     PUBLIC src/output/xdmf.cpp
     PUBLIC src/output/xdmf.hpp
   )
-  find_package(HDF5 REQUIRED)
-  target_link_libraries(idefix "${HDF5_LIBRARIES}")
+  # Prefer imported targets (matches older working behavior), then fall back.
+  find_package(HDF5 QUIET COMPONENTS C)
+  if(NOT HDF5_FOUND)
+    find_package(HDF5 REQUIRED MODULE COMPONENTS C)
+  endif()
+
+  set(_idefix_hdf5_link_items "${HDF5_LIBRARIES}")
+  if(TARGET hdf5::hdf5)
+    set(_idefix_hdf5_link_items hdf5::hdf5)
+  elseif(TARGET HDF5::HDF5)
+    set(_idefix_hdf5_link_items HDF5::HDF5)
+  endif()
+
+  target_link_libraries(idefix ${_idefix_hdf5_link_items})
   message(STATUS "Found HDF5 include directories: ${HDF5_INCLUDE_DIRS}")
   target_include_directories(idefix PUBLIC "${HDF5_INCLUDE_DIRS}")
   if(Idefix_MPI)
-    if(NOT HDF5_IS_PARALLEL)
+    include(CheckCSourceCompiles)
+
+    # Some CMake/HDF5 combinations do not reliably set HDF5_IS_PARALLEL.
+    # Combine metadata checks with a compile+link probe for MPI-IO symbols.
+    set(_idefix_hdf5_is_parallel FALSE)
+    if(HDF5_IS_PARALLEL OR HDF5_C_IS_PARALLEL)
+      set(_idefix_hdf5_is_parallel TRUE)
+    endif()
+
+    if(NOT _idefix_hdf5_is_parallel AND DEFINED HDF5_C_COMPILER_EXECUTABLE)
+      execute_process(
+        COMMAND ${HDF5_C_COMPILER_EXECUTABLE} -showconfig
+        OUTPUT_VARIABLE _idefix_hdf5_showconfig
+        ERROR_QUIET
+      )
+      if(_idefix_hdf5_showconfig MATCHES "Parallel HDF5:[ \t]*yes")
+        set(_idefix_hdf5_is_parallel TRUE)
+      endif()
+    endif()
+
+    set(CMAKE_REQUIRED_INCLUDES ${HDF5_INCLUDE_DIRS} ${MPI_C_INCLUDE_DIRS})
+    set(CMAKE_REQUIRED_LIBRARIES ${_idefix_hdf5_link_items} MPI::MPI_C)
+    check_c_source_compiles(
+      "#include <mpi.h>
+       #include <hdf5.h>
+       int main(void) {
+         hid_t plist = H5Pcreate(H5P_FILE_ACCESS);
+         H5Pset_fapl_mpio(plist, MPI_COMM_WORLD, MPI_INFO_NULL);
+         H5Pclose(plist);
+         return 0;
+       }"
+      IDEFIX_HDF5_HAS_MPI_IO
+    )
+    unset(CMAKE_REQUIRED_INCLUDES)
+    unset(CMAKE_REQUIRED_LIBRARIES)
+
+    if(IDEFIX_HDF5_HAS_MPI_IO)
+      set(_idefix_hdf5_is_parallel TRUE)
+    endif()
+
+    if(NOT _idefix_hdf5_is_parallel)
       message(FATAL_ERROR "Parallel HDF5 required for Idefix_MPI but the found HDF5 library does not support it")
     endif()
   endif()

doc/source/modules.rst

@@ -25,6 +25,10 @@ In this section, you will find a more detailed documentation about each module t
   The Braginskii module, models the anisotropic flux of heat and momentum
   taking place in weakly collisional, magnetised plasma (like the intracluster medium).
 
+:ref:`radiativeCoolingModule`
+  The optically thin radiative cooling module, handles the loss of internal thermal energy due to radiation in an
+  optically thin medium.
+
 :ref:`gridCoarseningModule`
   The grid coarsening module, that allows to derefine the grid in particular locations to speed up the computation.
 
@@ -42,5 +46,6 @@ In this section, you will find a more detailed documentation about each module t
    modules/eos.rst
    modules/selfGravity.rst
    modules/braginskii.rst
+   modules/radCooling.rst
    modules/gridCoarsening.rst
    modules/pydefix.rst

doc/source/modules/radCooling.rst

@@ -0,0 +1,57 @@
+.. _radiativeCoolingModule:
+
+Radiative Cooling module
+========================
+Equations solved and method
+---------------------------
+
+The ``RadiativeCooling`` module implements the computation of the loss of internal thermal energy
+due to radiation in an optically thin medium. Physically, it solves for :math:`\dot{e}=\mathcal{L}`,
+where we have used :math:`\mathcal{L}=-n_H^2 \Lambda(T)` (where :math:`T` is the gas temperature,
+:math:`n_H=\rho X_H/m_p` is the total hydrogen number density, and :math:`\Lambda(T)` is the
+radiative cooling rate computed separately from quantum mechanical calculations
+by other plasma modeling codes, for example, Cloudy (Ferland et. al, PASP 110, 749 (1998)).
+
+This computation becomes especially relevant for multiphase gas in astrophysical environments
+prevalent in the ISM, the CGM, and the ICM, for which this module has been designed.
+
+The ``RadiativeCooling`` module implemented in *Idefix* follows the algorithm of the Townsend
+to integrate the loss of internal thermal energy (Townsend, ApJS 181, 2 (2009)) at every timestep
+in an operator split manner. The cooling rate is read from a table at runtime where the `first` row
+is temperature (in :math:`\rm K`) and second  row is :math:`\Lambda (T)` (in :math:`\rm erg cm^3 s^{-1}`).
+
+.. note::
+    We assume a normalization of :math:`n_H`, the total hydrogen number density for the
+    of the cooling curve supplied by the rate table at runtime to *Idefix*. Different might cooling curves with
+    different normalisation is known to exist in literature and special attention must be given to
+    what is supplied to the code. Right now, this module has been tested only with the ideal gas equation of state.
+    We also assume a fixed composition with :math:`X_H=0.71` and mean particle mass :math:`\mu=0.609`, i.e., constant in the current implementation (appropriate
+    for fully ionized plasma).
+    It is recommended to include conversion factors between code and physical units in ``idefix.ini`` which is in the `[Units]`
+    block. This is to ensure that the cooling curve supplied in physical units is correctly interpreted by the code.
+    A test example of a uniform advecting box is in `test/HD/CoolBox`.
+Main parameters of the module
+-----------------------------
+
+The ``RadiativeCooling`` module is a submodule of the ``Hydro`` module to compute the loss of internal thermal energy (pressure)
+of the gas. The parameters specific to radiative cooling are to be set in a dedicated line starting with the word
+``Cooling`` in the ``[Hydro]`` block. An example is as follows succeded by a detailed explanation.
+
+.. code-block:: ini
+
+   Cooling   Tabulated    cooltable.dat    Townsend    TcoolFloor    1.0e+04
+
++----------------------+-------------------------+----------------------------------------------------------------------------------------------+
+|  Entry name          | Parameter type          | Comment                                                                                      |
++======================+=========================+==============================================================================================+
+| cooling mode         | string                  | | Type of radiative cooling. Only `Tabulated` is supported right now.                        |
++----------------------+-------------------------+----------------------------------------------------------------------------------------------+
+| table name           | string                  | | name/location of the cooling table w.r.t. *Idefix* binary to be loaded at runtime.         |
++----------------------+-------------------------+----------------------------------------------------------------------------------------------+
+| integration method   | string                  | | Integration method to calculate the internal thermal energy loss due to radiative cooling. |
+|                      |                         | | Only `Townsend` supported right now.                                                       |
++----------------------+-------------------------+----------------------------------------------------------------------------------------------+
+| TcoolFloor           | string (optional)       | | Floor temperature in K below which cooling is turned off.                                  |
++----------------------+-------------------------+----------------------------------------------------------------------------------------------+
+| temperature floor    | float (optional)        | | Default is 1.0e+04                                                                         |
++----------------------+-------------------------+----------------------------------------------------------------------------------------------+

doc/source/programmingguide.rst

@@ -98,7 +98,7 @@ A typical loop on three indices looks like
   idefix_for("LoopName",
              kbeg,kend,
              jbeg,jend,
-             ibeg,ieng,
+             ibeg,iend,
              KOKKOS_LAMBDA (int k, int j, int i) {
                 myArray(k,j,i) = 0.0;
               });
@@ -153,7 +153,7 @@ For instance, a sum over all of the elements would be done through:
   idefix_reduce("Sum",
                 kbeg,kend,
                 jbeg,jend,
-                ibeg,ieng,
+                ibeg,iend,
                 KOKKOS_LAMBDA (int k, int j, int i, real &localSum) {
                     localSum += myArray(k,j,i);
                 },
@@ -179,7 +179,7 @@ snippet:
   idefix_reduce("Minimum",
                 kbeg,kend,
                 jbeg,jend,
-                ibeg,ieng,
+                ibeg,iend,
                 KOKKOS_LAMBDA (int k, int j, int i, real &localMin) {
                     localMin = std::fmin(localMin, myArray(k,j,i));
                 },

pytools/idfx_test.py

@@ -130,13 +130,18 @@ def __init__ (self, current_test_file, name=""):
 
     # transform all arguments from args into attributes of this instance
     self.__dict__.update(vars(args))
+    idefix_dir = os.path.realpath(os.path.abspath(self.idefixDir))
+    current_test_file_abs = os.path.realpath(os.path.abspath(current_test_file))
+
     # store the full path of problem directory
-    self.currentTestFile = current_test_file
+    self.currentTestFile = current_test_file_abs
     self.currentTestName = name
-    self.problemDir=os.path.dirname(current_test_file)
-    self.referenceDirectory = os.path.join(idefix_dir_env,"reference")
+    self.problemDir=os.path.dirname(current_test_file_abs)
+    self.referenceDirectory = os.path.join(idefix_dir,"reference")
     # current directory relative to $IDEFIX_DIR/test (used to retrieve the path ot reference files)
-    self.testDir=os.path.relpath(self.problemDir,os.path.join(idefix_dir_env,"test"))
+    self.testDir=os.path.relpath(self.problemDir,os.path.join(idefix_dir,"test"))
+    # Default input file name used by Idefix when no -i argument is provided.
+    self.inifile = "idefix.ini"
     # build directory, currently inside the test named build-test
     self.buildDir=os.path.join(self.problemDir,"build-test")
     # remind what build we dit last
@@ -343,8 +348,11 @@ def run(self, inputFile="", np=2, nowrite=False, restart=-1):
 
       comm=[os.path.join(self.buildDir,"idefix")]
       if inputFile:
+          self.inifile = os.path.basename(inputFile)
           comm.append("-i")
           comm.append(inputFile)
+      else:
+          self.inifile = "idefix.ini"
       if self.mpi:
           if self.dec:
             np=1

src/fluid/CMakeLists.txt

@@ -4,6 +4,7 @@ add_subdirectory(constrainedTransport)
 add_subdirectory(eos)
 add_subdirectory(RiemannSolver)
 add_subdirectory(tracer)
+add_subdirectory(cooling)
 
 
 target_sources(idefix

src/fluid/addSourceTerms.hpp

@@ -19,7 +19,8 @@ struct Fluid_AddSourceTermsFunctor {
   //*****************************************************************
   // Functor constructor
   //*****************************************************************
-  explicit Fluid_AddSourceTermsFunctor(Fluid<Phys> *hydro, real dt) {
+  explicit Fluid_AddSourceTermsFunctor(Fluid<Phys> *hydro, real dt):
+    hydroin(hydro) {
     Uc = hydro->Uc;
     Vc = hydro->Vc;
     x1 = hydro->data->x[IDIR];
@@ -42,6 +43,13 @@ struct Fluid_AddSourceTermsFunctor {
     }
     // shearing box (only with fargo&cartesian)
     sbS = hydro->sbS;
+
+    // Radiative cooling
+    coolingOn = hydro->coolingOn;
+    if (coolingOn) {
+      hydro->radCooling->CalculateCoolingSource(dt);
+      this->delta_eng_cool = hydro->radCooling->delta_eng;
+    }
   }
 
   //*****************************************************************
@@ -52,7 +60,9 @@ struct Fluid_AddSourceTermsFunctor {
   IdefixArray1D<real> x1;
   IdefixArray1D<real> x2;
   IdefixArray3D<real> csIsoArr;
+  IdefixArray3D<real> delta_eng_cool;
 
+  Fluid<Phys> *hydroin;
   real dt;
 #if GEOMETRY == SPHERICAL
   IdefixArray1D<real> sinx2;
@@ -72,6 +82,9 @@ struct Fluid_AddSourceTermsFunctor {
   // shearing box (only with fargo&cartesian)
   real sbS;
 
+  // Radiative cooling
+  bool coolingOn;
+
   //*****************************************************************
   // Functor Operator
   //*****************************************************************
@@ -191,7 +204,12 @@ struct Fluid_AddSourceTermsFunctor {
       Uc(MX2,k,j,i) += dt*Sm / rt(i);
   #endif // COMPONENTS
 #endif
+#if HAVE_ENERGY
+    if (coolingOn) {
+      Uc(ENG,k,j,i) += delta_eng_cool(k,j,i); // energy per unit volume
     }
+#endif
+  }
 };
 
 

src/fluid/convertConsToPrim.hpp

@@ -108,7 +108,7 @@ void Fluid<Phys>::ConvertConsToPrim() {
 
   IdefixArray4D<real> Vc = this->Vc;
   IdefixArray4D<real> Uc = this->Uc;
-  EquationOfState eos;
+  EquationOfState eos{};
   if constexpr(Phys::eos) {
     eos = *(this->eos.get());
   }
@@ -155,7 +155,7 @@ void Fluid<Phys>::ConvertPrimToCons() {
 
   IdefixArray4D<real> Vc = this->Vc;
   IdefixArray4D<real> Uc = this->Uc;
-  EquationOfState eos;
+  EquationOfState eos{};
   if constexpr(Phys::eos) {
     eos = *(this->eos.get());
   }

src/fluid/cooling/CMakeLists.txt

@@ -0,0 +1,9 @@
+target_sources(idefix
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/cooling.hpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/coolingConstructor.hpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/cooling.cpp
+  )
+
+target_include_directories(idefix
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}
+  )