plasma temperature units are now all provided in keV.

bluemira/base/constants.py

@@ -62,7 +62,11 @@ def _add_contexts(self, contexts: Optional[List[Context]] = 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)
@@ -107,6 +111,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"""
@@ -248,15 +277,14 @@ def _transform(
 # Physical constants
 # =============================================================================
 
-# Speed of light
-C_LIGHT = ureg.Quantity("c").to_base_units().magnitude  # [m/s]
-
 # Vacuum permeability
 MU_0 = ureg.Quantity("mu_0").to_base_units().magnitude  # [T.m/A] or [V.s/(A.m)]
 
 # 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)]
@@ -337,15 +365,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/codes/process/_solver.py

@@ -11,6 +11,7 @@
 
 import numpy as np
 
+from bluemira.base.constants import raw_uc
 from bluemira.base.look_and_feel import bluemira_warn
 from bluemira.base.parameter_frame import ParameterFrame
 from bluemira.codes.error import CodesError
@@ -223,7 +224,7 @@
         Returns
         -------
         tref:
-            The temperature in eV.
+            The temperature in keV.
         l_ref:
             The loss function value $L_z(n_e, T_e)$ in W.m3.
         z_ref:
@@ -234,8 +235,10 @@
         t_ref = filter(lambda lz: lz.content == "Te[eV]", lz_ref)
         lz_ref = filter(lambda lz: f"{confinement_time_ms:.1f}" in lz.content, lz_ref)
         z_av_ref = filter(lambda z: f"{confinement_time_ms:.1f}" in z.content, z_ref)
-        return tuple(
-            np.array(next(ref).data, dtype=float) for ref in (t_ref, lz_ref, z_av_ref)
+        return (
+            raw_uc(np.array(next(t_ref).data, dtype=float), "eV", "keV"),
+            np.array(next(lz_ref).data, dtype=float),
+            np.array(next(z_av_ref).data, dtype=float),
         )
 
     def get_species_fraction(self, impurity: str) -> float:

bluemira/fuel_cycle/cycle.py

@@ -20,7 +20,7 @@
 import numpy as np
 from scipy.interpolate import interp1d
 
-from bluemira.base.constants import N_AVOGADRO, T_LAMBDA, T_MOLAR_MASS, YR_TO_S, raw_uc
+from bluemira.base.constants import N_AVOGADRO, T_LAMBDA, T_MOLAR_MASS, raw_uc
 from bluemira.base.look_and_feel import bluemira_print
 from bluemira.fuel_cycle.blocks import FuelCycleComponent, FuelCycleFlow
 from bluemira.fuel_cycle.tools import (
@@ -186,10 +186,11 @@
         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 = raw_uc(dt, "yr", "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 @@
 
         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(raw_uc(self.DEMO_t[-1], "yr", "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

@@ -17,7 +17,7 @@
 import numpy as np
 from matplotlib.lines import Line2D
 
-from bluemira.base.constants import S_TO_YR, YR_TO_S, raw_uc
+from bluemira.base.constants import raw_uc
 from bluemira.base.look_and_feel import bluemira_print, bluemira_warn
 from bluemira.fuel_cycle.timeline import Timeline
 from bluemira.utilities.tools import abs_rel_difference, is_num, json_writer
@@ -125,8 +125,8 @@
 
         self.n_blk_replace = 1  # HLR
         self.n_div_replace = ndivch_in1blk + ndivch_in2blk
-        m_short = self.maintenance_s * S_TO_YR
-        m_long = self.maintenance_l * S_TO_YR
+        m_short = raw_uc(self.maintenance_s, "s", "yr")
+        m_long = raw_uc(self.maintenance_l, "s", "yr")
         phases = []
         for i in range(ndivch_in1blk):
             p_str = "Phase P1." + str(i + 1)
@@ -156,10 +156,12 @@
                 fpy += phases[i][0]
         self.fpy = fpy
         # Irreplaceable components life checks
-        self.t_on_total = self.fpy * YR_TO_S  # [s] total fusion time
+        self.t_on_total = raw_uc(self.fpy, "yr", "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(
+                raw_uc(self.params.tf_fluence / tf_ins_nflux, "s", "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"
@@ -174,7 +176,7 @@
                 f" years, or {vvlifeperc:.2f} % of neutron budget."
             )
             # TODO: treat output parameter
-        self.n_cycles = self.fpy * YR_TO_S / self.t_flattop
+        self.n_cycles = raw_uc(self.fpy, "yr", "s") / self.t_flattop
 
     def set_availabilities(self, load_factor: float):
         """
@@ -214,7 +216,7 @@
         Calculate the number of pulses per phase.
         """
         self.n_pulse_p = [
-            int(YR_TO_S * phases[i][0] // self.t_flattop)
+            int(raw_uc(phases[i][0], "yr", "s") // self.t_flattop)
             for i in range(len(phases))
             if phases[i][1].startswith("Phase P")
         ]
@@ -276,12 +278,14 @@
         results that violate the tolerances.
         """
         life = self.fpy / self.params.A_global
-        actual_life = S_TO_YR * (
+        actual_life = raw_uc(
             self.t_on_total
             + self.total_ramptime
             + self.t_interdown
             + self.total_planned_maintenance
-            + self.t_unplanned_m
+            + self.t_unplanned_m,
+            "s",
+            "yr",
         )
         actual_lf = self.fpy / actual_life
         delt = abs_rel_difference(actual_life, life)
@@ -318,7 +322,9 @@
                 self.inputs,
             )  # Phoenix
 
-        if self.params.A_global > self.fpy / (self.fpy + S_TO_YR * self.min_downtime):
+        if self.params.A_global > self.fpy / (
+            self.fpy + raw_uc(self.min_downtime, "s", "yr")
+        ):
             bluemira_warn("FuelCycle::Lifecyle: Input availability is unachievable.")
         # Re-assign A
         self.params.A_global = actual_lf

bluemira/fuel_cycle/timeline.py

@@ -13,7 +13,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from bluemira.base.constants import S_TO_YR, YR_TO_S
+from bluemira.base.constants import raw_uc
 from bluemira.fuel_cycle.timeline_tools import (
     LogNormalAvailabilityStrategy,
     OperationalAvailabilityStrategy,
@@ -354,7 +354,7 @@
                 )
                 j += 1
             elif "Phase M" in name:
-                p = MaintenancePhase(name, duration * YR_TO_S, t_start=t_start)
+                p = MaintenancePhase(name, raw_uc(duration, "yr", "s"), t_start=t_start)
             phases.append(p)
         self.phases = phases
         self.build_arrays(phases)
@@ -385,9 +385,8 @@
         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.t *= S_TO_YR
+        self.ft = raw_uc(np.cumsum(self.ft), "s", "yr")
+        self.t = raw_uc(self.t, "s", "yr")
         self.plant_life = self.t[-1]  # total plant lifetime [calendar]
 
     def to_dict(self) -> Dict[str, Union[np.ndarray, int]]:

bluemira/fuel_cycle/tools.py

@@ -263,8 +263,8 @@
         The time vector
     m_t_flow:
         The mass flow vector
-    t_delay:
-        The delay duration [s]
+    tt_delay:
+        The delay duration [yr]
 
     Returns
     -------
@@ -575,7 +575,7 @@
     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
@@ -731,13 +731,12 @@
     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
@@ -825,13 +824,12 @@
     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 @@
     -------
     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 @@
     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 @@
     -------
     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 @@
     Z_eff:
         Effective charge [a.m.u.]
     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 @@
 
     \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)
+    )