Add plotting and documentation for cyclotron radiation sink.

Includes new post-processed variable, leading to regeneration of all sim tests.

PiperOrigin-RevId: 712888033

docs/configuration.rst

@@ -931,7 +931,7 @@ Time dependent Gaussian pellet source. No first-principle-based model is yet imp
   Total particle source in units of particles/s
 
 generic_particle_source
-^^^^^^^^^^^^^^^^^^^
+^^^^^^^^^^^^^^^^^^^^^^^
 
 Time dependent Gaussian particle source. No first-principle-based model is yet implemented in TORAX.
 
@@ -988,6 +988,26 @@ Bremsstrahlung model from Wesson, with an optional correction for relativistic e
 
 ``use_relativistic_correction`` (bool = False)
 
+cyclotron_radiation_heat_sink
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Cyclotron radiation model from Albajar NF 2001 with a deposition profile from Artaud NF 2018.
+
+``mode`` (str = 'model')
+
+``wall_reflection_coeff`` (float = 0.9)
+  Machine-dependent dimensionless parameter corresponding to the fraction of
+  cyclotron radiation reflected off the wall and reabsorbed by the plasma.
+
+``beta_min`` (float = 0.5)
+
+``beta_max`` (float = 8.0)
+
+``beta_grid_size`` (int = 32)
+  beta in this context is a variable in the temperature profile parameterization used
+  in the Albajar model. The parameter is fit with simple grid search performed over
+  the range ``[beta_min, beta_max]``, with ``beta_grid_size`` uniformly spaced steps.
+
 electron_cyclotron_source
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 Electron-cyclotron heating and current drive, based on the local efficiency model in `Lin-Liu et al., 2003 <https://doi.org/10.1063/1.1610472>`_.

docs/physics_models.rst

@@ -409,17 +409,25 @@ Presently, TORAX provides three built-in formula-based particle sources for the
 
 Radiation
 ---------
-Currently, TORAX only has a dedicated model for Bremsstrahlung. Models for cyclotron radiation,
-recombination, and line radiation are still left for future work.
+Currently, TORAX has dedicated models for Bremsstrahlung and cyclotron radiation.
+Models for line radiation are left for future work.
 
 Bremsstrahlung
 ^^^^^^^^^^^^^^
 
 Uses the model from Wesson, John, and David J. Campbell. Tokamaks. Vol. 149.
 An optional correction for relativistic effects from Stott PPCF 2005 can be enabled with the flag "use_relativistic_correction".
 
-Ion Cyclotron Resonance
-^^^^^^^^^^^^^^^^^^^^^^^
+Cyclotron Radiation
+^^^^^^^^^^^^^^^^^^^
+
+Uses the total radiation power from `Albajar NF 2001 <https://doi.org/10.1088/0029-5515/41/6/301>`_
+with a deposition profile from `Artaud NF 2018 <https://doi.org/10.1088/1741-4326/aad5b1>`_.
+The Albajar model includes a parameterization of the temperature profile which in TORAX is fit by simple
+grid search for computational efficiency.
+
+Ion Cyclotron Resonance Heating
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Presently this source is implemented for a SPARC specific ICRH scenario.
 

torax/plotting/configs/default_plot_config.py

