Modification to the Units module

inputs/hydro/athinput.unit_test

@@ -0,0 +1,124 @@
+<comment>
+problem    = unit test
+configure  = --prob=unit_test
+
+<job>
+problem_id = unit_test      # problem ID: basename of output filenames
+
+<output1>
+file_type  = hdf5      # Binary data dump
+variable   = prim      # variables to be output
+dt         = 1.0       # time increment between outputs
+
+<time>
+cfl_number = 0.4        # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = 0          # cycle limit
+tlim       = 4.0        # time limit
+
+<mesh>
+nx1        = 64         # Number of zones in X1-direction
+x1min      = -15.0      # minimum value of X1
+x1max      = 15.0       # maximum value of X1
+ix1_bc     = periodic   # inner-X1 boundary flag
+ox1_bc     = periodic   # outer-X1 boundary flag
+
+nx2        = 64         # Number of zones in X2-direction
+x2min      = -12.5      # minimum value of X2
+x2max      = 12.5       # maximum value of X2
+ix2_bc     = periodic   # inner-X2 boundary flag
+ox2_bc     = periodic   # outer-X2 boundary flag
+
+nx3        = 64         # Number of zones in X3-direction
+x3min      = -10.0      # minimum value of X3
+x3max      = 10.0       # maximum value of X3
+ix3_bc     = periodic   # inner-X3 boundary flag
+ox3_bc     = periodic   # outer-X3 boundary flag
+
+<hydro>
+gamma      = 1.666666666667 # gamma = C_p/C_v
+
+<problem>
+# Values of length, time, mass, ndensity, and velocity can be in
+# the units of the basis unit system.
+# All other values (e.g. energy, pressure, etc.) will have to be
+# in cgs units.
+pamb       = 1.1   # ambient pressure in **cgs** units
+energy     = 2.2   # energy value in **cgs** units
+ndamb      = 3.3   # ambient number density (ndensity) in physical units
+ndin       = 4.4   # number density (ndensity) inside radius in physical units
+radius     = 5.5   # radius of high density in physical units
+mass_val   = 6.6   # mass value in physical units
+vel        = 7.7   # velocity value in physical units
+
+<units>
+#Default value for the mean atomic weight is 1.4 
+mean_weight = 1.4
+
+#Preset unit systems
+unit_system = ism
+#unit_system = galaxy
+#unit_system = galaxypc
+#unit_system = ism_SR
+#unit_system = cgs
+#unit_system = SI
+
+#Examples of custom_basis unit systems
+#Default units are:
+#  length_unit   = pc
+#  time_unit     = Myr
+#  velocity_unit = km/s
+#  ndensity_unit = n/cm^3
+#  mass_unit     = Msun
+
+#Recreates "ism" unit system
+#unit_system = custom_basis
+#length      = 1.0    # pc
+#velocity    = 1.0    # km/s
+#ndensity    = 1.0    # n/cm^3
+
+#Recreates "galaxy" unit system
+#unit_system = custom_basis
+#length      = 1.0    # kpc
+#length_unit = kpc
+#time        = 1.0    # Myr
+#ndensity    = 1.0    # n/cm^3
+
+#Recreates "galaxypc" unit system
+#unit_system = custom_basis
+#length      = 1.0    # pc
+#time        = 1.0    # Myr
+#ndensity    = 1.0    # n/cm^3
+
+#Recreates "ism_SR" unit system
+#unit_system = custom_basis
+#length      = 1.0          # pc
+#velocity    = 2.99792458e5 # km/s
+#ndensity    = 1.0          # n/cm^3
+
+#Recreates cgs unit system
+#unit_system = custom_basis
+#length      = 1.0
+#length_unit = cm
+#time        = 1.0
+#time_unit   = s 
+#mass        = 1.0
+#mass_unit   = g
+#velocity_unit = cm/s
+
+#Recreates SI unit system
+#unit_system = custom_basis
+#length      = 1.0
+#length_unit = m
+#time        = 1.0
+#time_unit   = s 
+#mass        = 1.0
+#mass_unit   = kg
+#velocity_unit = m/s
+#ndensity_unit = n/m^3
+
+#Example of a "custom" unit system
+#Recreates "ism" unit system
+#unit_system = custom
+#mass_cgs    = 6.882618e+31   # 1 n/cm^3*pc^3
+#length_cgs  = 3.08567758e+18 # 1 pc 
+#time_cgs    = 3.085678e+13   # 1 pc/1 km/s

src/chemistry/network/gow17.cpp

