Add units

doc/source/reference/idefix.ini.rst

@@ -265,8 +265,10 @@ section as followed:
 +----------------+-------------------------+---------------------------------------------------------------------------------------------+
 | Mcentral       | real                    | | Mass of the central object when a central potential is enabled (see above). Default is 1. |
 +----------------+-------------------------+---------------------------------------------------------------------------------------------+
-| gravCst        | real                    | | Set the value of the gravitational constant :math:`G_c` used by the central               |
+| gravCst        | real or string          | | Set the value of the gravitational constant :math:`G_c` used by the central               |
 |                |                         | | mass potential and self-gravitational potential (when enabled) ). Default is 1.           |
+|                |                         | | If set to "units", Idefix will compute the gravitational constant from the user units     |
+|                |                         | | defined in the units section (see :ref:`unitsSection`).                                   |
 +----------------+-------------------------+---------------------------------------------------------------------------------------------+
 | bodyForce      | string                  | | Adds an acceleration vector to each cell of the domain. The only value allowed            |
 |                |                         | | is ``userdef``. The ``Gravity`` class then expects a user-defined bodyforce function to   |
@@ -310,7 +312,28 @@ of gravitational potential in the ``Gravity`` section (see :ref:`gravitySection`
 |                |                         | | Default is 1 (i.e. self-gravity is computed at every cycle).                              |
 +----------------+-------------------------+---------------------------------------------------------------------------------------------+
 
+.. _unitsSection:
 
+``Units`` section
+------------------
+
+This optional section controls how the code handles units. By default, Idefix work "unit free" and this section is not useful.
+However, some physical processes like self-gravity or radiative transfer require the explicit use of units to define natural constants.
+These entries define how one code units should be converted to CGS units. When used, the following values should be set:
+
++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+
+|  Entry name    | Parameter type     | Comment                                                                                                   |
++================+====================+===========================================================================================================+
+| length         | float              | The code unit length: 1 code unit = `length` cm. 1 by default.                                            |
++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+
+| velocity       | float              | The code unit velocity: 1 code unit = `velocity` cm/s. 1 by default.                                      |
++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+
+| density        | float              | The code unit density: 1 code unit = `density` g/cm^3. 1 by default.                                      |
++----------------+--------------------+-----------------------------------------------------------------------------------------------------------+
+
+.. note::
+    The units module automatically reconstruct all of the units (time, temperature, magnetic field, etc.) from the above 3 fondamental constants. These can
+    all be obtained from the global object ``idfx::units``. Note however that the code outputs remain in code units, even if `[Units]` are defined.
 
 ``RKL`` section
 ------------------

src/CMakeLists.txt

@@ -28,6 +28,8 @@ target_sources(idefix
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/reduce.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/setup.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/setup.cpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/units.hpp
+  PUBLIC ${CMAKE_CURRENT_LIST_DIR}/units.cpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/timeIntegrator.hpp
   PUBLIC ${CMAKE_CURRENT_LIST_DIR}/timeIntegrator.cpp
   )

src/global.cpp

@@ -10,6 +10,7 @@
 #include "idefix.hpp"
 #include "global.hpp"
 #include "profiler.hpp"
+#include "units.hpp"
 
 #ifdef WITH_MPI
 #include "mpi.hpp"
@@ -28,6 +29,7 @@ IdefixOutStream cout;
 IdefixErrStream cerr;
 Profiler prof;
 LoopPattern defaultLoopPattern;
+Units units;
 
 #ifdef DEBUG
 static int regionIndent = 0;

src/global.hpp

@@ -18,6 +18,7 @@ void safeExit(int );       // Exit the code
 class IdefixOutStream;
 class IdefixErrStream;
 class Profiler;
+class Units;
 
 extern int prank;                       //< parallel rank
 extern int psize;
@@ -27,6 +28,7 @@ extern Profiler prof;                   //< profiler (for memory & performance u
 extern double mpiCallsTimer;            //< time significant MPI calls
 extern LoopPattern defaultLoopPattern;  //< default loop patterns (for idefix_for loops)
 extern bool warningsAreErrors;    //< whether warnings should be considered as errors
+extern Units units;               //< Units for the run
 
 void pushRegion(const std::string&);
 void popRegion();

src/gravity/gravity.cpp

@@ -11,13 +11,27 @@
 #include "planetarySystem.hpp"
 #include "dataBlock.hpp"
 #include "input.hpp"
+#include "units.hpp"
 
 Gravity::Gravity(Input &input, DataBlock *datain) {
   idfx::pushRegion("Gravity::Gravity");
   this->data = datain;
 
   // Gravitational constant G
-  this->gravCst = input.GetOrSet<real>("Gravity","gravCst",0, 1.0);
+  // should we compute it from the units?
+  if(input.CheckEntry("Gravity","gravCst")>=0) {
+    if(input.Get<std::string>("Gravity","gravCst",0).compare("units") == 0) {
+      // User ask us to compute the gravitational constants from the units
+
+      this->gravCst = idfx::units.GetDensity() * idfx::units.GetTime() * idfx::units.GetTime()
+                      * idfx::units.G;
+    } else {
+      this->gravCst = input.Get<real>("Gravity","gravCst",0);
+    }
+  } else { // default value to 1.0
+    this->gravCst = 1.0;
+  }
+
   // Gravitational potential
   int nPotential = input.CheckEntry("Gravity","potential");
   if(nPotential >=0) {

src/main.cpp

@@ -29,6 +29,7 @@
 #include "idefix.hpp"
 #include "profiler.hpp"
 #include "input.hpp"
+#include "units.hpp"
 #include "grid.hpp"
 #include "gridHost.hpp"
 #include "fluid.hpp"
@@ -78,6 +79,9 @@ int main( int argc, char* argv[] ) {
     input.PrintLogo();
     idfx::cout << "Main: initialization stage." << std::endl;
 
+    // Init the units when needed
+    idfx::units.Init(input);
+
     // Allocate the grid on device
     Grid grid(input);
     // Allocate the grid image on host
@@ -107,6 +111,7 @@ int main( int argc, char* argv[] ) {
                  << std::endl;
     }
     input.ShowConfig();
+    idfx::units.ShowConfig();
     grid.ShowConfig();
     data.ShowConfig();
     Tint.ShowConfig();

src/units.cpp

@@ -0,0 +1,55 @@
+// ***********************************************************************************
+// 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
+// ***********************************************************************************
+
+#include "units.hpp"
+#include "idefix.hpp"
+#include "input.hpp"
+
+
+void idfx::Units::Init(Input &input) {
+  if(input.CheckBlock("Units")) {
+    this->_length = input.GetOrSet<real>("Units","length",0,1.0);
+    this->_velocity = input.GetOrSet<real>("Units","velocity",0,1.0);
+    this->_density = input.GetOrSet<real>("Units","density",0,1.0);
+    this->ComputeUnits();
+  }
+}
+
+void idfx::Units::SetDensity(const real in) {
+  this->_density = in;
+  this->ComputeUnits();
+}
+
+void idfx::Units::SetLength(const real in) {
+  this->_length = in;
+  this->ComputeUnits();
+}
+
+void idfx::Units::SetVelocity(const real in) {
+  this->_velocity = in;
+  this->ComputeUnits();
+}
+
+void idfx::Units::ShowConfig() {
+  if(_isInitialized) {
+    idfx::cout << "Units: [Length]      = " << this->_length << " cm" << std::endl;
+    idfx::cout << "Units: [Velocity]    = " << this->_velocity << " cm/s" << std::endl;
+    idfx::cout << "Units: [Density]     = " << this->_density << " g/cm3" << std::endl;
+    idfx::cout << "Units: [Energy]     = " << this->_energy << " erg/cm3" << std::endl;
+    idfx::cout << "Units: [Time]     = " << this->_time << " s" << std::endl;
+    idfx::cout << "Units: [Temperature] = " << this->_Kelvin << " K" << std::endl;
+    idfx::cout << "Units: [Mag. Field]  = " << this->_magField << " G" << std::endl;
+  }
+}
+
+void idfx::Units::ComputeUnits() {
+  this->_isInitialized = true;
+  this->_magField = std::sqrt(4*M_PI*_density*_velocity*_velocity);
+  this->_Kelvin = u * _velocity*_velocity/k_B;
+  this->_energy = _density*_velocity*_velocity;
+  this->_time = _length/_velocity;
+}

src/units.hpp

@@ -0,0 +1,87 @@
+// ***********************************************************************************
+// 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
+// ***********************************************************************************
+
+#ifndef UNITS_HPP_
+#define UNITS_HPP_
+
+#include "idefix.hpp"
+
+// Units class implemented in CGS
+// Constants taken from Astropy
+class Input;
+namespace idfx {
+class Units {
+ public:
+  void Init(Input &input);
+
+
+  const real u{1.6605390666e-24};                   // Atomic mass unit (g)
+  const real m_p{1.67262192369e-24};                // Proton mass unit (g)
+  const real m_n{1.67492749804e-24};                // neutron mass unit (g)
+  const real m_e{9.1093837015e-28};                 // electron mass unit (g)
+  const real k_B{1.380649e-16};                     // Boltzmann constant (erg/K)
+  const real sigma_sb{5.6703744191844314e-05};      // Stephan Boltzmann constant (g/(K^4 s^3))
+  const real ar{7.5646e-15};                        // Radiation constant
+                                                    //      = 4*sigma_sb/c (g/(K^4 s^2 cm))
+  const real c{29979245800.0};                      // Speed of light (cm/s)
+  const real M_sun{1.988409870698051e+33};          // Solar mass (g)
+  const real R_sun{69570000000.0};                  // Solar radius (cm)
+  const real M_earth{5.972167867791379e+27};        // Earth mass (g)
+  const real R_earth{637810000.0};                  // Earth radius (cm)
+  const real G{6.674299999999999e-8};               // Gravitatonal constant  (cm3 / (g s2))
+  const real h{6.62607015e-27};                     // Planck constant (erg.s)
+  const real pc{3.08568e+18};                       // Parsec (cm)
+  const real au{1.49597892e13};                     // Astronomical unit (cm)
+  const real e{4.80320425e-10};                     // Elementary charge (cm^3/2 g^1/2 s^-1)
+
+
+  // User-defined units, non user-modifiable
+  // L (cm)  = L (code) * Units::length
+  KOKKOS_INLINE_FUNCTION real GetLength() const {return _length;}
+  // V(cm/s) = V(code) * Units::velocity
+  KOKKOS_INLINE_FUNCTION real GetVelocity() const {return _velocity;}
+  // density (g/cm^3) = density(code) * Units::density
+  KOKKOS_INLINE_FUNCTION real GetDensity() const {return _density;}
+
+  // Deduced units from user-defined units
+  // T(K) = P(code)/rho(code) * mu * Units::Kelvin
+  KOKKOS_INLINE_FUNCTION  real GetKelvin() const {return _Kelvin;}
+  // B(G) = B(code) * Units::MagField
+  KOKKOS_INLINE_FUNCTION  real GetMagField() const {return _magField;}
+  // energy(erg/cm3) = energy(code) * Units::energy
+  KOKKOS_INLINE_FUNCTION  real GetEnergy() const {return _energy;}
+  // time(s) = time(code) * Units::time
+  KOKKOS_INLINE_FUNCTION  real GetTime() const {return _time;}
+
+  KOKKOS_INLINE_FUNCTION  bool GetIsInitialized() const {return _isInitialized;}
+
+  // code-style change of the units
+  void SetLength(const real);
+  void SetVelocity(const real);
+  void SetDensity(const real);
+
+  // Show the configuration
+  void ShowConfig();
+
+ private:
+  bool _isInitialized{false};
+  // read-write variables
+  real _length{1.0};
+  real _velocity{1.0};
+  real _density{1.0};
+
+  // Deduced units from user-defined units
+  real _Kelvin{1.0};
+  real _magField{1.0};
+  real _time{1.0};
+  real _energy{1.0};
+
+  // Recompute deduced user-units
+  void ComputeUnits();
+};
+} // namespace idfx
+#endif // UNITS_HPP_