Tabulated optically thin radiative cooling module with Townsend integration

CMakeLists.txt

@@ -253,7 +253,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
   //*****************************************************************
@@ -191,7 +204,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); // energy per unit volume
     }
+  }
 };
 
 

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}
+  )

src/fluid/cooling/cooling.cpp

@@ -0,0 +1,210 @@
+// ***********************************************************************************
+// Idefix MHD astrophysical code
+// Copyright(C) Geoffroy R. J. Lesur <geoffroy.lesur@univ-grenoble-alpes.fr>
+// and other code contributors
+// Licensed under CeCILL 2.1 License, see COPYING for more information
+// ***********************************************************************************
+
+// This source code is largely inspired from the viscous_flux of Pluto4.2
+// ((c) P. Tzeferacos & A. Mignone)
+
+// Implementation of monotonicity-preserving viscous flux following ZuHone et al.,
+// ApJ
+
+#include <string>
+#include <sstream>
+
+#include "idefix.hpp"
+#include "cooling.hpp"
+#include "dataBlock.hpp"
+#include "fluid.hpp"
+
+void RadCooling::InitUnits() {
+  #if defined(UNIT_DENSITY) || defined(UNIT_MASS)
+  #ifdef UNIT_DENSITY
+  rho_unit = UNIT_DENSITY;
+  #endif
+  #ifdef UNIT_MASS
+  mass_unit = UNIT_MASS;
+  #endif
+  #else
+  #define UNIT_DENSITY    1.0
+  rho_unit = 1.0;
+  IDEFIX_WARNING("UNIT_DENSITY or UNIT_MASS "
+                 "missing in definitions.hpp");
+  #endif
+
+  #if defined(UNIT_VELOCITY) || defined(UNIT_TIME)
+  #ifdef UNIT_VELOCITY
+  vel_unit = UNIT_VELOCITY;
+  #endif
+  #ifdef UNIT_TIME
+  time_unit = UNIT_TIME;
+  #endif
+  #else
+  #define UNIT_VELOCITY    1.0
+  vel_unit = 1.0;
+  IDEFIX_WARNING("UNIT_VELOCITY or UNIT_TIME "
+                 "missing in definitions.hpp");
+  #endif
+
+  #ifdef UNIT_LENGTH
+  len_unit = UNIT_LENGTH;
+  #else
+  #define UNIT_LENGTH    1.0
+  len_unit = 1.0;
+  IDEFIX_WARNING("UNIT_LENGTH "
+                 "missing in definitions.hpp");
+  #endif
+
+  #ifndef UNIT_DENSITY
+  rho_unit = mass_unit/pow(len_unit,3);
+  #endif
+  #ifndef UNIT_MASS
+  mass_unit = rho_unit*pow(len_unit,3);
+  #endif
+
+  #ifndef UNIT_VELOCITY
+  vel_unit = len_unit/time_unit;
+  #endif
+  #ifndef UNIT_TIME
+  time_unit = len_unit/vel_unit;
+  #endif
+}
+
+void RadCooling::InitArrays() {
+  // Allocate and fill arrays when neede
+  cooltable_data = LookupTable<1>(cooltable_filename,' ');
+  temperature_data = cooltable_data.xinDev;
+  Lambda_cool_data = cooltable_data.dataDev;
+
+  delta_eng = IdefixArray3D<real>("delta_eng_cooling_source", data->np_tot[KDIR],
+                                                              data->np_tot[JDIR],
+                                                              data->np_tot[IDIR]);
+}
+void RadCooling::ShowConfig() {
+  if(cooling_type.compare("tabulated")==0) {
+    idfx::cout << "Optically thin radiative cooling: ENABLED"
+               << std::endl;
+    if(cooling_integration.compare("townsend") == 0) {
+      idfx::cout << "Radiative cooling Integration: Townsend"
+                 << std::endl;
+    } else {
+      IDEFIX_ERROR("Unknown radiative cooling integration mode");
+    }
+  } else {
+    IDEFIX_ERROR("Unknown radiative cooling mode");
+  }
+  /*
+    IDEFIX_WARNING("A sample warning "
+                   "to be used with Idefix.");
+  */
+}
+
+// This function computes the pressure drop due to optically thin radiative cooling
+// using Townsend integration
+void RadCooling::TownsendIntegration(real dt) {
+  idfx::pushRegion("RadCooling::TownsendIntegration");
+  IdefixArray1D<real> x1 = this->data->x[IDIR];
+  IdefixArray1D<real> x2 = this->data->x[JDIR];
+  IdefixArray1D<real> x3 = this->data->x[KDIR];
+  IdefixArray3D<real> dV = this->data->dV;
+
+  real kB = 1.380649e-16;
+  real mu = 0.609;
+  real mp = 1.6726e-24;
+  real XH = 0.71;
+
+  int ibeg, iend, jbeg, jend, kbeg, kend;
+  ibeg = this->data->beg[IDIR];
+  iend = this->data->end[IDIR];
+  jbeg = this->data->beg[JDIR];
+  jend = this->data->end[JDIR];
+  kbeg = this->data->beg[KDIR];
+  kend = this->data->end[KDIR];
+
+  idefix_for("RadCoolingLoop",kbeg, kend, jbeg, jend, ibeg, iend,
+    KOKKOS_LAMBDA (int k, int j, int i) {
+      int T_indx_lo = 0, T_indx_hi=temperature_data.extent(0)-1;
+      int T_indx_mid;
+      // ideal gas eos is used
+      real temperature = Vc(PRS,k,j,i)/Vc(RHO,k,j,i)*(mu*mp/kB)*pow(vel_unit,2);
+
+      if ( (temperature<temperature_data(0)) ||
+           (temperature>temperature_data(temperature_data.extent(0)-1)) ) {
+        std::ostringstream tmp;
+        tmp << std::scientific <<
+            "Temperature out of range: T="
+            << temperature << " K" << std::endl;
+        IDEFIX_ERROR(tmp.str());
+      }
+
+      while (T_indx_lo<=T_indx_hi) {
+        T_indx_mid = (T_indx_lo + T_indx_hi)/2;
+        if (temperature < temperature_data(T_indx_mid)) {
+          T_indx_hi = T_indx_mid-1;
+        } else if (temperature > temperature_data(T_indx_mid)) {
+          T_indx_lo = T_indx_mid+1;
+        } else {
+          T_indx_lo = T_indx_mid;
+          T_indx_hi = T_indx_mid;
+          break;
+        }
+      }
+      if (T_indx_lo!=T_indx_hi) {
+        // swap
+        T_indx_mid = T_indx_lo;
+        T_indx_lo = T_indx_hi;
+        T_indx_hi = T_indx_mid;
+      }
+
+      real temperature_lo = temperature_data(T_indx_lo);
+      real temperature_hi = temperature_data(T_indx_hi);
+      real Lambda_lo = Lambda_cool_data(T_indx_lo);
+      real Lambda_hi = Lambda_cool_data(T_indx_hi);
+      // T_ref = temperature_hi
+      real alpha = std::log(Lambda_hi/Lambda_lo)/std::log(temperature_hi/temperature_lo);
+      real Lambda_T = Lambda_lo * pow((temperature/temperature_lo), alpha);
+      // real gamma = eos.GetGamma(Vc(PRS,k,j,i),Vc(RHO,k,j,i));
+      real eint = eos.GetInternalEnergy(Vc(PRS,k,j,i),
+                                        Vc(RHO,k,j,i))
+                                       *rho_unit*pow(vel_unit,2); // cgs
+      real edot = pow(Vc(RHO,k,j,i)*rho_unit*XH/mp,2)*Lambda_T; // cgs
+      real t_cool = eint/edot; // cgs
+
+      real Y, del_prs;
+      if (alpha!=1.0) {
+        Y = (1.0-pow(temperature_hi/temperature, alpha-1.0))/(1.0-alpha);
+      } else {
+        Y = std::log(temperature_hi/temperature);
+      }
+      real inv_Y_arg = Y + (temperature/temperature_hi)
+                         * (Lambda_hi/Lambda_T)
+                         * (dt/(t_cool/(len_unit/vel_unit)));
+      real T_fin;
+      if (alpha!=1.0) {
+        T_fin = temperature_hi*pow((1.0 - (1.0-alpha)*inv_Y_arg), 1.0/(1.0-alpha));
+      } else {
+        T_fin = temperature_hi * std::exp(-inv_Y_arg);
+      }
+      // ideal gas eos is used
+      if (T_fin>=TcoolFloor) {
+        del_prs = -Vc(RHO,k,j,i)/(mu*mp/kB)*(temperature-T_fin)/pow(vel_unit,2);
+        delta_eng(k,j,i) = eos.GetInternalEnergy(del_prs, Vc(RHO,k,j,i));
+      }
+      else {
+        del_prs = ZERO_F;
+        delta_eng(k,j,i) = ZERO_F;
+      }
+    });
+  idfx::popRegion();
+}
+
+// This function computes the energy loss in the conservative variable due
+// to optically thin radiative cooling
+void RadCooling::CalculateCoolingSource(real dt) {
+  idfx::pushRegion("RadCooling::CalculateCoolingSource");
+  // calculate the change in internal thermal energy
+  rad_cooling_int_func(dt);
+  idfx::popRegion();
+}