Merge pull request #122 from parthenon-hpc-lab/BenWibking/fix-magneti…

…c-units

add sqrt(4pi) factor to magnetic field unit conversion

CHANGELOG.md

@@ -7,6 +7,7 @@
 - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 122]](https://github.com/parthenon-hpc-lab/athenapk/pull/122) Fixed sqrt(4pi) factor in CGS Gauss unit and add unit doc
 - [[PR 119]](https://github.com/parthenon-hpc-lab/athenapk/pull/119) Fixed Athena++ paper test case for KHI pgen. Added turbulence pgen doc.
 - [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97) Fixed Schure cooling curve. Removed SD one. Added description of cooling function conventions.
 - [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)

README.md

@@ -131,21 +131,8 @@ the `file_type = hdf5` format, see
 [VisIt](https://wci.llnl.gov/simulation/computer-codes/visit/).
 In ParaView, select the "XDMF Reader" when prompted.
 
-2. With [yt](https://yt-project.org/) -- though currently through a custom frontend
-that is not yet part of the main yt branch and, thus, has to be installed manually, e.g.,
-as follows:
-```bash
-cd ~/src # or any other folder of choice
-git clone https://github.com/forrestglines/yt.git
-cd yt
-git checkout parthenon-frontend
-
-# If you're using conda or virtualenv
-pip install -e .
-# OR alternatively, if you using the plain Python environment
-pip install --user -e .
-```
-Afterwards, `*.phdf` files can be read as usual with `yt.load()`.
+2. With [yt](https://yt-project.org/)
+As of versions >=4.4 `*.phdf` files can be read as usual with `yt.load()`.
 
 3. Using [Ascent](https://github.com/Alpine-DAV/ascent) (for in situ visualization and analysis).
 This requires Ascent to be installed/available at compile time of AthenaPK.

docs/README.md

@@ -13,6 +13,7 @@ The documentation currently includes
   - [Notebooks to calculate cooling tables from literature](cooling)
 - [Brief notes on developing code for AthenaPK](development.md)
 - [How to add a custom/user problem generator](pgen.md)
+- [Units](units.md)
 - Detailed descriptions of more complex problem generators
   - [Galaxy Cluster and Cluster-like Problem Setup](cluster.md)
   - [Driven turbulence](turbulence.md)

docs/input.md

@@ -67,6 +67,10 @@ To control the floors, following parameters can be set in the `<hydro>` block:
 *Note* the pressure floor will take precedence over the temperature floor in the
 conserved to primitive conversion if both are defined.
 
+#### Units
+
+See(here)[units.md].
+
 #### Diffusive processes
 
 Diffusive processes in AthenaPK can be configured in the `<diffusion>` block of the input file.

docs/units.md

@@ -0,0 +1,58 @@
+# AthenaPK units
+
+## General unit system
+
+Internally, all calculations are done in "code units" and there are no conversions between
+code and physical units during runtime (with excetion like the temperature when cooling is
+being used).
+Therefore, in general no units need to be prescribed to run a simulation.
+
+If units are required (e.g., if cooling is used and, thus, a conversion between internal energy
+in code units and physical temperature is required) they are configured in the input block
+as follows:
+
+```
+<units>
+code_length_cgs = 3.085677580962325e+24 # 1 Mpc in cm
+code_mass_cgs = 1.98841586e+47          # 1e14 Msun in g
+code_time_cgs = 3.15576e+16             # 1 Gyr in s
+```
+
+This information will also be used by postprocessing tools (like yt) to convert between
+code units and a physical unit system (like cgs).
+
+Moreover, internally a set of factors from code to cgs units are available to process conversions
+if required (e.g., from the input file).
+
+For example, for an input parameter (in the input file) like
+
+```
+<problem/cloud>
+r0_cgs = 3.085677580962325e+20          # 100 pc
+```
+
+the conversion should happen in the problem generator lik
+
+```c++
+  r_cloud = pin->GetReal("problem/cloud", "r0_cgs") / units.code_length_cgs();
+```
+
+so that the resulting quantity is internally in code units (here code length).
+
+It is highly recommended to be *very* explicit/specific about units everywhere (as it is
+a common source of confusion) like adding the `_cgs` suffix to the parameter in the
+input file above.
+
+## Magnetic units
+
+Internally, AthenaPK (and almost all MHD codes) use
+[Heaviside-Lorentz units](https://en.wikipedia.org/wiki/Heaviside%E2%80%93Lorentz_units),
+where the magnetic field is transformed from $B \rightarrow B / \sqrt{4 \pi}$.
+(See also the note in the
+[Castro documentation](https://amrex-astro.github.io/Castro/docs/mhd.html) about this.)
+
+So when converting from CGS-Gaussian units to code units, it is necessary to divide
+by $\sqrt{4 \pi}$ (in addition to the base dimensional factors).
+This is automatically handled by the `units.code_magnetic_cgs()` factors.
+
+

src/units.hpp

@@ -90,7 +90,8 @@ class Units {
     return code_energy_cgs() / kev_cgs * code_length_cgs() * code_length_cgs();
   }
   parthenon::Real code_magnetic_cgs() const {
-    return sqrt(code_mass_cgs()) / sqrt(code_length_cgs()) / code_time_cgs();
+    return std::sqrt(4.0 * M_PI) * sqrt(code_mass_cgs()) / sqrt(code_length_cgs()) /
+           code_time_cgs();
   }
 
   // Physical Constants in code units