@@ -145,19 +145,32 @@
             suppress_zero_values=True,  # Do not plot all-zero data
         ),
         plotruns_lib.PlotProperties(
-            attrs=('q_brems', 'q_rad'),
-            labels=(r'$Q_\mathrm{brems}$', r'$Q_\mathrm{rad}$'),
+            attrs=('q_brems', 'q_rad', 'q_cycl'),
+            labels=(
+                r'$Q_\mathrm{brems}$',
+                r'$Q_\mathrm{rad}$',
+                r'$Q_\mathrm{cycl}$',
+            ),
             ylabel=r'Heat sink density $[MW~m^{-3}]$',
             suppress_zero_values=True,  # Do not plot all-zero data
         ),
         plotruns_lib.PlotProperties(
             plot_type=plotruns_lib.PlotType.TIME_SERIES,
-            attrs=('p_auxiliary', 'p_ohmic', 'p_alpha', 'p_sink'),
+            attrs=(
+                'p_auxiliary',
+                'p_ohmic',
+                'p_alpha',
+                'p_brems',
+                'p_rad',
+                'p_cycl',
+            ),
             labels=(
                 r'$P_\mathrm{aux}$',
                 r'$P_\mathrm{ohm}$',
                 r'$P_\mathrm{\alpha}$',
-                r'$P_\mathrm{sink}$',
+                r'$P_\mathrm{brems}$',
+                r'$P_\mathrm{rad}$',
+                r'$P_\mathrm{cycl}$',
             ),
             ylabel=r'Total heating/sink powers $[MW]$',
             legend_fontsize=6,  # Smaller fontsize for this plot

torax/plotting/configs/sources_plot_config.py

@@ -58,12 +58,21 @@
         ),
         plotruns_lib.PlotProperties(
             plot_type=plotruns_lib.PlotType.TIME_SERIES,
-            attrs=('p_auxiliary', 'p_ohmic', 'p_alpha', 'p_sink'),
+            attrs=(
+                'p_auxiliary',
+                'p_ohmic',
+                'p_alpha',
+                'p_rad',
+                'p_brems',
+                'p_cycl',
+            ),
             labels=(
                 r'$P_\mathrm{aux}$',
                 r'$P_\mathrm{ohm}$',
                 r'$P_\mathrm{\alpha}$',
-                r'$P_\mathrm{sink}$',
+                r'$P_\mathrm{rad}$',
+                r'$P_\mathrm{brems}$',
+                r'$P_\mathrm{cycl}$',
             ),
             ylabel=r'Total heating/sink powers $[MW]$',
             legend_fontsize=8,  # Smaller fontsize for this plot
@@ -99,8 +108,12 @@
             suppress_zero_values=True,  # Do not plot all-zero data
         ),
         plotruns_lib.PlotProperties(
-            attrs=('q_brems', 'q_rad'),
-            labels=(r'$Q_\mathrm{brems}$', r'$Q_\mathrm{rad}$'),
+            attrs=('q_brems', 'q_rad', 'q_cycl'),
+            labels=(
+                r'$Q_\mathrm{brems}$',
+                r'$Q_\mathrm{rad}$',
+                r'$Q_\mathrm{cycl}$',
+            ),
             ylabel=r'Heat sink density $[MW~m^{-3}]$',
             suppress_zero_values=True,  # Do not plot all-zero data
         ),

torax/plotting/plotruns_lib.py

@@ -112,6 +112,7 @@ class PlotData:
   q_alpha_e: np.ndarray  # [MW/m^3]
   q_ohmic: np.ndarray  # [MW/m^3]
   q_brems: np.ndarray  # [MW/m^3]
+  q_cycl: np.ndarray  # [MW/m^3]
   q_ei: np.ndarray  # [MW/m^3]
   q_rad: np.ndarray  # [MW/m^3]
   Q_fusion: np.ndarray  # pylint: disable=invalid-name  # Dimensionless
@@ -126,6 +127,9 @@ class PlotData:
   p_ohmic: np.ndarray  # [MW]
   p_alpha: np.ndarray  # [MW]
   p_sink: np.ndarray  # [MW]
+  p_brems: np.ndarray  # [MW]
+  p_cycl: np.ndarray  # [MW]
+  p_rad: np.ndarray  # [MW]
   t: np.ndarray  # [s]
   rho_cell_coord: np.ndarray  # Normalized toroidal flux coordinate
   rho_face_coord: np.ndarray  # Normalized toroidal flux coordinate
@@ -180,12 +184,14 @@ def _transform_data(ds: xr.Dataset):
         'fusion_heat_source_el': 1e6,  # W/m^3 to MW/m^3
         'ohmic_heat_source': 1e6,  # W/m^3 to MW/m^3
         'bremsstrahlung_heat_sink': 1e6,  # W/m^3 to MW/m^3
+        'cyclotron_radiation_heat_sink': 1e6,  # W/m^3 to MW/m^3
         'impurity_radiation_heat_sink': 1e6,  # W/m^3 to MW/m^3
         'qei_source': 1e6,  # W/m^3 to MW/m^3
         'P_ohmic': 1e6,  # W to MW
         'P_external_tot': 1e6,  # W to MW
         'P_alpha_tot': 1e6,  # W to MW
         'P_brems': 1e6,  # W to MW
+        'P_cycl': 1e6,  # W to MW
         'P_ecrh': 1e6,  # W to MW
         'P_rad': 1e6,  # W to MW
         'I_ecrh': 1e6,  # A to MA
@@ -258,6 +264,9 @@ def _transform_data(ds: xr.Dataset):
       q_brems=get_optional_data(
           core_sources_dataset, 'bremsstrahlung_heat_sink', 'cell'
       ),
+      q_cycl=get_optional_data(
+          core_sources_dataset, 'cyclotron_radiation_heat_sink', 'cell'
+      ),
       q_rad=get_optional_data(
           core_sources_dataset, 'impurity_radiation_heat_sink', 'cell'
       ),
@@ -279,7 +288,11 @@ def _transform_data(ds: xr.Dataset):
       ).to_numpy(),
       p_alpha=post_processed_outputs_dataset['P_alpha_tot'].to_numpy(),
       p_sink=post_processed_outputs_dataset['P_brems'].to_numpy()
