Add H89 L-mode confinement scaling to post-processing.

torax/physics.py

@@ -518,13 +518,26 @@ def calculate_scaling_law_confinement_time(
   Args:
     geo: Torus geometry.
     core_profiles: Core plasma profiles.
-    Ploss: Plasma power loss in MW.
+    Ploss: Plasma power loss in W.
     scaling_law: Scaling law to use.
 
   Returns:
     Thermal energy confinement time in s.
   """
   scaling_params = {
+      'H89P': {
+          # From Yushmanov et al, Nuclear Fusion, vol. 30, no. 10, pp. 4-6, 1990
+          'prefactor': 0.038128,
+          'Ip_exponent': 0.85,
+          'B_exponent': 0.2,
+          'line_avg_ne_exponent': 0.1,
+          'Ploss_exponent': -0.5,
+          'R_exponent': 1.5,
+          'inverse_aspect_ratio_exponent': 0.3,
+          'elongation_exponent': 0.5,
+          'effective_mass_exponent': 0.50,
+          'triangularity_exponent': 0.0,
+      },
       'H98': {
           # H98 empirical confinement scaling law:
           # ITER Physics Expert Groups on Confinement and Transport and
@@ -576,7 +589,8 @@ def calculate_scaling_law_confinement_time(
 
   params = scaling_params[scaling_law]
 
-  Ip = core_profiles.currents.Ip_profile_face[-1] / 1e6  # in MA
+  scaled_Ip = core_profiles.currents.Ip_profile_face[-1] / 1e6  # convert to MA
+  scaled_Ploss = Ploss / 1e6  # convert to MW
   B = geo.B0
   line_avg_ne = _calculate_line_avg_density(geo, core_profiles) / 1e19
   R = geo.Rmaj
@@ -591,10 +605,10 @@ def calculate_scaling_law_confinement_time(
 
   tau_scaling = (
       params['prefactor']
-      * Ip ** params['Ip_exponent']
+      * scaled_Ip ** params['Ip_exponent']
       * B ** params['B_exponent']
       * line_avg_ne ** params['line_avg_ne_exponent']
-      * Ploss ** params['Ploss_exponent']
+      * scaled_Ploss ** params['Ploss_exponent']
       * R ** params['R_exponent']
       * inverse_aspect_ratio ** params['inverse_aspect_ratio_exponent']
       * elongation ** params['elongation_exponent']

torax/post_processing.py

@@ -398,16 +398,20 @@ def make_outputs(
   # TODO(b/380848256): include dW/dt term
   tauE = W_thermal_tot / Ploss
 
+  tauH89P = physics.calculate_scaling_law_confinement_time(
+      geo, sim_state.core_profiles, Ploss, 'H89P'
+  )
   tauH98 = physics.calculate_scaling_law_confinement_time(
-      geo, sim_state.core_profiles, Ploss / 1e6, 'H98'
+      geo, sim_state.core_profiles, Ploss, 'H98'
   )
   tauH97L = physics.calculate_scaling_law_confinement_time(
-      geo, sim_state.core_profiles, Ploss / 1e6, 'H97L'
+      geo, sim_state.core_profiles, Ploss, 'H97L'
   )
   tauH20 = physics.calculate_scaling_law_confinement_time(
-      geo, sim_state.core_profiles, Ploss / 1e6, 'H20'
+      geo, sim_state.core_profiles, Ploss, 'H20'
   )
 
+  H89P = tauE / tauH89P
   H98 = tauE / tauH98
   H97L = tauE / tauH97L
   H20 = tauE / tauH20
@@ -487,6 +491,7 @@ def make_outputs(
       W_thermal_el=W_thermal_el,
       W_thermal_tot=W_thermal_tot,
       tauE=tauE,
+      H89P=H89P,
       H98=H98,
       H97L=H97L,
       H20=H20,

torax/state.py

@@ -276,10 +276,12 @@ class PostProcessedOutputs:
     W_thermal_el: Electron thermal stored energy [J]
     W_thermal_tot: Total thermal stored energy [J]
     tauE: Thermal energy confinement time [s]
+    H89P: L-mode confinement quality factor with respect to the ITER89P scaling
+      law derived from the ITER L-mode confinement database
     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
+      law derived from the ITER L-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
@@ -337,6 +339,7 @@ class PostProcessedOutputs:
   W_thermal_el: array_typing.ScalarFloat
   W_thermal_tot: array_typing.ScalarFloat
   tauE: array_typing.ScalarFloat
+  H89P: array_typing.ScalarFloat
   H98: array_typing.ScalarFloat
   H97L: array_typing.ScalarFloat
   H20: array_typing.ScalarFloat
@@ -395,6 +398,7 @@ def zeros(cls, geo: geometry.Geometry) -> PostProcessedOutputs:
         W_thermal_el=jnp.array(0.0),
         W_thermal_tot=jnp.array(0.0),
         tauE=jnp.array(0.0),
+        H89P=jnp.array(0.0),
         H98=jnp.array(0.0),
         H97L=jnp.array(0.0),
         H20=jnp.array(0.0),

torax/tests/physics.py

@@ -488,8 +488,11 @@ def test_calculate_scaling_law_confinement_time(self, elongation_LCFS):
             Ip_profile_face=jnp.ones_like(geo.rho_face_norm) * 10e6,
         ),
     )
-    Ploss = jnp.array(50.0)
+    Ploss = jnp.array(50e6)
 
+    H89P = physics.calculate_scaling_law_confinement_time(
+        geo, core_profiles, Ploss, 'H89P'
+    )
     H98 = physics.calculate_scaling_law_confinement_time(
         geo, core_profiles, Ploss, 'H98'
     )
@@ -499,6 +502,19 @@ def test_calculate_scaling_law_confinement_time(self, elongation_LCFS):
     H20 = physics.calculate_scaling_law_confinement_time(
         geo, core_profiles, Ploss, 'H20'
     )
+
+    expected_H89P = (
+        0.038128
+        * 10**0.85
+        * 5**0.2
+        * 20**0.1
+        * 50**-0.5
+        * 6**1.5
+        * (1 / 3) ** 0.3
+        * 3**0.50
+        * elongation_LCFS**0.50
+    )
+
     expected_H98 = (
         0.0562
         * 10**0.93
@@ -535,6 +551,7 @@ def test_calculate_scaling_law_confinement_time(self, elongation_LCFS):
         * elongation_LCFS**0.80
     )
     # pylint: enable=invalid-name
+    np.testing.assert_allclose(H89P, expected_H89P)
     np.testing.assert_allclose(H98, expected_H98)
     np.testing.assert_allclose(H97L, expected_H97L)
     np.testing.assert_allclose(H20, expected_H20)