@@ -318,7 +318,7 @@ ChemNetwork::ChemNetwork(MeshBlock *pmb, ParameterInput *pin) {
     // isothermal, temperature is calculated from a fixed mean molecular weight
     const Real mu_iso = pin->GetReal("chemistry", "mu_iso");
     const Real cs = pin->GetReal("hydro", "iso_sound_speed");
-    temperature_ = cs * cs * pmy_mb_->pmy_mesh->punit->code_temperature_mu_cgs * mu_iso;
+    temperature_ = cs * cs * pmy_mb_->pmy_mesh->punit->code_temp_cgs * mu_iso;
     std::cout << "isothermal temperature = " << temperature_ << " K" << std::endl;
   }
   // CR shielding
@@ -653,7 +653,7 @@ Real ChemNetwork::Edot(const Real t, const Real y[NSPECIES], const Real ED) {
     //Thermo::CoolingHotGas(nH_,  T, zdg_);
     // CO rotational lines
     // Calculate effective CO column density
-    vth = std::sqrt(2. * Constants::k_boltzmann_cgs * Tcool_nm / ChemistryUtility::mCO);
+    vth = std::sqrt(2. * Constants::k_B_cgs * Tcool_nm / ChemistryUtility::mCO);
     nCO = nH_ * yprev[iCO_];
     grad_small_ = vth/Leff_CO_max_;
     gradeff = std::max(gradv_, grad_small_);

src/chemistry/network/kida.cpp

@@ -105,7 +105,7 @@ ChemNetwork::ChemNetwork(MeshBlock *pmb, ParameterInput *pin) :
     // isothermal, temperature is calculated from a fixed mean molecular weight
     const Real mu_iso = pin->GetReal("chemistry", "mu_iso");
     const Real cs = pin->GetReal("hydro", "iso_sound_speed");
-    temperature_ = cs * cs * pmy_mb_->pmy_mesh->punit->code_temperature_mu_cgs * mu_iso;
+    temperature_ = cs * cs * pmy_mb_->pmy_mesh->punit->code_temp_cgs * mu_iso;
     std::cout << "isothermal temperature = " << temperature_ << " K" << std::endl;
   }
   // whether to cap temperature if the reaction is outside of the temperature range