-      + post_processed_outputs_dataset['P_rad'].to_numpy(),
+      + post_processed_outputs_dataset['P_rad'].to_numpy()
+      + post_processed_outputs_dataset['P_cycl'].to_numpy(),
+      p_brems=post_processed_outputs_dataset['P_brems'].to_numpy(),
+      p_rad=post_processed_outputs_dataset['P_rad'].to_numpy(),
+      p_cycl=post_processed_outputs_dataset['P_cycl'].to_numpy(),
       te_volume_avg=post_processed_outputs_dataset['te_volume_avg'].to_numpy(),
       ti_volume_avg=post_processed_outputs_dataset['ti_volume_avg'].to_numpy(),
       ne_volume_avg=post_processed_outputs_dataset['ne_volume_avg'].to_numpy(),

torax/post_processing.py

@@ -39,6 +39,7 @@
 EL_HEAT_SOURCE_TRANSFORMATIONS = {
     'ohmic_heat_source': 'P_ohmic',
     'bremsstrahlung_heat_sink': 'P_brems',
+    'cyclotron_radiation_heat_sink': 'P_cycl',
     'electron_cyclotron_source': 'P_ecrh',
     'impurity_radiation_heat_sink': 'P_rad',
 }

torax/sources/cyclotron_radiation_heat_sink.py

@@ -45,11 +45,11 @@ class RuntimeParams(runtime_params_lib.RuntimeParams):
   wall_reflection_coeff: float = 0.9
 
   # The beta parameter is used in the parameterized function for the
-  # temperature fit. beta_min, beta_max, and beta_grid_length are used for a
+  # temperature fit. beta_min, beta_max, and beta_grid_size are used for a
   # grid search to find the best fit.
   beta_min: float = 0.5
   beta_max: float = 8.0
-  beta_grid_length: int = 32
+  beta_grid_size: int = 32
   mode: runtime_params_lib.Mode = runtime_params_lib.Mode.MODEL_BASED
 
   def make_provider(
@@ -64,7 +64,7 @@ def build_static_params(self) -> 'StaticRuntimeParams':
         is_explicit=self.is_explicit,
         beta_min=self.beta_min,
         beta_max=self.beta_max,
-        beta_grid_length=self.beta_grid_length,
+        beta_grid_size=self.beta_grid_size,
     )
 
 
@@ -87,7 +87,7 @@ def build_dynamic_params(
 class StaticRuntimeParams(runtime_params_lib.StaticRuntimeParams):
   beta_min: float
   beta_max: float
-  beta_grid_length: int
+  beta_grid_size: int
 
 
 @chex.dataclass(frozen=True)
@@ -249,7 +249,7 @@ def _solve_alpha_t_beta_t_grid_search(
   Args:
     rho_norm: Normalized toroidal flux coordinate.
     te_data: The temperature data to be fit, assumed to be on the face grid.
-    beta_scan_parameters: A tuple of (beta_min, beta_max, beta_grid_length)
+    beta_scan_parameters: A tuple of (beta_min, beta_max, beta_grid_size)
       parameters for the grid search.
 
   Returns:
@@ -352,7 +352,7 @@ def cyclotron_radiation_albajar(
   beta_scan_parameters = (
       static_source_runtime_params.beta_min,
       static_source_runtime_params.beta_max,
-      static_source_runtime_params.beta_grid_length,
+      static_source_runtime_params.beta_grid_size,
   )
   alpha_t, beta_t = _solve_alpha_t_beta_t_grid_search(
       rho_norm=static_runtime_params_slice.torax_mesh.face_centers,

torax/state.py

@@ -315,6 +315,7 @@ class PostProcessedOutputs:
     P_alpha_tot: Total fusion power to plasma [W]
     P_ohmic: Ohmic heating power to electrons [W]
     P_brems: Bremsstrahlung electron heat sink [W]
+    P_cycl: Cyclotron radiation electron heat sink [W]
     P_ecrh: Total electron cyclotron source power [W]
     P_rad: Impurity radiation heat sink [W]
     I_ecrh: Total electron cyclotron source current [A]
@@ -370,6 +371,7 @@ class PostProcessedOutputs:
   P_alpha_tot: array_typing.ScalarFloat
   P_ohmic: array_typing.ScalarFloat
   P_brems: array_typing.ScalarFloat
+  P_cycl: array_typing.ScalarFloat
   P_ecrh: array_typing.ScalarFloat
   P_rad: array_typing.ScalarFloat
   I_ecrh: array_typing.ScalarFloat
@@ -424,6 +426,7 @@ def zeros(cls, geo: geometry.Geometry) -> PostProcessedOutputs:
         P_alpha_tot=jnp.array(0.0),
         P_ohmic=jnp.array(0.0),
         P_brems=jnp.array(0.0),
+        P_cycl=jnp.array(0.0),
         P_ecrh=jnp.array(0.0),
         P_rad=jnp.array(0.0),
         I_ecrh=jnp.array(0.0),