Add H98, H20, H98L confinement quality factors.

docs/output.rst

@@ -241,6 +241,31 @@ analysis and inspection.
 ``W_thermal_tot`` (time) [J]:
   Total thermal stored energy.
 
+``tauE`` (time) [s]:
+  Thermal confinement time defined as ``W_thermal_tot`` / ``P_heating``, where
+  ``P_heating`` is the total heating power into the plasma, including external
+  contributions and fusion heating. Radiative losses are not subtracted from
+  heating power.
+
+``H98`` (time) [dimensionless]:
+  H-mode confinement quality factor with respect to the ITER98y2 scaling law,
+  defined as ``tauE`` / ``tau98_scaling``, where ``tau98_scaling`` is the
+  confinement time defined by the ITER98y2 scaling law, derived from the ITER
+  H-mode confinement database. As for ``tauE``, radiative losses are not
+  subtracted from the ``P_loss`` term used to calculate the empirical scaling
+  law confinement time.
+
+``H97L`` (time) [dimensionless]:
+  L-mode confinement quality factor with respect to the ITER97L scaling law
+  derived from the ITER L-mode confinement database. Defined similarly to ``H98``
+  above, but using the ITER97L scaling law for the confinement time.
+
+``H20`` (time) [dimensionless]:
+  H-mode confinement quality factor with respect to the ITER20 scaling law
+  derived from the updated (2020) ITER confinement database. Defined similarly
+  to ``H98`` above, but using the updated ITER20 scaling law law for the
+  confinement time.
+
 ``FFprime_face`` (time, rho_face) [m^2 T^2 / Wb]:
   :math:`FF'` on the face grid, where F is the toroidal flux function, and
   F' is its derivative with respect to the poloidal flux.

torax/physics.py

