Tabulated optically thin radiative cooling module with Townsend integration

CMakeLists.txt

@@ -98,7 +98,7 @@ endif()
 if(Idefix_MPI)
   add_compile_definitions("WITH_MPI")
   find_package(MPI REQUIRED)
-  target_link_libraries(idefix MPI::MPI_CXX)
+  target_link_libraries(idefix PRIVATE MPI::MPI_CXX)
   target_sources(idefix
     PUBLIC src/mpi.cpp
     PUBLIC src/mpi.hpp
@@ -114,9 +114,9 @@ if(Idefix_HDF5)
     PUBLIC src/output/xdmf.cpp
     PUBLIC src/output/xdmf.hpp
   )
-  find_package(HDF5 REQUIRED)
-  target_link_libraries(idefix "${HDF5_LIBRARIES}")
-  target_include_directories(idefix "${HDF5_INCLUDE_DIRS}")
+  find_package(HDF5 REQUIRED COMPONENTS C)
+  target_link_libraries(idefix PRIVATE HDF5::HDF5)
+  # target_include_directories(idefix PRIVATE ${HDF5_INCLUDE_DIRS})
   message(STATUS "XDMF (hdf5+xmf) dumps enabled")
 else()
   set(Idefix_HDF5 OFF)
@@ -237,7 +237,7 @@ target_include_directories(idefix PUBLIC
                            src
                            )
 
-target_link_libraries(idefix Kokkos::kokkos)
+target_link_libraries(idefix PRIVATE Kokkos::kokkos)
 
 message(STATUS "Idefix final configuration")
 if(Idefix_EVOLVE_VECTOR_POTENTIAL)

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,62 @@
+.. _radiativeCoolingModule:
+
+Radiative Cooling module
+===================
+
+Equations solved and method
+---------------------------
+
+The ``RadiativeCooling`` module implements the computation of the loss of internal thermal energy
+due 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 seperately 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 the 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 ``definitions.hpp``. For example,
+    ``
+    #define     UNIT_LENGTH     3.086e+18
+    #define     UNIT_DENSITY    (1.0e-02*0.609*1.673e-24)
+    #define     UNIT_VELOCITY   1.000e+05
+    ``
+
+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.
+
+``
+Cooling   Tabulated    cooltable.dat    Townsend    TcoolFloor    1.0e+04
+``
+
++----------------------+-------------------------+----------------------------------------------------------------------------------------------+
+|  Entry name          | Parameter type          | Comment                                                                                      |
++======================+=========================+==============================================================================================+
+| cooling mode         | string                  | | Type of radiative cooling. Only `Tabulated` 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 (skip)    | string                  | | 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));
                 },

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
   //*****************************************************************
@@ -194,7 +207,10 @@ struct Fluid_AddSourceTermsFunctor {
       Uc(MX2,k,j,i) += dt*Sm / rt(i);
   #endif // COMPONENTS
 #endif
+    if (coolingOn) {
+      Uc(ENG,k,j,i) += (delta_eng_cool(k,j,i)/Uc(RHO,k,j,i)); // specific energy
     }
+  }
 };
 
 

src/fluid/cooling/CMakeLists.txt

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