Add stored thermal energy calcs

documentation/source/development/standards.md

@@ -350,6 +350,12 @@ This should be used for units of $\text{kg} \cdot \text{m}^{-2}\text{s}^{-1}$
 
 ---------------------
 
+##### Energy densities
+
+- Energy densities should start with the `eden_` prefix
+
+---------------------
+
 ##### Particle fluence
 
 - Fluences should start with the `flu_` prefix

process/data_structure/physics_variables.py

@@ -182,6 +182,24 @@ class PhysicsData:
     e_plasma_beta_thermal: float = 0.0
     """Plasma thermal energy derived from thermal beta"""
 
+    e_plasma_thermal_total: float = 0.0
+    """Plasma total stored thermal energy (J)"""
+
+    eden_plasma_thermal_vol_avg: float = 0.0
+    """Plasma volume averaged thermal energy density (J/m³)"""
+
+    e_plasma_electrons_thermal: float = 0.0
+    """Plasma thermal energy in electrons (J)"""
+
+    eden_plasma_electrons_thermal_vol_avg: float = 0.0
+    """Plasma volume averaged thermal energy density in electrons (J/m³)"""
+
+    e_plasma_ions_thermal: float = 0.0
+    """Plasma thermal energy in ions (J)"""
+
+    eden_plasma_ions_thermal_vol_avg: float = 0.0
+    """Plasma volume averaged thermal energy density in ions (J/m³)"""
+
     betbm0: float = 1.5
     """leading coefficient for NB beta fraction"""
 
@@ -998,7 +1016,7 @@ class PhysicsData:
     """Plasma central ion temperature (Tᵢ₀) [keV]"""
 
     temp_plasma_ion_density_weighted_kev: float = 0.0
-    """Plasma density weighted average ion temperature (⟨Tᵢ⟩_n) [keV]"""
+    """Plasma density weighted average ion temperature (⟨Tᵢ⟩ₙ) [keV]"""
 
     f_temp_plasma_ion_electron: float = 1.0
     """Plasma ratio of ion temperature to electron temperature (used to calculate temp_plasma_ion_vol_avg_kev if `f_temp_plasma_ion_electron > 0.0`)"""

process/models/physics/confinement_time.py