@@ -471,19 +471,7 @@ def calculate_plh_scaling_factor(
       density corresponding to the P_LH minimum.
   """
 
-  # Line-averaged electron density is poorly defined. In general, the definition
-  # is machine-dependent and even shot-dependent since it depends on the usage
-  # of a specific interferometry chord. Furthermore, even if we knew the
-  # specific chord used, its calculation would depend on magnetic geometry
-  # information beyond what is available in StandardGeometry.
-  # In lieu of a better solution, we use line-averaged electron density defined
-  # on the outer midplane.
-  Rmin_out = geo.Rout_face[-1] - geo.Rout_face[0]
-  line_avg_ne = (
-      core_profiles.nref
-      * _trapz(core_profiles.ne.face_value(), geo.Rout_face)
-      / Rmin_out
-  )
+  line_avg_ne = _calculate_line_avg_density(geo, core_profiles)
   # LH transition power for deuterium, in W. Eq 3 from Martin 2008.
   P_LH_hi_dens_D = (
       2.15
@@ -518,4 +506,130 @@ def calculate_plh_scaling_factor(
   return P_LH_hi_dens, P_LH_low_dens, ne_min_P_LH
 
 
+def calculate_scaling_law_confinement_time(
+    geo: Geometry,
+    core_profiles: state.CoreProfiles,
+    Ploss: jax.Array,
+    scaling_law: str,
+) -> jax.Array:
+  """Calculates the thermal energy confinement time for a given empirical scaling law.
+
+  Args:
+    geo: Torus geometry.
+    core_profiles: Core plasma profiles.
+    Ploss: Plasma power loss in MW.
+    scaling_law: Scaling law to use.
+
+  Returns:
+    Thermal energy confinement time in s.
+  """
+  scaling_params = {
+      'H98': {
+          # H98 empirical confinement scaling law:
+          # ITER Physics Expert Groups on Confinement and Transport and
+          # Confinement Modelling and Database, Nucl. Fusion 39 2175, 1999
+          # Doyle et al, Nucl. Fusion 47 (2007) S18–S127, Eq 30
+          'prefactor': 0.0562,
+          'Ip_exponent': 0.93,
+          'B_exponent': 0.15,
+          'line_avg_ne_exponent': 0.41,
+          'Ploss_exponent': -0.69,
+          'R_exponent': 1.97,
+          'inverse_aspect_ratio_exponent': 0.58,
+          'elongation_exponent': 0.78,
+          'effective_mass_exponent': 0.19,
+          'triangularity_exponent': 0.0,
+      },
+      'H97L': {
+          # From the ITER L-mode confinement database.
+          # S.M. Kaye et al 1997 Nucl. Fusion 37 1303, Eq 7
+          'prefactor': 0.023,
+          'Ip_exponent': 0.96,
+          'B_exponent': 0.03,
+          'line_avg_ne_exponent': 0.4,
+          'Ploss_exponent': -0.73,
+          'R_exponent': 1.83,
+          'inverse_aspect_ratio_exponent': -0.06,
+          'elongation_exponent': 0.64,
+          'effective_mass_exponent': 0.20,
+          'triangularity_exponent': 0.0,
+      },
+      'H20': {
+          # Updated ITER H-mode confinement database, using full dataset.
+          # G. Verdoolaege et al 2021 Nucl. Fusion 61 076006, Eq 7
+          'prefactor': 0.053,
+          'Ip_exponent': 0.98,
+          'B_exponent': 0.22,
+          'line_avg_ne_exponent': 0.24,
+          'Ploss_exponent': -0.669,
+          'R_exponent': 1.71,
+          'inverse_aspect_ratio_exponent': 0.35,
+          'elongation_exponent': 0.80,
+          'effective_mass_exponent': 0.20,
+          'triangularity_exponent': 0.36,  # (1+delta)^exponent
+      },
+  }
+
+  if scaling_law not in scaling_params:
+    raise ValueError(f'Unknown scaling law: {scaling_law}')
+
+  params = scaling_params[scaling_law]
+
+  Ip = core_profiles.currents.Ip_profile_face[-1] / 1e6  # in MA
+  B = geo.B0
+  line_avg_ne = _calculate_line_avg_density(geo, core_profiles) / 1e19
+  R = geo.Rmaj
+  inverse_aspect_ratio = geo.Rmin / geo.Rmaj
+
+  # Effective elongation definition. This is a different definition than
+  # the standard definition used in geo.elongation.
+  elongation = geo.area_face[-1] / (jnp.pi * geo.Rmin**2)
+  # TODO(b/317360834): extend when multiple ions are supported.
+  effective_mass = core_profiles.Ai
+  triangularity = geo.delta_face[-1]
+
+  tau_scaling = (
+      params['prefactor']
+      * Ip**params['Ip_exponent']
+      * B**params['B_exponent']
+      * line_avg_ne**params['line_avg_ne_exponent']
+      * Ploss**params['Ploss_exponent']
+      * R**params['R_exponent']
+      * inverse_aspect_ratio**params['inverse_aspect_ratio_exponent']
+      * elongation**params['elongation_exponent']
+      * effective_mass**params['effective_mass_exponent']
+      * (1 + triangularity)**params['triangularity_exponent']
+  )
+  return tau_scaling
+
+
+def _calculate_line_avg_density(
+    geo: Geometry,
+    core_profiles: state.CoreProfiles,
+) -> jax.Array:
+  """Calculates line-averaged electron density.
+
+  Line-averaged electron density is poorly defined. In general, the definition
+  is machine-dependent and even shot-dependent since it depends on the usage of
+  a specific interferometry chord. Furthermore, even if we knew the specific
+  chord used, its calculation would depend on magnetic geometry information
+  beyond what is available in StandardGeometry. In lieu of a better solution, we
+  use line-averaged electron density defined on the outer midplane.
+
+  Args:
+    geo: Torus geometry.
+    core_profiles: Core plasma profiles.
+
+  Returns:
+    Line-averaged electron density.
+  """
+  Rmin_out = geo.Rout_face[-1] - geo.Rout_face[0]
+  line_avg_ne = (
+      core_profiles.nref
+      * _trapz(core_profiles.ne.face_value(), geo.Rout_face)
+      / Rmin_out
+  )
+  return line_avg_ne
+
+
 # pylint: enable=invalid-name

torax/post_processing.py

@@ -366,6 +366,33 @@ def make_outputs(
       physics.calculate_plh_scaling_factor(geo, sim_state.core_profiles)
   )
 
+  # Thermal energy confinement time is the stored energy divided by the total
+  # input power into the plasma.
+
+  # Ploss term here does not include the reduction of radiated power. Most
+  # analysis of confinement times from databases have not included this term.
+  # Therefore highly radiative scenarios can lead to skewed results.
+
+  Ploss = (
+      integrated_sources['P_alpha_tot'] + integrated_sources['P_external_tot']
+  )
+  # TODO(b/380848256): include dW/dt term
+  tauE = W_thermal_tot / Ploss
+
+  tauH98 = physics.calculate_scaling_law_confinement_time(
+      geo, sim_state.core_profiles, Ploss/1e6, 'H98'
+  )
+  tauH97L = physics.calculate_scaling_law_confinement_time(
+      geo, sim_state.core_profiles, Ploss/1e6, 'H97L'
+  )
+  tauH20 = physics.calculate_scaling_law_confinement_time(
+      geo, sim_state.core_profiles, Ploss/1e6, 'H20'
+  )
+
+  H98 = tauE / tauH98
+  H97L = tauE / tauH97L
+  H20 = tauE / tauH20
+
   # Calculate total external (injected) and fusion (generated) energies based on
   # interval average.
   if previous_sim_state is not None:
@@ -408,6 +435,10 @@ def make_outputs(
       W_thermal_ion=W_thermal_ion,
       W_thermal_el=W_thermal_el,
       W_thermal_tot=W_thermal_tot,
+      tauE=tauE,
+      H98=H98,
+      H97L=H97L,
+      H20=H20,
       FFprime_face=FFprime_face,
       psi_norm_face=psi_norm_face,
       psi_face=sim_state.core_profiles.psi.face_value(),

torax/state.py

@@ -283,6 +283,13 @@ class PostProcessedOutputs:
     W_thermal_ion: Ion thermal stored energy [J]
     W_thermal_el: Electron thermal stored energy [J]
     W_thermal_tot: Total thermal stored energy [J]
+    tauE: Thermal energy confinement time [s]
+    H98: H-mode confinement quality factor with respect to the ITER98y2 scaling
+      law derived from the ITER H-mode confinement database
+    H97L: L-mode confinement quality factor with respect to the ITER97L scaling
+      law derived from the ITER H-mode confinement database
+    H20: H-mode confinement quality factor with respect to the ITER20 scaling
+      law derived from the updated (2020) ITER H-mode confinement database
     FFprime_face: FF' on the face grid, where F is the toroidal flux function
     psi_norm_face: Normalized poloidal flux on the face grid [Wb]
     psi_face: Poloidal flux on the face grid [Wb]
@@ -330,6 +337,10 @@ class PostProcessedOutputs:
   W_thermal_ion: array_typing.ScalarFloat
   W_thermal_el: array_typing.ScalarFloat
   W_thermal_tot: array_typing.ScalarFloat
+  tauE: array_typing.ScalarFloat
+  H98: array_typing.ScalarFloat
+  H97L: array_typing.ScalarFloat
+  H20: array_typing.ScalarFloat
   FFprime_face: array_typing.ArrayFloat
   psi_norm_face: array_typing.ArrayFloat
   # psi_face included in post_processed output for convenience, since the
@@ -377,6 +388,10 @@ def zeros(cls, geo: geometry.Geometry) -> PostProcessedOutputs:
         W_thermal_ion=jnp.array(0.0),
         W_thermal_el=jnp.array(0.0),
         W_thermal_tot=jnp.array(0.0),
+        tauE=jnp.array(0.0),
+        H98=jnp.array(0.0),
+        H97L=jnp.array(0.0),
+        H20=jnp.array(0.0),
         FFprime_face=jnp.zeros(geo.rho_face.shape),
         psi_norm_face=jnp.zeros(geo.rho_face.shape),
         psi_face=jnp.zeros(geo.rho_face.shape),