@@ -698,7 +698,7 @@ Real ChemNetwork::Edot(const Real t, const Real *y, const Real ED) {
     if (ispec_map_.find("CO") != ispec_map_.end()) {
       // Calculate effective CO column density
       Real y_CO = y0[ispec_map_["CO"]];
-      vth = std::sqrt(2. * Constants::k_boltzmann_cgs * Tcool_nm / ChemistryUtility::mCO);
+      vth = std::sqrt(2. * Constants::k_B_cgs * Tcool_nm / ChemistryUtility::mCO);
       nCO = nH_ * y_CO;
       grad_small = vth/Leff_CO_max_;
       gradeff = std::max(gradv_, grad_small);
@@ -1291,7 +1291,7 @@ void ChemNetwork::InitializeReactions() {
         frml_gr_(igr) = pr->formula_;
         kgr_(igr) = 0.;
         TDgr_(igr) = pr->gamma_;
-        nu0gr_(igr) = std::sqrt( 3.0e15*pr->gamma_*Constants::k_boltzmann_cgs
+        nu0gr_(igr) = std::sqrt( 3.0e15*pr->gamma_*Constants::k_B_cgs
                                  /(M_PI*M_PI*mi) );
         igr++;
       } else {
@@ -1357,8 +1357,8 @@ void ChemNetwork::InitializeReactions() {
         }
         nu_gc_(igc) = species_[ispec_map_[pr->reactants_[0]]].charge_ / q_charge;
         r1_gc_(igc) = br*se* M_PI *ag*ag
-                      * std::sqrt(8.*Constants::k_boltzmann_cgs/(M_PI*mi));
-        t1_gc_(igc) = ag * Constants::k_boltzmann_cgs / (qi*qi);
+                      * std::sqrt(8.*Constants::k_B_cgs/(M_PI*mi));
+        t1_gc_(igc) = ag * Constants::k_B_cgs / (qi*qi);
         igc++;
       } else if (pr->formula_ == 9) { // neutral freeze-out
         const Real ag = pr->gamma_;
@@ -1379,7 +1379,7 @@ void ChemNetwork::InitializeReactions() {
         }
         kgc_(igc) = 0.;
         nu_gc_(igc) = 9; // flag for freeze-out reaction
-        r1_gc_(igc) = M_PI*ag*ag*std::sqrt( 8.*Constants::k_boltzmann_cgs/(M_PI*mi) );
+        r1_gc_(igc) = M_PI*ag*ag*std::sqrt( 8.*Constants::k_B_cgs/(M_PI*mi) );
         t1_gc_(igc) = 0.;
         igc++;
       } else {

src/chemistry/utils/kida_species.cpp

@@ -41,7 +41,7 @@ KidaSpecies::KidaSpecies(std::string line, int index) :
   for (int i=0; i < natom_; i++) {
     atom_count_[i] = std::stoi(fields[i+2]);
     mass_ += static_cast<float>(atom_count_[i]) * ma_atom_[i]
-      * Constants::hydrogen_mass_cgs;
+      * Constants::H_mass_cgs;
   }
   //electron mass
   mass_ += static_cast<float>(-charge_) * ChemistryUtility::me;

src/chemistry/utils/thermo.cpp

@@ -18,7 +18,7 @@
 
 //physical constants
 const Real Thermo::eV_ = 1.602e-12;
-const Real Thermo::kb_ = Constants::k_boltzmann_cgs;
+const Real Thermo::kb_ = Constants::k_B_cgs;
 const Real Thermo::ca_ = 2.27e-4;
 const Real Thermo::TCMB_ = 2.73;
 //ortho to para ratio of H2

src/pgen/chem_G14Sod.cpp

@@ -51,8 +51,8 @@
 void MeshBlock::ProblemGenerator(ParameterInput *pin) {
   // Define the constants
   const Real mu = pin->GetReal("problem", "mu");
-  const Real k_b = Constants::k_boltzmann_cgs;
-  const Real m_h = Constants::hydrogen_mass_cgs;
+  const Real k_b = Constants::k_B_cgs;
+  const Real m_h = Constants::H_mass_cgs;
 
   // read density and temperature
   const Real dl_cgs = pin->GetReal("problem", "LHS_rho_cgs");

src/pgen/unit_test.cpp

@@ -0,0 +1,171 @@
+//========================================================================================
+// Athena++ astrophysical MHD code
+// Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+//! \file unit_test.cpp
+//! \brief Problem file to demonstrate unit module
+//!
+
+// C headers
+
+// C++ headers
+#include <algorithm>
+#include <cmath>
+#include <cstdio>    // fopen(), fprintf(), freopen()
+#include <cstring>   // strcmp()
+#include <sstream>
+#include <stdexcept>
+#include <string>
+
+// Athena++ headers
+#include "../athena.hpp"
+#include "../athena_arrays.hpp"
+#include "../coordinates/coordinates.hpp"
+#include "../eos/eos.hpp"
+#include "../field/field.hpp"
+#include "../globals.hpp"
+#include "../hydro/hydro.hpp"
+#include "../mesh/mesh.hpp"
+#include "../parameter_input.hpp"
+
+void Mesh::InitUserMeshData(ParameterInput *pin) {
+  // Get units pointer
+  Units *punit = Mesh::punit;
+  // Print units info
+  punit->PrintBasisUnits();
+  punit->PrintCodeUnits();
+  punit->PrintConstantsInCodeUnits();
+
+  std::string lunit = std::get<1>(punit->basis_length);
+  std::string tunit = std::get<1>(punit->basis_time);
+  std::string munit = std::get<1>(punit->basis_mass);
+  std::string vunit = std::get<1>(punit->basis_velocity);
+  std::string nunit = std::get<1>(punit->basis_ndensity);
+
+  // Convert time and mesh values from input file
+  Real tlim0  = pin->GetReal("time", "tlim");
+  Real dt0    = pin->GetReal("output1", "dt");
+  Real x1max0 = pin->GetReal("mesh", "x1max");
+  punit->ConvertInputFile(pin);
+  Real tlim1  = pin->GetReal("time", "tlim");
+  Real dt1    = pin->GetReal("output1", "dt");
+  Real x1max1 = pin->GetReal("mesh", "x1max");
+
+  std::cout << "=== Values from input file converted ===" << std::endl;
+  std::cout << "From the <time> block:" << std::endl;
+  std::cout << "  tlim before conversion  = " << tlim0 << " " << tunit << std::endl;
+  std::cout << "  tlim after conversion   = " << tlim1 << " code time" << std::endl;
+  std::cout << std::endl;
+  std::cout << "From the <output1> block:" << std::endl;
+  std::cout << "  dt before conversion    = " << dt0 << " " << tunit << std::endl;
+  std::cout << "  dt after conversion     = " << dt1 << " code time" << std::endl;
+  std::cout << std::endl;
+  std::cout << "From the <mesh> block:" << std::endl;
+  std::cout << "  x1max before conversion = " << x1max0 << " " << lunit << std::endl;
+  std::cout << "  x1max after conversion  = " << x1max1 << " code length" << std::endl;
+  std::cout << std::endl;
+
+
+  // Pressure
+  Real pres = pin->GetReal("problem", "pamb");
+  Real pres_code = pres/punit->code_pressure_cgs;
+  std::cout << "Pressure value in cgs units  = " << pres << " dyne" << std::endl;
+  std::cout << "Pressure value in code units = " << pres_code
+            << " code pressure" << std::endl;
+  std::cout << std::endl;
+
+  // Energy
+  Real en = pin->GetReal("problem", "energy");
+  Real en_code = en/punit->code_energy_cgs;
+  std::cout << "Energy value in cgs units  = " << en << " erg" << std::endl;
+  std::cout << "Energy value in code units = " << en_code << " code energy" << std::endl;
+  std::cout << std::endl;
+
+  // Mass
+  Real mass = pin->GetReal("problem", "mass_val");
+  Real mass_code = mass/std::get<0>(punit->basis_mass);
+  std::cout << "Mass value in basis units = " << mass << " " << munit << std::endl;
+  std::cout << "Mass value in code units  = " << mass_code << " code mass" << std::endl;
+  std::cout << std::endl;
+
+  // Length
+  Real rad = pin->GetReal("problem", "radius");
+  Real rad_code = rad/std::get<0>(punit->basis_length);
+  std::cout << "Length value in basis units = " << rad << " " << lunit << std::endl;
+  std::cout << "Length value in code units  = " << rad_code
+            << " code length" << std::endl;
+  std::cout << std::endl;
+
+  // Velocity
+  Real vel = pin->GetReal("problem", "vel");
+  Real vel_code = vel/std::get<0>(punit->basis_velocity);
+  std::cout << "Velocity value in basis units = " << vel << " " << vunit << std::endl;
+  std::cout << "Velocity value in code units  = " << vel_code
+            << " code velocity" << std::endl;
+  std::cout << std::endl;
+
+  // ndensity
+  Real ndin = pin->GetReal("problem", "ndin");
+  Real ndin_code = ndin/std::get<0>(punit->basis_ndensity);
+  std::cout << "ndensity value in basis units = " << ndin << " " << nunit << std::endl;
+  std::cout << "ndensity value in code units  = " << ndin_code
+            << " code ndensity" << std::endl;
+  std::cout << std::endl;
+
+  std::cout << "Gravitational constant in code units = "
+            << punit->Gconst_code << std::endl;
+  return;
+}
+
+//========================================================================================
+//! \fn void MeshBlock::ProblemGenerator(ParameterInput *pin)
+//! \brief
+//========================================================================================
+
+void MeshBlock::ProblemGenerator(ParameterInput *pin) {
+  // Get units pointer
+  Units *punit = MeshBlock::pmy_mesh->punit;
+
+  Real pres = pin->GetReal("problem", "pamb")/punit->code_pressure_cgs;
+  Real da   = pin->GetReal("problem", "ndamb")/std::get<0>(punit->basis_ndensity);
+  Real din  = pin->GetReal("problem", "ndin")/std::get<0>(punit->basis_ndensity);
+  Real rin  = pin->GetReal("problem", "radius")/std::get<0>(punit->basis_length);
+  Real vel  = pin->GetReal("problem", "vel")/std::get<0>(punit->basis_velocity);
+
+  Real gamma = peos->GetGamma();
+  Real gm1 = gamma - 1.0;
+
+  // setup uniform ambient medium with spherical over-pressured region
+  for (int k=ks; k<=ke; k++) {
+    for (int j=js; j<=je; j++) {
+      for (int i=is; i<=ie; i++) {
+        Real x = pcoord->x1v(i);
+        Real y = pcoord->x2v(j);
+        Real z = pcoord->x3v(k);
+        Real rad = std::sqrt(SQR(x) + SQR(y) + SQR(z));
+        Real den = da;
+        if (rad < rin) {
+          den = din;
+        }
+        phydro->u(IDN,k,j,i) = den;
+        phydro->u(IM1,k,j,i) = vel;
+        phydro->u(IM2,k,j,i) = vel/2;
+        phydro->u(IM3,k,j,i) = vel/3;
+        phydro->u(IEN,k,j,i) = pres/gm1;
+      }
+    }
+  }
+}
+
+//========================================================================================
+//! \fn void Mesh::UserWorkAfterLoop(ParameterInput *pin)
+//! \brief Print yt units override
+//========================================================================================
+
+void Mesh::UserWorkAfterLoop(ParameterInput *pin) {
+  Units *punit = Mesh::punit;
+  // Prints a Python dictionary to use units_override
+  //  when loading hdf5 data into yt.
+  punit->PrintytUnitsOverride();
+}