@@ -139,7 +139,6 @@ def calculate_confinement_time(
         p_alpha_total_mw: float,
         aspect: float,
         b_plasma_toroidal_on_axis: float,
-        nd_plasma_ions_total_vol_avg: float,
         nd_plasma_electrons_vol_avg: float,
         nd_plasma_electron_line: float,
         eps: float,
@@ -155,11 +154,12 @@ def calculate_confinement_time(
         rmajor: float,
         rminor: float,
         temp_plasma_electron_density_weighted_kev: float,
-        temp_plasma_ion_density_weighted_kev: float,
         q95: float,
         qstar: float,
         vol_plasma: float,
         zeff: float,
+        eden_plasma_electrons_thermal_vol_avg: float,
+        eden_plasma_ions_thermal_vol_avg: float,
     ) -> tuple[float, float, float, float, float, float, float]:
         """Calculate the confinement times and the transport power loss terms.
 
@@ -173,8 +173,6 @@ def calculate_confinement_time(
             Aspect ratio
         b_plasma_toroidal_on_axis :
             Toroidal field on axis (T)
-        nd_plasma_ions_total_vol_avg :
-            Total ion density (/m3)
         nd_plasma_electrons_vol_avg :
             Volume averaged electron density (/m3)
         nd_plasma_electron_line :
@@ -209,19 +207,21 @@ def calculate_confinement_time(
             Plasma minor radius (m)
         temp_plasma_electron_density_weighted_kev :
             Density weighted average electron temperature (keV)
-        temp_plasma_ion_density_weighted_kev :
-            Density weighted average ion temperature (keV)
         vol_plasma :
             Plasma volume (m3)
         zeff :
             Plasma effective charge
+        eden_plasma_electrons_thermal_vol_avg :
+            Volume averaged electron thermal energy density [J/m³]
+        eden_plasma_ions_thermal_vol_avg :
+            Volume averaged ion thermal energy density [J/m³]
 
         Returns
         -------
         type
             Tuple containing:
-            - pden_electron_transport_loss_mw (float): Electron transport power (MW/m3)
-            - pden_ion_transport_loss_mw (float): Ion transport power (MW/m3)
+            - pden_electron_transport_loss_mw (float): Electron transport power (MW/m³)
+            - pden_ion_transport_loss_mw (float): Ion transport power (MW/m³)
             - t_electron_energy_confinement (float): Electron energy confinement time (s)
             - t_ion_energy_confinement (float): Ion energy confinement time (s)
             - t_energy_confinement (float): Global energy confinement time (s)
@@ -1060,24 +1060,13 @@ def calculate_confinement_time(
         # The transport losses is just the electron and ion thermal energies divided by
         # the confinement time.
         pden_ion_transport_loss_mw = (
-            (3 / 2)
-            * (constants.ELECTRON_CHARGE / 1e3)
-            * nd_plasma_ions_total_vol_avg
-            * temp_plasma_ion_density_weighted_kev
-            / t_ion_energy_confinement
-        )
+            eden_plasma_ions_thermal_vol_avg / t_ion_energy_confinement
+        ) / 1e6  # Convert from W/m³ to MW/m³
         pden_electron_transport_loss_mw = (
-            (3 / 2)
-            * (constants.ELECTRON_CHARGE / 1e3)
-            * nd_plasma_electrons_vol_avg
-            * temp_plasma_electron_density_weighted_kev
-            / t_electron_energy_confinement
-        )
+            eden_plasma_electrons_thermal_vol_avg / t_electron_energy_confinement
+        ) / 1e6  # Convert from W/m³ to MW/m³
 
-        ratio = (nd_plasma_ions_total_vol_avg / nd_plasma_electrons_vol_avg) * (
-            temp_plasma_ion_density_weighted_kev
-            / temp_plasma_electron_density_weighted_kev
-        )
+        ratio = eden_plasma_ions_thermal_vol_avg / eden_plasma_electrons_thermal_vol_avg
 
         # Global energy confinement time
 
@@ -1175,7 +1164,6 @@ def fhz(hfact: float) -> float:
                 p_alpha_total_mw=self.data.physics.p_alpha_total_mw,
                 aspect=self.data.physics.aspect,
                 b_plasma_toroidal_on_axis=self.data.physics.b_plasma_toroidal_on_axis,
-                nd_plasma_ions_total_vol_avg=self.data.physics.nd_plasma_ions_total_vol_avg,
                 nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg,
                 nd_plasma_electron_line=self.data.physics.nd_plasma_electron_line,
                 eps=self.data.physics.eps,
@@ -1191,11 +1179,12 @@ def fhz(hfact: float) -> float:
                 rmajor=self.data.physics.rmajor,
                 rminor=self.data.physics.rminor,
                 temp_plasma_electron_density_weighted_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,
-                temp_plasma_ion_density_weighted_kev=self.data.physics.temp_plasma_ion_density_weighted_kev,
                 q95=self.data.physics.q95,
                 qstar=self.data.physics.qstar,
                 vol_plasma=self.data.physics.vol_plasma,
                 zeff=self.data.physics.n_charge_plasma_effective_vol_avg,
+                eden_plasma_electrons_thermal_vol_avg=self.data.physics.eden_plasma_electrons_thermal_vol_avg,
+                eden_plasma_ions_thermal_vol_avg=self.data.physics.eden_plasma_ions_thermal_vol_avg,
             )
 
             # At power balance, fhz is zero.
@@ -1461,7 +1450,6 @@ def output_confinement_comparison(self, istell: int):
                 p_alpha_total_mw=self.data.physics.p_alpha_total_mw,
                 aspect=self.data.physics.aspect,
                 b_plasma_toroidal_on_axis=self.data.physics.b_plasma_toroidal_on_axis,
-                nd_plasma_ions_total_vol_avg=self.data.physics.nd_plasma_ions_total_vol_avg,
                 nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg,
                 nd_plasma_electron_line=self.data.physics.nd_plasma_electron_line,
                 eps=self.data.physics.eps,
@@ -1477,11 +1465,12 @@ def output_confinement_comparison(self, istell: int):
                 rmajor=self.data.physics.rmajor,
                 rminor=self.data.physics.rminor,
                 temp_plasma_electron_density_weighted_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,
-                temp_plasma_ion_density_weighted_kev=self.data.physics.temp_plasma_ion_density_weighted_kev,
                 q95=self.data.physics.q95,
                 qstar=self.data.physics.qstar,
                 vol_plasma=self.data.physics.vol_plasma,
                 zeff=self.data.physics.n_charge_plasma_effective_vol_avg,
+                eden_plasma_electrons_thermal_vol_avg=self.data.physics.eden_plasma_electrons_thermal_vol_avg,
+                eden_plasma_ions_thermal_vol_avg=self.data.physics.eden_plasma_ions_thermal_vol_avg,
             )
 
             try:

process/models/physics/physics.py

@@ -404,6 +404,36 @@ def run(self):
 
         self.beta.run()
 
+        # Stored thermal energies in the plasma
+
+        (
+            self.data.physics.eden_plasma_electrons_thermal_vol_avg,
+            self.data.physics.e_plasma_electrons_thermal,
+        ) = self.calaculate_stored_thermal_energy(
+            vol_plasma=self.data.physics.vol_plasma,
+            nd_plasma_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg,
+            temp_plasma_density_weighted_vol_avg_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,
+        )
+
+        (
+            self.data.physics.eden_plasma_ions_thermal_vol_avg,
+            self.data.physics.e_plasma_ions_thermal,
+        ) = self.calaculate_stored_thermal_energy(
+            vol_plasma=self.data.physics.vol_plasma,
+            nd_plasma_vol_avg=self.data.physics.nd_plasma_ions_total_vol_avg,
+            temp_plasma_density_weighted_vol_avg_kev=self.data.physics.temp_plasma_ion_density_weighted_kev,
+        )
+
+        self.data.physics.eden_plasma_thermal_vol_avg = (
+            self.data.physics.eden_plasma_electrons_thermal_vol_avg
+            + self.data.physics.eden_plasma_ions_thermal_vol_avg
+        )
+
+        self.data.physics.e_plasma_thermal_total = (
+            self.data.physics.e_plasma_electrons_thermal
+            + self.data.physics.e_plasma_ions_thermal
+        )
+
         # =======================================================
 
         # Set PF coil ramp times
@@ -825,7 +855,6 @@ def run(self):
             p_alpha_total_mw=self.data.physics.p_alpha_total_mw,
             aspect=self.data.physics.aspect,
             b_plasma_toroidal_on_axis=self.data.physics.b_plasma_toroidal_on_axis,
-            nd_plasma_ions_total_vol_avg=self.data.physics.nd_plasma_ions_total_vol_avg,
             nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg,
             nd_plasma_electron_line=self.data.physics.nd_plasma_electron_line,
             eps=self.data.physics.eps,
@@ -841,11 +870,12 @@ def run(self):
             rmajor=self.data.physics.rmajor,
             rminor=self.data.physics.rminor,
             temp_plasma_electron_density_weighted_kev=self.data.physics.temp_plasma_electron_density_weighted_kev,
-            temp_plasma_ion_density_weighted_kev=self.data.physics.temp_plasma_ion_density_weighted_kev,
             q95=self.data.physics.q95,
             qstar=self.data.physics.qstar,
             vol_plasma=self.data.physics.vol_plasma,
             zeff=self.data.physics.n_charge_plasma_effective_vol_avg,
+            eden_plasma_electrons_thermal_vol_avg=self.data.physics.eden_plasma_electrons_thermal_vol_avg,
+            eden_plasma_ions_thermal_vol_avg=self.data.physics.eden_plasma_ions_thermal_vol_avg,
         )
 
         # Total transport power from scaling law (MW)
@@ -2751,6 +2781,57 @@ def output_temperature_density_profile_info(self) -> None:
                 f"(pres_plasma_fuel_profile{i})",
                 self.data.physics.pres_plasma_fuel_profile[i],
             )
+
+        po.oblnkl(self.outfile)
+        po.ocmmnt(self.outfile, "----------------------------")
+
+        po.osubhd(self.outfile, "Stored thermal energies (Wₜₕ):")
+
+        po.ovarre(
+            self.outfile,
+            "Total plasma thermal energy [W]",
+            "(e_plasma_thermal_total)",
+            self.data.physics.e_plasma_thermal_total,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Plasma thermal energy in electrons [W]",
+            "(e_plasma_electrons_thermal)",
+            self.data.physics.e_plasma_electrons_thermal,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Plasma thermal energy in ions [W]",
+            "(e_plasma_ions_thermal)",
+            self.data.physics.e_plasma_ions_thermal,
+            "OP ",
+        )
+        po.oblnkl(self.outfile)
+
+        po.ovarre(
+            self.outfile,
+            "Volume averaged plasma thermal energy density [J/m³]",
+            "(eden_plasma_thermal_vol_avg)",
+            self.data.physics.eden_plasma_thermal_vol_avg,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Volume averaged plasma thermal energy density in electrons [J/m³]",
+            "(eden_plasma_electrons_thermal_vol_avg)",
+            self.data.physics.eden_plasma_electrons_thermal_vol_avg,
+            "OP ",
+        )
+        po.ovarre(
+            self.outfile,
+            "Volume averaged plasma thermal energy density in ions [J/m³]",
+            "(eden_plasma_ions_thermal_vol_avg)",
+            self.data.physics.eden_plasma_ions_thermal_vol_avg,
+            "OP ",
+        )
+
         po.oblnkl(self.outfile)
         po.ocmmnt(self.outfile, "----------------------------")
         po.oblnkl(self.outfile)
@@ -3011,6 +3092,42 @@ def calculate_total_plasma_heating_power(
             + p_hcd_injected_total_mw
         )
 
+    @staticmethod
+    def calaculate_stored_thermal_energy(
+        vol_plasma: float,
+        nd_plasma_vol_avg: float,
+        temp_plasma_density_weighted_vol_avg_kev: float,
+    ) -> tuple[float, float]:
+        """Calculate the stored thermal energy in the (W) plasma.
+
+        Parameters
+        ----------
+        vol_plasma : float
+            Plasma volume [m³].
+        nd_plasma_vol_avg : float
+            Volume averaged plasma density (⟨n⟩) [/m³].
+        temp_plasma_density_weighted_vol_avg_kev : float
+            Volume averaged density weighted plasma temperature (⟨T⟩ₙ) [keV].
+
+        Returns
+        -------
+        tuple[float, float]
+            A tuple containing:
+            - Volume averaged thermal energy density (J/m³).
+            - Total stored thermal energy in the plasma (J).
+
+        """
+        eden_plasma_thermal_vol_avg = (
+            1.5
+            * constants.KILOELECTRON_VOLT
+            * nd_plasma_vol_avg
+            * temp_plasma_density_weighted_vol_avg_kev
+        )
+
+        e_plasma_thermal = eden_plasma_thermal_vol_avg * vol_plasma
+
+        return eden_plasma_thermal_vol_avg, e_plasma_thermal
+
 
 def res_diff_time(rmajor, res_plasma, kappa95):
     """Calculates resistive diffusion time

process/models/stellarator/stellarator.py

@@ -2235,7 +2235,6 @@ def st_phys(self, output):
             self.data.physics.p_alpha_total_mw,
             self.data.physics.aspect,
             self.data.physics.b_plasma_toroidal_on_axis,
-            self.data.physics.nd_plasma_ions_total_vol_avg,
             self.data.physics.nd_plasma_electrons_vol_avg,
             self.data.physics.nd_plasma_electron_line,
             self.data.physics.eps,
@@ -2251,11 +2250,12 @@ def st_phys(self, output):
             self.data.physics.rmajor,
             self.data.physics.rminor,
             self.data.physics.temp_plasma_electron_density_weighted_kev,
-            self.data.physics.temp_plasma_ion_density_weighted_kev,
             self.data.stellarator.iotabar,
             self.data.physics.qstar,
             self.data.physics.vol_plasma,
             self.data.physics.n_charge_plasma_effective_vol_avg,
+            eden_plasma_electrons_thermal_vol_avg=self.data.physics.eden_plasma_electrons_thermal_vol_avg,
+            eden_plasma_ions_thermal_vol_avg=self.data.physics.eden_plasma_ions_thermal_vol_avg,
         )
 
         self.data.physics.ntau, self.data.physics.nTtau = (