Fix units and unit-conversion factors (#2891)

* Removed redundant units

* Removed all 'years' variable and YR_TO_S variables. Found two possible errors.

* Made units more understandable

* Undid all the raw_uc conversions in tools.py as they are not compatible with numba.jit's optimization. Using the constants YR_TO_S and S_TO_YR.

* Removed as many S_TO_YR and YR_TO_S as possible.

* Fixed the problem of missing elementary charge (e) constant in base.constants

* Undid all S_TO_YR/YR_TO_S changes

* Minor fixes to the YR_TO_S

* Added back in C_LIGHT

bluemira/base/constants.py

@@ -92,7 +92,11 @@ def _add_contexts(self, contexts: list[Context] | None = None):
         Add new contexts to registry
         """
         if not self._contexts_added:
-            self.contexts = [self._energy_temperature_context(), self._flow_context()]
+            self.contexts = [
+                self._energy_temperature_context(),
+                self._mass_energy_context(),
+                self._flow_context(),
+            ]
 
             for c in self.contexts:
                 self.add_context(c)
@@ -137,6 +141,31 @@ def _energy_temperature_context(self):
             lambda _, x: x / conversion,
         )
 
+    def _mass_energy_context(self):
+        """
+        Converter between mass and energy
+
+        energy = mass * speed-of-light^2
+
+        Returns
+        -------
+        pint context
+        """
+        m_to_e = Context("Mass_to_Energy")
+
+        m_units = "[mass]"
+        e_units = "[energy]"
+
+        conversion = self.Quantity("c^2")
+
+        return self._transform(
+            m_to_e,
+            m_units,
+            e_units,
+            lambda _, x: x * conversion,
+            lambda _, x: x / conversion,
+        )
+
     @property
     def flow_conversion(self):
         """Gas flow conversion factor R * T"""
@@ -283,6 +312,8 @@ def _transform(
 # Vacuum permittivity
 EPS_0 = ureg.Quantity("eps_0").to_base_units().magnitude  # [A^2.s^4/kg/m^3]
 
+# absolute charge of an electron
+ELEMENTARY_CHARGE = ureg.Quantity("e").to_base_units().magnitude  # [e]
 
 # Commonly used..
 MU_0_2PI = 2e-7  # [T.m/A] or [V.s/(A.m)]
@@ -366,15 +397,6 @@ def _transform(
 # Conversions
 # =============================================================================
 
-# Electron-volts to Joules
-EV_TO_J = ureg.Quantity(1, ureg.eV).to(ureg.joule).magnitude
-
-# Joules to Electron-volts
-J_TO_EV = ureg.Quantity(1, ureg.joule).to(ureg.eV).magnitude
-
-# Atomic mass units to kilograms
-AMU_TO_KG = ureg.Quantity(1, ureg.amu).to(ureg.kg).magnitude
-
 # Years to seconds
 YR_TO_S = ureg.Quantity(1, ureg.year).to(ureg.second).magnitude
 

bluemira/fuel_cycle/cycle.py

@@ -186,10 +186,11 @@ def tbreed(self, TBR: float, m_T_0: float) -> np.ndarray:
         m_T = m_T_0 * np.ones(len(self.DEMO_t))
         for i in range(1, len(self.DEMO_t)):
             dt = self.DEMO_t[i] - self.DEMO_t[i - 1]
-            t_bred = TBR * self.brate[i] * (YR_TO_S * dt)
-            t_bred += self.prate[i] * (YR_TO_S * dt)
-            t_burnt = self.brate[i] * (YR_TO_S * dt)
-            t_DD = self.prate[i] * (YR_TO_S * dt)
+            dts = dt * YR_TO_S
+            t_bred = TBR * self.brate[i] * dts
+            t_bred += self.prate[i] * dts
+            t_burnt = self.brate[i] * dts
+            t_DD = self.prate[i] * dts
             m_T[i] = (m_T[i - 1]) * np.exp(-T_LAMBDA * dt) - t_burnt + t_bred + t_DD
         return m_T
 
@@ -350,7 +351,7 @@ def recycle(self):
 
         m_T_out = self.plasma(self.params.eta_iv, self.params.I_miv, flows=flows)
         # Resolution - Not used everywhere for speed
-        n_ts = int(round((YR_TO_S * self.DEMO_t[-1]) / self.timestep))
+        n_ts = int(round(self.DEMO_t[-1] * YR_TO_S / self.timestep))
         self.t, m_pellet_in = discretise_1d(self.DEMO_t, self.m_T_in, n_ts)
 
         # Flow out of the vacuum vessel

bluemira/fuel_cycle/lifecycle.py

@@ -175,7 +175,7 @@ def life_neutronics(self):
         self.t_on_total = self.fpy * YR_TO_S  # [s] total fusion time
         tf_ins_life_dose = tf_ins_nflux * self.t_on_total / self.params.tf_fluence
         if tf_ins_life_dose > 1:
-            self.tf_lifeend = round(self.params.tf_fluence / tf_ins_nflux / YR_TO_S, 2)
+            self.tf_lifeend = round(self.params.tf_fluence / tf_ins_nflux * S_TO_YR, 2)
             tflifeperc = round(100 * self.tf_lifeend / self.fpy, 1)
             bluemira_warn(
                 f"TF coil insulation fried after {self.tf_lifeend:.2f} full-power years"

bluemira/fuel_cycle/timeline.py

@@ -386,8 +386,7 @@ def concatenate(p_phases, k_key):
         self.ft = np.zeros(len(self.t))
         for i in fuse_indices[1::2]:
             self.ft[i] = self.t[i] - self.t[i - 1]
-        self.ft = np.cumsum(self.ft)
-        self.ft *= S_TO_YR
+        self.ft = np.cumsum(self.ft) * S_TO_YR
         self.t *= S_TO_YR
         self.plant_life = self.t[-1]  # total plant lifetime [calendar]
 

bluemira/fuel_cycle/tools.py

@@ -282,8 +282,8 @@ def delay_decay(t: np.ndarray, m_t_flow: np.ndarray, tt_delay: float) -> np.ndar
         The time vector
     m_t_flow:
         The mass flow vector
-    t_delay:
-        The delay duration [s]
+    tt_delay:
+        The delay duration [yr]
 
     Returns
     -------
@@ -594,7 +594,7 @@ def _fountain_linear_sink(
     if dt == 0:
         return m_flow, inventory, sum_in, decayed
 
-    m_in = m_flow * YR_TO_S  # kg/yr
+    m_in = m_flow * YR_TO_S  # converts kg/s to kg/yr
     dts = dt * YR_TO_S
     mass_in = m_flow * dts
     sum_in += mass_in
@@ -747,13 +747,12 @@ def _linear_thresh_sink(
     decayed:
         Accountancy parameter to calculate the total value of decayed T in a sink
     """
-    years = 365 * 24 * 3600
     dt = t_out - t_in
     if dt == 0:
         return m_flow, inventory, sum_in, decayed
 
-    m_in = m_flow * years  # kg/yr
-    dts = dt * years
+    m_in = m_flow * YR_TO_S  # converts kg/s to kg/yr
+    dts = dt * YR_TO_S
     mass_in = m_flow * dts
     sum_in += mass_in
     j_inv0 = inventory
@@ -842,13 +841,12 @@ def _sqrt_thresh_sink(
     not accounted for in this function. We have to add decay in the sink and
     ensure this is handled when calculation the absorbtion and out-flow.
     """
-    years = 365 * 24 * 3600
     dt = t_out - t_in
     if dt == 0:
         # Nothing can happen if time is zero
         return m_flow, inventory, sum_in, decayed
 
-    dts = dt * years
+    dts = dt * YR_TO_S
     mass_in = m_flow * dts
     sum_in += mass_in
 

bluemira/plasma_physics/collisions.py

@@ -12,11 +12,11 @@
 
 from bluemira.base.constants import (
     ELECTRON_MASS,
+    ELEMENTARY_CHARGE,
     EPS_0,
-    EV_TO_J,
     H_PLANCK,
-    K_BOLTZMANN,
     PROTON_MASS,
+    raw_uc,
 )
 
 
@@ -35,7 +35,9 @@ def debye_length(temperature: float, density: float) -> float:
     -------
     Debye length [m]
     """
-    return np.sqrt(EPS_0 * K_BOLTZMANN * temperature / (EV_TO_J**2 * density))
+    return np.sqrt(
+        EPS_0 * raw_uc(temperature, "K", "J") / (ELEMENTARY_CHARGE**2 * density)
+    )
 
 
 def reduced_mass(mass_1: float, mass_2: float) -> float:
@@ -74,7 +76,10 @@ def thermal_velocity(temperature: float, mass: float) -> float:
     The sqrt(2) term is for a 3-dimensional system and the most probable velocity in
     the particle velocity distribution.
     """
-    return np.sqrt(2) * np.sqrt(K_BOLTZMANN * temperature / mass)
+    return np.sqrt(2) * np.sqrt(
+        raw_uc(temperature, "K", "J")  # = Joule = kg*m^2/s^2
+        / mass
+    )  # sqrt(m^2/s^2) = m/s
 
 
 def de_broglie_length(velocity: float, mu_12: float) -> float:
@@ -110,7 +115,7 @@ def impact_parameter_perp(velocity: float, mu_12: float) -> float:
     -------
     Perpendicular impact parameter [m]
     """
-    return EV_TO_J**2 / (4 * np.pi * EPS_0 * mu_12 * velocity**2)
+    return ELEMENTARY_CHARGE**2 / (4 * np.pi * EPS_0 * mu_12 * velocity**2)
 
 
 def coulomb_logarithm(temperature: float, density: float) -> float:
@@ -146,7 +151,8 @@ def spitzer_conductivity(Z_eff: float, T_e: float, ln_lambda: float) -> float:
     Z_eff:
         Effective charge [dimensionless]
     T_e:
-        Electron temperature on axis [eV]
+        Electron temperature on axis [keV]
+        The equation takes in temperature as [eV], so an in-line conversion is used here.
     ln_lambda:
         Coulomb logarithm value
 
@@ -160,4 +166,9 @@ def spitzer_conductivity(Z_eff: float, T_e: float, ln_lambda: float) -> float:
 
     \t:math:`\\sigma = 1.92e4 (2-Z_{eff}^{-1/3}) \\dfrac{T_{e}^{3/2}}{Z_{eff}ln\\Lambda}`
     """
-    return 1.92e4 * (2 - Z_eff ** (-1 / 3)) * T_e**1.5 / (Z_eff * ln_lambda)
+    return (
+        1.92e4
+        * (2 - Z_eff ** (-1 / 3))
+        * raw_uc(T_e, "keV", "eV") ** 1.5
+        / (Z_eff * ln_lambda)
+    )

bluemira/plasma_physics/reactions.py

@@ -15,14 +15,10 @@
 import numpy.typing as npt
 
 from bluemira.base.constants import (
-    AMU_TO_KG,
-    C_LIGHT,
     D_MOLAR_MASS,
     ELECTRON_MOLAR_MASS,
-    EV_TO_J,
     HE3_MOLAR_MASS,
     HE_MOLAR_MASS,
-    J_TO_EV,
     NEUTRON_MOLAR_MASS,
     N_AVOGADRO,
     PROTON_MOLAR_MASS,
@@ -58,7 +54,7 @@ def E_DT_fusion() -> float:
         \\Delta E = \\Delta m c^2
     """
     delta_m = (D_MOLAR_MASS + T_MOLAR_MASS) - (HE_MOLAR_MASS + NEUTRON_MOLAR_MASS)
-    return delta_m * C_LIGHT**2 * AMU_TO_KG * J_TO_EV
+    return raw_uc(delta_m, "amu", "eV")
 
 
 def E_DD_fusion() -> float:
@@ -85,7 +81,7 @@ def E_DD_fusion() -> float:
         (D_MOLAR_MASS + D_MOLAR_MASS) - (HE3_MOLAR_MASS + NEUTRON_MOLAR_MASS),
     ])
     delta_m = np.average(delta_m)
