From 4308f6d4d74df47d8416f584a32063754d3beb0a Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Fri, 14 Jun 2024 15:58:48 +0200 Subject: [PATCH 01/11] add unit system to idefix --- src/CMakeLists.txt | 2 ++ src/global.cpp | 2 ++ src/global.hpp | 2 ++ src/main.cpp | 5 ++++ src/units.cpp | 49 +++++++++++++++++++++++++++++++ src/units.hpp | 73 ++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 133 insertions(+) create mode 100644 src/units.cpp create mode 100644 src/units.hpp 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/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..3d9eb2084 --- /dev/null +++ b/src/units.cpp @@ -0,0 +1,49 @@ +// *********************************************************************************** +// 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) { + 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(is_initialized) { + 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: [Temperature] = " << this->_Kelvin << " K" << std::endl; + idfx::cout << "Units: [Mag. Field] = " << this->_magField << " G" << std::endl; + } +} + +void idfx::Units::ComputeUnits() { + this->_is_initialized = true; + this->_magField = std::sqrt(4*M_PI*density*velocity*velocity); + this->_Kelvin = u*velocity*velocity/k_B; +} diff --git a/src/units.hpp b/src/units.hpp new file mode 100644 index 000000000..ad9f0fd2f --- /dev/null +++ b/src/units.hpp @@ -0,0 +1,73 @@ +// *********************************************************************************** +// 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 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.49598e+13}; // Astronomical unit (cm) + + + // User-defined units, non user-modifiable + const real &length{_length}; // L (cm) = L (code) * Units::length + const real &velocity{_velocity}; // V(cm/s) = V(code) * Units::velocity + const real &density{_density}; // density (g/cm^3) = density(code) + // * Units::density + + // Deduced units from user-defined units + const real &Kelvin{_Kelvin}; // T(K) = P(code)/rho(code) * mu * Units::Kelvin + const real &magField{_magField}; // B(G) = B(code) * Units::MagField + + bool &is_initialized{_is_initialized}; + + // 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 _is_initialized{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}; + + // Recompute deduced user-units + void ComputeUnits(); +}; +} // namespace idfx +#endif // UNITS_HPP_ From 155107e7c5f04554ad18c9ab074af11bdc1d4f83 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 19 Dec 2024 17:39:04 +0100 Subject: [PATCH 02/11] add proper copy constructor to units (tentative) --- src/units.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/units.hpp b/src/units.hpp index ad9f0fd2f..bbad738df 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -16,7 +16,24 @@ class Input; namespace idfx { class Units { public: + Units() = default; + // Copy constructor make reference to new object + Units(const Units &u) : length(_length), + velocity(_velocity), + density(_density), + Kelvin(_Kelvin), + magField(_magField), + is_initialized(_is_initialized), + _length(u._length), + _velocity(u._velocity), + _density(u._density), + _Kelvin(u._Kelvin), + _magField(u._Kelvin), + _is_initialized(u._is_initialized) + {} + 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) From 3f02fdefe4f701b913733d537240fcd105934124 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 19 Dec 2024 17:42:32 +0100 Subject: [PATCH 03/11] add decoration for Kokkos --- src/units.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/units.hpp b/src/units.hpp index bbad738df..dd8686e4f 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -18,7 +18,7 @@ class Units { public: Units() = default; // Copy constructor make reference to new object - Units(const Units &u) : length(_length), + KOKKOS_INLINE_FUNCTION (const Units &u) : length(_length), velocity(_velocity), density(_density), Kelvin(_Kelvin), From 557505df25c0f867d1e53c819b58aaadbf24dbc6 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 19 Dec 2024 17:44:22 +0100 Subject: [PATCH 04/11] better decoration --- src/units.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/units.hpp b/src/units.hpp index dd8686e4f..4b7a73060 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -18,7 +18,7 @@ class Units { public: Units() = default; // Copy constructor make reference to new object - KOKKOS_INLINE_FUNCTION (const Units &u) : length(_length), + KOKKOS_FUNCTION (const Units &u) : length(_length), velocity(_velocity), density(_density), Kelvin(_Kelvin), From c743f90e262f55376c63db6cbe8f2f59e0c08e1f Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Thu, 19 Dec 2024 18:09:18 +0100 Subject: [PATCH 05/11] use accessors in place of references --- src/units.cpp | 6 +++--- src/units.hpp | 28 ++++++---------------------- 2 files changed, 9 insertions(+), 25 deletions(-) diff --git a/src/units.cpp b/src/units.cpp index 3d9eb2084..90a8de5f5 100644 --- a/src/units.cpp +++ b/src/units.cpp @@ -33,7 +33,7 @@ void idfx::Units::SetVelocity(const real in) { } void idfx::Units::ShowConfig() { - if(is_initialized) { + if(_is_initialized) { 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; @@ -44,6 +44,6 @@ void idfx::Units::ShowConfig() { void idfx::Units::ComputeUnits() { this->_is_initialized = true; - this->_magField = std::sqrt(4*M_PI*density*velocity*velocity); - this->_Kelvin = u*velocity*velocity/k_B; + this->_magField = std::sqrt(4*M_PI*_density*_velocity*_velocity); + this->_Kelvin = u * _velocity*_velocity/k_B; } diff --git a/src/units.hpp b/src/units.hpp index 4b7a73060..67f16d083 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -16,22 +16,6 @@ class Input; namespace idfx { class Units { public: - Units() = default; - // Copy constructor make reference to new object - KOKKOS_FUNCTION (const Units &u) : length(_length), - velocity(_velocity), - density(_density), - Kelvin(_Kelvin), - magField(_magField), - is_initialized(_is_initialized), - _length(u._length), - _velocity(u._velocity), - _density(u._density), - _Kelvin(u._Kelvin), - _magField(u._Kelvin), - _is_initialized(u._is_initialized) - {} - void Init(Input &input); @@ -53,16 +37,16 @@ class Units { // User-defined units, non user-modifiable - const real &length{_length}; // L (cm) = L (code) * Units::length - const real &velocity{_velocity}; // V(cm/s) = V(code) * Units::velocity - const real &density{_density}; // density (g/cm^3) = density(code) + KOKKOS_INLINE_FUNCTION real length() const {return _density;} // L (cm) = L (code) * Units::length + KOKKOS_INLINE_FUNCTION real velocity() const {return _velocity;} // V(cm/s) = V(code) * Units::velocity + KOKKOS_INLINE_FUNCTION real density() const {return _density;} // density (g/cm^3) = density(code) // * Units::density // Deduced units from user-defined units - const real &Kelvin{_Kelvin}; // T(K) = P(code)/rho(code) * mu * Units::Kelvin - const real &magField{_magField}; // B(G) = B(code) * Units::MagField + KOKKOS_INLINE_FUNCTION real Kelvin() const {return _Kelvin;} // T(K) = P(code)/rho(code) * mu * Units::Kelvin + KOKKOS_INLINE_FUNCTION real magField() const {return _magField;} // B(G) = B(code) * Units::MagField - bool &is_initialized{_is_initialized}; + KOKKOS_INLINE_FUNCTION bool is_initialized() const {return _is_initialized;} // code-style change of the units void SetLength(const real); From ba5eac839d2c9c53e5ca26fe28cc775595cec729 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Fri, 20 Dec 2024 08:55:54 +0100 Subject: [PATCH 06/11] use proper getter and conform to idefix coding style --- src/units.cpp | 4 ++-- src/units.hpp | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/units.cpp b/src/units.cpp index 90a8de5f5..28879deee 100644 --- a/src/units.cpp +++ b/src/units.cpp @@ -33,7 +33,7 @@ void idfx::Units::SetVelocity(const real in) { } void idfx::Units::ShowConfig() { - if(_is_initialized) { + 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; @@ -43,7 +43,7 @@ void idfx::Units::ShowConfig() { } void idfx::Units::ComputeUnits() { - this->_is_initialized = true; + this->_isInitialized = true; this->_magField = std::sqrt(4*M_PI*_density*_velocity*_velocity); this->_Kelvin = u * _velocity*_velocity/k_B; } diff --git a/src/units.hpp b/src/units.hpp index 67f16d083..48502eb5d 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -37,16 +37,16 @@ class Units { // User-defined units, non user-modifiable - KOKKOS_INLINE_FUNCTION real length() const {return _density;} // L (cm) = L (code) * Units::length - KOKKOS_INLINE_FUNCTION real velocity() const {return _velocity;} // V(cm/s) = V(code) * Units::velocity - KOKKOS_INLINE_FUNCTION real density() const {return _density;} // density (g/cm^3) = density(code) + KOKKOS_INLINE_FUNCTION real GetLength() const {return _density;} // L (cm) = L (code) * Units::length + KOKKOS_INLINE_FUNCTION real GetVelocity() const {return _velocity;} // V(cm/s) = V(code) * Units::velocity + KOKKOS_INLINE_FUNCTION real GetDensity() const {return _density;} // density (g/cm^3) = density(code) // * Units::density // Deduced units from user-defined units - KOKKOS_INLINE_FUNCTION real Kelvin() const {return _Kelvin;} // T(K) = P(code)/rho(code) * mu * Units::Kelvin - KOKKOS_INLINE_FUNCTION real magField() const {return _magField;} // B(G) = B(code) * Units::MagField + KOKKOS_INLINE_FUNCTION real GetKelvin() const {return _Kelvin;} // T(K) = P(code)/rho(code) * mu * Units::Kelvin + KOKKOS_INLINE_FUNCTION real GetMagField() const {return _magField;} // B(G) = B(code) * Units::MagField - KOKKOS_INLINE_FUNCTION bool is_initialized() const {return _is_initialized;} + KOKKOS_INLINE_FUNCTION bool GetIsInitialized() const {return _isInitialized;} // code-style change of the units void SetLength(const real); @@ -57,7 +57,7 @@ class Units { void ShowConfig(); private: - bool _is_initialized{false}; + bool _isInitialized{false}; // read-write variables real _length{1.0}; real _velocity{1.0}; From 15e8c98c16b4a771a394e25b34f56d205d82f3a8 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 23 Apr 2025 11:39:56 +0200 Subject: [PATCH 07/11] add documentation for units use units in gravity module --- doc/source/reference/idefix.ini.rst | 25 ++++++++++++++++++++++++- src/gravity/gravity.cpp | 15 ++++++++++++++- 2 files changed, 38 insertions(+), 2 deletions(-) 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/gravity/gravity.cpp b/src/gravity/gravity.cpp index 2492ce090..5d7eb28c3 100644 --- a/src/gravity/gravity.cpp +++ b/src/gravity/gravity.cpp @@ -17,7 +17,20 @@ Gravity::Gravity(Input &input, DataBlock *datain) { 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").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) { From 5440af32e561ff9ce9685bd13c7b51bae47ecddd Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 23 Apr 2025 15:02:54 +0200 Subject: [PATCH 08/11] initialise units on request --- src/units.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/units.cpp b/src/units.cpp index 28879deee..01011a23e 100644 --- a/src/units.cpp +++ b/src/units.cpp @@ -11,10 +11,12 @@ void idfx::Units::Init(Input &input) { - 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(); + 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) { From c9891133cb071c341d30423d0613d8b11b29dfb9 Mon Sep 17 00:00:00 2001 From: nscepi <63753204+nscepi@users.noreply.github.com> Date: Wed, 23 Apr 2025 15:03:51 +0200 Subject: [PATCH 09/11] Add units for energy and time and fixed typo (#304) --- src/units.cpp | 5 +++++ src/units.hpp | 12 +++++++++--- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/src/units.cpp b/src/units.cpp index 01011a23e..32bd320ee 100644 --- a/src/units.cpp +++ b/src/units.cpp @@ -39,6 +39,8 @@ void idfx::Units::ShowConfig() { 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; } @@ -48,4 +50,7 @@ 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 index 48502eb5d..b6c9b4a8d 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -19,12 +19,13 @@ class Units { void Init(Input &input); - const real u {1.6605390666e-24}; // Atomic mass unit (g) + 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) @@ -33,11 +34,11 @@ class Units { 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.49598e+13}; // Astronomical unit (cm) + const real au{1.49597892e13}; // Astronomical unit (cm) // User-defined units, non user-modifiable - KOKKOS_INLINE_FUNCTION real GetLength() const {return _density;} // L (cm) = L (code) * Units::length + KOKKOS_INLINE_FUNCTION real GetLength() const {return _length;} // L (cm) = L (code) * Units::length KOKKOS_INLINE_FUNCTION real GetVelocity() const {return _velocity;} // V(cm/s) = V(code) * Units::velocity KOKKOS_INLINE_FUNCTION real GetDensity() const {return _density;} // density (g/cm^3) = density(code) // * Units::density @@ -45,6 +46,8 @@ class Units { // Deduced units from user-defined units KOKKOS_INLINE_FUNCTION real GetKelvin() const {return _Kelvin;} // T(K) = P(code)/rho(code) * mu * Units::Kelvin KOKKOS_INLINE_FUNCTION real GetMagField() const {return _magField;} // B(G) = B(code) * Units::MagField + KOKKOS_INLINE_FUNCTION real GetEnergy() const {return _energy;} // energy(erg/cm3) = energy(code) * Units::energy + KOKKOS_INLINE_FUNCTION real GetTime() const {return _time;} // time(s) = time(code) * Units::time KOKKOS_INLINE_FUNCTION bool GetIsInitialized() const {return _isInitialized;} @@ -66,9 +69,12 @@ class Units { // 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_ + From 4729ad203dc3cc7fc18d31e691fef8c3a973c5b2 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 23 Apr 2025 16:32:23 +0200 Subject: [PATCH 10/11] fix linter --- src/units.cpp | 1 - src/units.hpp | 31 +++++++++++++++++++------------ 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/src/units.cpp b/src/units.cpp index 32bd320ee..2baf97b1e 100644 --- a/src/units.cpp +++ b/src/units.cpp @@ -53,4 +53,3 @@ void idfx::Units::ComputeUnits() { this->_energy = _density*_velocity*_velocity; this->_time = _length/_velocity; } - diff --git a/src/units.hpp b/src/units.hpp index b6c9b4a8d..ffafdd4f1 100644 --- a/src/units.hpp +++ b/src/units.hpp @@ -17,7 +17,7 @@ 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) @@ -25,7 +25,8 @@ class Units { 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 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) @@ -34,20 +35,27 @@ class Units { 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 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 - KOKKOS_INLINE_FUNCTION real GetLength() const {return _length;} // L (cm) = L (code) * Units::length - KOKKOS_INLINE_FUNCTION real GetVelocity() const {return _velocity;} // V(cm/s) = V(code) * Units::velocity - KOKKOS_INLINE_FUNCTION real GetDensity() const {return _density;} // density (g/cm^3) = density(code) - // * Units::density + // 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 - KOKKOS_INLINE_FUNCTION real GetKelvin() const {return _Kelvin;} // T(K) = P(code)/rho(code) * mu * Units::Kelvin - KOKKOS_INLINE_FUNCTION real GetMagField() const {return _magField;} // B(G) = B(code) * Units::MagField - KOKKOS_INLINE_FUNCTION real GetEnergy() const {return _energy;} // energy(erg/cm3) = energy(code) * Units::energy - KOKKOS_INLINE_FUNCTION real GetTime() const {return _time;} // time(s) = time(code) * Units::time + // 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;} @@ -77,4 +85,3 @@ class Units { }; } // namespace idfx #endif // UNITS_HPP_ - From 5c577459c9e7652f8b4e629e587e435d97230433 Mon Sep 17 00:00:00 2001 From: Geoffroy Lesur Date: Wed, 23 Apr 2025 16:52:57 +0200 Subject: [PATCH 11/11] fix foward declaration --- src/gravity/gravity.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/gravity/gravity.cpp b/src/gravity/gravity.cpp index 5d7eb28c3..ae4764047 100644 --- a/src/gravity/gravity.cpp +++ b/src/gravity/gravity.cpp @@ -11,6 +11,7 @@ #include "planetarySystem.hpp" #include "dataBlock.hpp" #include "input.hpp" +#include "units.hpp" Gravity::Gravity(Input &input, DataBlock *datain) { idfx::pushRegion("Gravity::Gravity"); @@ -19,7 +20,7 @@ Gravity::Gravity(Input &input, DataBlock *datain) { // Gravitational constant G // should we compute it from the units? if(input.CheckEntry("Gravity","gravCst")>=0) { - if(input.Get("Gravity","gravCst").compare("units") == 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()