Set electron convectivity proportional to diffusivity

docs/physics_models.rst

@@ -138,8 +138,7 @@ TORAX currently offers three transport models:
     The particle convection velocity, :math:`V_e`, is user-defined.
 
   - **Bohm-GyroBohm:** A widely used semi-empirical model summing terms proportional to Bohm
-    and gyro-Bohm scaling factors. We use the implementation from `Tholerus et al <https://doi.org/10.1088/1741-4326/ad6ea2>`_,
-    Section 3.3.
+    and gyro-Bohm scaling factors (`Erba et al, 1998 <doi.org/10.1088/0029-5515/38/7/305>`_).
 
     The heat diffusivities for electrons and ions are given by:
 
@@ -182,7 +181,7 @@ TORAX currently offers three transport models:
     coordinate, :math:`p_e` is the electron pressure, and :math:`T_e` is the
     electron temperature.
 
-    The electron diffusivity is given by:
+    The electron diffusivity is given by `Garzotti et al, 2003 <doi.org/10.1088/0029-5515/43/12/025>`_:
 
     .. math::
       D_e = \eta \frac{\chi_e \chi_i}{\chi_e + \chi_i}
@@ -193,15 +192,15 @@ TORAX currently offers three transport models:
 
       \eta = c_1 + (c_2 - c_1) \rho_{\text{tor}}
 
-    where :math:`c_1` and :math:`c_2` are constants.
+    where :math:`c_1` and :math:`c_2` are user-defined parameters.
 
-    The electron convectivity is given by:
+    There is little discussion in the literature about setting the electron convectivity from the Bohm/gyro-Bohm model.
+    Following RAPTOR's `vpdn_chiescal` method, in TORAX, we set the electron convectivity proportional to the diffusivity,
 
     .. math::
-      v_e = \frac{1}{2} \frac{D_e A^2}{V \frac{dV}{d\rho}}
+      v_e = c_v D_e
 
-    where :math:`A` and :math:`V` are the area and volume of the flux surface
-    respectively.
+    where :math:`c_v` is a user-defined parameter.
 
     The default values for the model parameters are as follows:
 
@@ -213,7 +212,9 @@ TORAX currently offers three transport models:
 
     :math:`c_2 = 0.3`
 
-    We note that the Bohm-GyroBohm model TORAX implementation is presently
+    :math:`c_v = 0.1`
+
+    Please note that the Bohm-GyroBohm model TORAX implementation is presently
     experimental and subject to ongoing verification against established simulations.
 
   - **QLKNN:** This is a ML-surrogate model trained on a large dataset of the `QuaLiKiz <https://gitlab.com/qualikiz-group/QuaLiKiz>`_

torax/transport_model/bohm_gyrobohm.py

@@ -52,6 +52,9 @@ class RuntimeParams(runtime_params_lib.RuntimeParams):
   # Constants for the electron diffusivity weighting factor.
   d_face_c1: runtime_params_lib.TimeInterpolatedInput = 1.0
   d_face_c2: runtime_params_lib.TimeInterpolatedInput = 0.3
+  # Prefactor for the electron convectivity
+  v_face_coeff: runtime_params_lib.TimeInterpolatedInput = 0.1
+
 
   def make_provider(
       self, torax_mesh: geometry.Grid1D | None = None
@@ -70,6 +73,7 @@ class RuntimeParamsProvider(runtime_params_lib.RuntimeParamsProvider):
   chi_i_gyrobohm_coeff: interpolated_param.InterpolatedVarSingleAxis
   d_face_c1: interpolated_param.InterpolatedVarSingleAxis
   d_face_c2: interpolated_param.InterpolatedVarSingleAxis
+  v_face_coeff: interpolated_param.InterpolatedVarSingleAxis
 
   def build_dynamic_params(self, t: chex.Numeric) -> DynamicRuntimeParams:
     return DynamicRuntimeParams(**self.get_dynamic_params_kwargs(t))
@@ -85,6 +89,7 @@ class DynamicRuntimeParams(runtime_params_lib.DynamicRuntimeParams):
   chi_i_gyrobohm_coeff: array_typing.ScalarFloat
   d_face_c1: array_typing.ScalarFloat
   d_face_c2: array_typing.ScalarFloat
+  v_face_coeff: array_typing.ScalarFloat
 
   def sanity_check(self):
     runtime_params_lib.DynamicRuntimeParams.sanity_check(self)
@@ -210,19 +215,8 @@ def _call_implementation(
         / (chi_e[1:] + chi_i[1:] + constants_module.CONSTANTS.eps),
     ])
 
-    # Convection
-    # v_face_el is also zero on-axis by definition
-    # To avoid 0/0, we set the first element to 0 manually
-    # This definition for pinch velocity is from Tholerus et al.
-    # Pinch velocity = inward = -ve in TORAX
-    v_face_el = -jnp.concatenate([
-        jnp.zeros(1),
-        0.5
-        * d_face_el[1:]
-        * geo.area_face[1:] ** 2
-        * geo.rho_b
-        / (geo.volume_face[1:] * geo.vpr_face[1:]),
-    ])
+    # Electron convectivity set proportional to the electron diffusivity
+    v_face_el = dynamic_runtime_params_slice.transport.v_face_coeff * d_face_el
 
     return state.CoreTransport(
         chi_face_ion=chi_i,