-    return delta_m * C_LIGHT**2 * AMU_TO_KG * J_TO_EV
+    return raw_uc(delta_m, "amu", "eV")
 
 
 def n_DT_reactions(p_fus: float) -> float:
@@ -105,7 +101,7 @@ def n_DT_reactions(p_fus: float) -> float:
     Number of D-T reactions per second [1/s]
     """
     e_dt = E_DT_fusion()
-    return raw_uc(p_fus, "MW", "W") / (e_dt * EV_TO_J)
+    return raw_uc(p_fus, "MW", "W") / raw_uc(e_dt, "eV", "J")
 
 
 def n_DD_reactions(p_fus: float) -> float:
@@ -125,7 +121,7 @@ def n_DD_reactions(p_fus: float) -> float:
     Number of D-D reactions per second [1/s]
     """
     e_dd = E_DD_fusion()
-    return p_fus / (e_dd * EV_TO_J)
+    return p_fus / raw_uc(e_dd, "eV", "J")
 
 
 def r_T_burn(p_fus: float) -> float:  # noqa: N802

bluemira/plasma_physics/rules_of_thumb.py

@@ -10,7 +10,7 @@
 
 import numpy as np
 
-from bluemira.base.constants import EV_TO_J, K_BOLTZMANN, MU_0
+from bluemira.base.constants import MU_0, raw_uc
 from bluemira.plasma_physics.collisions import coulomb_logarithm, spitzer_conductivity
 
 
@@ -29,7 +29,7 @@ def estimate_loop_voltage(
     Z_eff:
         Effective charge [dimensionless]
     T_e:
-        Electron temperature on axis [eV]
+        Electron temperature on axis [keV]
     n_e:
         Electron density [1/m^3]
     q_0:
@@ -54,7 +54,7 @@ def estimate_loop_voltage(
 
     There is no neo-classical resistivity on axis because there are no trapped particles
     """  # noqa: W505, E501
-    ln_lambda = coulomb_logarithm(T_e * EV_TO_J / K_BOLTZMANN, n_e)
+    ln_lambda = coulomb_logarithm(raw_uc(T_e, "keV", "K"), n_e)
     sigma = spitzer_conductivity(Z_eff, T_e, ln_lambda)
 
     # Current density on axis