diff --git a/doc/source/reference/idefix.ini.rst b/doc/source/reference/idefix.ini.rst index ae2b8e5a9..1921b4265 100644 --- a/doc/source/reference/idefix.ini.rst +++ b/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 ------------------ diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 727f21504..08dd93c61 100644 --- a/src/CMakeLists.txt +++ b/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 ) diff --git a/src/global.cpp b/src/global.cpp index 1e9c11eef..806717978 100644 --- a/src/global.cpp +++ b/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; diff --git a/src/global.hpp b/src/global.hpp index 47fca42ee..33478f728 100644 --- a/src/global.hpp +++ b/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(); diff --git a/src/gravity/gravity.cpp b/src/gravity/gravity.cpp index 2492ce090..ae4764047 100644 --- a/src/gravity/gravity.cpp +++ b/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("Gravity","gravCst",0, 1.0); + // should we compute it from the units? + if(input.CheckEntry("Gravity","gravCst")>=0) { + if(input.Get("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("Gravity","gravCst",0); + } + } else { // default value to 1.0 + this->gravCst = 1.0; + } + // Gravitational potential int nPotential = input.CheckEntry("Gravity","potential"); if(nPotential >=0) { diff --git a/src/main.cpp b/src/main.cpp index 0346b1c60..5a24043f4 100644 --- a/src/main.cpp +++ b/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(); diff --git a/src/units.cpp b/src/units.cpp new file mode 100644 index 000000000..2baf97b1e --- /dev/null +++ b/src/units.cpp @@ -0,0 +1,55 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// Copyright(C) Geoffroy R. J. Lesur +// 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("Units","length",0,1.0); + this->_velocity = input.GetOrSet("Units","velocity",0,1.0); + this->_density = input.GetOrSet("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; +} diff --git a/src/units.hpp b/src/units.hpp new file mode 100644 index 000000000..ffafdd4f1 --- /dev/null +++ b/src/units.hpp @@ -0,0 +1,87 @@ +// *********************************************************************************** +// Idefix MHD astrophysical code +// Copyright(C) Geoffroy R. J. Lesur +// 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_