add conv_vel cs limit to mlt/tdc

docs/source/changelog.rst

@@ -42,6 +42,14 @@ Previous MESA versions exclusively adopted the latter.
 By user request, a radiation diffusion flux limiter has been reintroduced with the 
 control ``use_flux_limiting_with_dPrad_dm_form`` for use with ``use_dPrad_dm_form_of_T_gradient_eqn``.
 
+By user request, an option for limiting the convective velocity predicted by mixing length theories has been introduced 
+allowing users to limit the convective velocity to some fraction of the local sound speed using the
+controls `max_conv_vel_div_csound` and `max_conv_vel_div_csound_maxq`.
+
+By user request, and motivated by the underestimation of line opacities from expanding material by the `Ferguson et al. (2005) <https://ui.adsabs.harvard.edu/abs/1994ApJ...437..879A/abstract>`_ tables, 
+see also section 2.2 in `Morozova et al. (2015) <https://ui.adsabs.harvard.edu/abs/2015ApJ...814...63M/abstract>`_. An optional control for an opacity floor, ``opacity_min``, has
+been introduced.
+
 .. _Bug Fixes main:
 
 Bug Fixes
@@ -290,7 +298,7 @@ shmesa
 ~~~~~~
 
 We have introduced a new set of command line utilities for interacting with MESA.
-See the README in ``$MESA_DIR/scripts/shmesa``, or online `here <https://github.com/MESAHub/mesa/tree/main/scripts/shmesa>`__.
+See the README in ``$MESA_DIR/scripts/shmesa``, or online `here <https://github.com/MESAHub/mesa/tree/main/scripts/shmesa>`_.
 
 These utilities provide functionality such as changing inlist parameters (``shmesa change``) or filling in the full
 ``run_star_extras.f90`` template (``shmesa extras``).

star/defaults/controls.defaults

@@ -781,16 +781,6 @@
     burn_min2 = 1000
 
 
-      ! max_conv_vel_div_csound_maxq
-      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      ! only consider from center out to this location
-
-      ! ::
-
-    max_conv_vel_div_csound_maxq = 1
-
-
       ! width_for_limit_conv_vel
       ! ~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -3378,13 +3368,22 @@
       ! max_conv_vel_div_csound
       ! ~~~~~~~~~~~~~~~~~~~~~~~
 
-      ! convective velocities are limited to local sound speed times this factor
+      ! convective velocities are limited to local sound speed times this factor.
 
       ! ::
 
     max_conv_vel_div_csound = 1d99
 
 
+      ! max_conv_vel_div_csound_maxq
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! only consider ``max_conv_vel_div_csound`` from center out to this location
+
+      ! ::
+
+    max_conv_vel_div_csound_maxq = 1d99
+
       ! max_v_for_convection
       ! ~~~~~~~~~~~~~~~~~~~~
 
@@ -7408,10 +7407,19 @@
     use_starting_composition_for_kap = .false.
 
 
+      ! opacity_min
+      ! ~~~~~~~~~~~
+
+      ! limit minimum opacities to this value (ignore this if value is < 0)
+
+      ! ::
+
+    opacity_min = -1
+
       ! opacity_max
       ! ~~~~~~~~~~~
 
-      ! limit opacities to this value (ignore this is value is < 0)
+      ! limit opacities to this value (ignore this if value is < 0)
 
       ! ::
 

star/job/run_star_support.f90

@@ -2833,6 +2833,9 @@ subroutine do_star_job_controls_after(id, s, restart, pgstar_ok, ierr)
          if (s% opacity_max > 0) &
             write(*,1) 'opacity_max', s% opacity_max
 
+         if (s% opacity_min > 0) &
+            write(*,1) 'opacity_min', s% opacity_min
+
          if (s% job% show_net_reactions_info) then
             write(*,'(a)') ' net reactions '
             call show_net_reactions_and_info(s% net_handle, 6, ierr)

star/private/ctrls_io.f90

@@ -336,7 +336,7 @@ module ctrls_io
 
     ! hydro parameters
     energy_eqn_option, &
-    opacity_factor, opacity_max, min_logT_for_opacity_factor_off, min_logT_for_opacity_factor_on, &
+    opacity_factor, opacity_min, opacity_max, min_logT_for_opacity_factor_off, min_logT_for_opacity_factor_on, &
     max_logT_for_opacity_factor_on, max_logT_for_opacity_factor_off, &
     non_nuc_neu_factor, &
     use_time_centered_eps_grav, &
@@ -1821,6 +1821,7 @@ subroutine store_controls(s, ierr)
  ! hydro parameters
  s% energy_eqn_option = energy_eqn_option
  s% opacity_factor = opacity_factor
+ s% opacity_min = opacity_min
  s% opacity_max = opacity_max
  s% min_logT_for_opacity_factor_off = min_logT_for_opacity_factor_off
  s% min_logT_for_opacity_factor_on = min_logT_for_opacity_factor_on
@@ -3503,6 +3504,7 @@ subroutine set_controls_for_writing(s, ierr)
 
  ! hydro parameters
  energy_eqn_option = s% energy_eqn_option
+ opacity_min = s% opacity_min
  opacity_max = s% opacity_max
  opacity_factor = s% opacity_factor
  min_logT_for_opacity_factor_off = s% min_logT_for_opacity_factor_off

star/private/micro.f90

@@ -653,6 +653,12 @@ subroutine do_kap_for_cell(s,k,ierr)
        s% d_opacity_dlnT(k) = 0
     end if
 
+    if (s% opacity(k) < s% opacity_min .and. s% opacity_min > 0) then
+       s% opacity(k) = s% opacity_min
+       s% d_opacity_dlnd(k) = 0
+       s% d_opacity_dlnT(k) = 0
+    end if
+
     if (is_bad_num(s% opacity(k))) then
        if (s% stop_for_bad_nums) then
           !$omp critical (star_kap_get_bad_num)

star/private/pre_ms_model.f90

@@ -242,11 +242,13 @@ subroutine build_pre_ms_model(id, s, nvar_hydro, species, ierr)
 
          mstar1 = rpar(1)
 
+! these pointers are set, but none of these vars are set.
          xh => s% xh
          q => s% q
          dq => s% dq
          nz = s% nz
 
+! this debug statement will cause a backtrace because the above vars are not set.
          if (dbg) then
             write(*,'(A)')
             write(*,*) 'finished pre-MS model'

star/private/turb_support.f90

@@ -184,7 +184,7 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
       integer, intent(out) :: ierr
 
       type(auto_diff_real_star_order1) :: Pr, Pg, grav, Lambda, gradL, beta
-      real(dp) :: conv_vel_start, scale
+      real(dp) :: conv_vel_start, scale, max_conv_vel
 
       ! these are used by use_superad_reduction
       real(dp) :: Gamma_limit, scale_value1, scale_value2, diff_grads_limit, reduction_limit, lambda_limit
@@ -201,12 +201,23 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
       beta = Pg / P
       Lambda = mixing_length_alpha*scale_height
       grav = cgrav*m/pow2(r)
+      max_conv_vel = 1d99
       if (s% use_Ledoux_criterion) then
          gradL = grada + gradL_composition_term  ! Ledoux temperature gradient
       else
          gradL = grada
       end if
 
+      ! maximum convection velocity.
+      if (k>=1) then
+         if (s% q(k) <= s% max_conv_vel_div_csound_maxq) then
+             max_conv_vel = s% csound_face(k) * s% max_conv_vel_div_csound
+         else
+            max_conv_vel = 1d99
+         end if
+      end if
+
+
       ! Initialize with no mixing
       mixing_type = no_mixing
       gradT = gradr
@@ -264,7 +275,7 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
          call set_TDC(&
             conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, &
             mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
-            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), ierr)
+            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, ierr)
          s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
 
             if (ierr /= 0) then
@@ -282,7 +293,7 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
                call set_TDC(&
                   conv_vel_start, mixing_length_alpha, s% alpha_TDC_DAMP, s%alpha_TDC_DAMPR, s%alpha_TDC_PtdVdt, s%dt, cgrav, m, report, &
                   mixing_type, scale, chiT, chiRho, gradr_scaled, r, P, T, rho, dV, Cp, opacity, &
-                  scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), ierr)
+                  scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, ierr)
                s% dvc_dt_TDC(k) = (conv_vel%val - conv_vel_start) / s%dt
                if (ierr /= 0) then
                   if (s% report_ierr) write(*,*) 'ierr from set_TDC when using superad_reduction'
@@ -296,7 +307,8 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
          call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, &
                         chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
                         gradr, grada, gradL, &
-                        Gamma, gradT, Y_face, conv_vel, D, mixing_type, ierr)
+                        Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
+
 
          if (ierr /= 0) then
             if (s% report_ierr) write(*,*) 'ierr from set_MLT'
@@ -313,7 +325,8 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
                call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, &
                               chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
                               gradr_scaled, grada, gradL, &
-                              Gamma, gradT, Y_face, conv_vel, D, mixing_type, ierr)
+                              Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
+
                if (ierr /= 0) then
                   if (s% report_ierr) write(*,*) 'ierr from set_MLT when using superad_reduction'
                   return

star_data/private/star_controls.inc

@@ -796,7 +796,7 @@
          logical :: use_simple_es_for_kap, use_starting_composition_for_kap
 
          real(dp) :: min_kap_for_dPrad_dm_eqn
-         real(dp) :: opacity_max, opacity_factor
+         real(dp) :: opacity_min, opacity_max, opacity_factor
          real(dp) :: min_logT_for_opacity_factor_off
          real(dp) :: min_logT_for_opacity_factor_on
          real(dp) :: max_logT_for_opacity_factor_on

turb/private/mlt.f90

@@ -57,25 +57,25 @@ module MLT
    !! @param D The chemical diffusion coefficient (cm^2/s).
    !! @param mixing_type Set to convective if convection operates (output).
    !! @param ierr Tracks errors (output).
-   subroutine calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
-                     chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
-                     gradr, grada, gradL, &
-                     Gamma, gradT, Y_face, conv_vel, D, mixing_type, ierr)
+   subroutine calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, &
+                     Henyey_MLT_y_param, chiT, chiRho, Cp, grav, Lambda, rho, &
+                     P, T, opacity, gradr, grada, gradL, Gamma, gradT, Y_face, &
+                     conv_vel, D, mixing_type, max_conv_vel, ierr)
       use const_def, only: dp, clight, convective_mixing, crad, two_13, four_13, one_third
       use num_lib
       use utils_lib
       use auto_diff
       type(auto_diff_real_star_order1), intent(in) :: chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL
       character(len=*), intent(in) :: MLT_option
-      real(dp), intent(in) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param
+      real(dp), intent(in) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
       type(auto_diff_real_star_order1), intent(out) :: Gamma, gradT, Y_face, conv_vel, D
       integer, intent(out) :: mixing_type, ierr
 
       real(dp) :: ff1, ff2, ff3, ff4
       type(auto_diff_real_star_order1) :: &
          Q, omega, a0, ff4_omega2_plus_1, A_1, A_2, &
          A_numerator, A_denom, A, Bcubed, delta, Zeta, &
-         f, f0, f1, f2, radiative_conductivity, convective_conductivity
+         f, f0, f1, f2, radiative_conductivity, convective_conductivity, csound
       include 'formats'
       if (gradr > gradL) then
          ! Convection zone
@@ -153,6 +153,10 @@ subroutine calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey
 
          ! average convection velocity   C&G 14.86b
          conv_vel = mixing_length_alpha*sqrt(Q*P/(8d0*rho))*Gamma / A
+
+         ! convective velocity limiter
+         if (conv_vel%val > max_conv_vel) conv_vel%val = max_conv_vel
+
          D = conv_vel*Lambda/3d0     ! diffusion coefficient [cm^2/sec]
 
          !Zeta = pow3(Gamma)/Bcubed  ! C&G 14.80

turb/public/turb.f90

@@ -92,10 +92,10 @@ end subroutine set_thermohaline
    subroutine set_TDC( &
             conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, report, &
             mixing_type, scale, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, &
-            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, ierr)
+            scale_height, gradL, grada, conv_vel, D, Y_face, gradT, tdc_num_iters, max_conv_vel, ierr)
       use tdc
       use tdc_support
-      real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale
+      real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, alpha_TDC_DAMP, alpha_TDC_DAMPR, alpha_TDC_PtdVdt, dt, cgrav, m, scale, max_conv_vel
       type(auto_diff_real_star_order1), intent(in) :: &
          chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada
       logical, intent(in) :: report
@@ -106,7 +106,7 @@ subroutine set_TDC( &
       real(dp), parameter :: alpha_c = (1d0/2d0)*sqrt_2_div_3
       real(dp), parameter :: lower_bound_Z = -1d2
       real(dp), parameter :: upper_bound_Z = 1d2
-      real(dp), parameter :: eps = 1d-2  ! Threshold in logY for separating multiple solutions.
+      real(dp), parameter :: eps = 1d-2 ! Threshold in logY for separating multiple solutions.
       type(auto_diff_real_tdc) :: Zub, Zlb
       include 'formats'
 
@@ -117,7 +117,8 @@ subroutine set_TDC( &
       call set_MLT('Cox', mixing_length_alpha, 0d0, 0d0, &
                      chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
                      gradr, grada, gradL, &
-                     Gamma, gradT, Y_face, conv_vel, D, mixing_type, ierr)
+                     Gamma, gradT, Y_face, conv_vel, D, mixing_type,1d99, ierr)
+
 
       ! Pack TDC info
       info%report = report
@@ -145,6 +146,16 @@ subroutine set_TDC( &
       Zlb = lower_bound_Z
       call get_TDC_solution(info, scale, Zlb, Zub, conv_vel, Y_face, tdc_num_iters, ierr)
 
+      ! Cap conv_vel at max_conv_vel_div_csound*cs
+      if (conv_vel%val > max_conv_vel) then
+         conv_vel = max_conv_vel
+         ! if max_conv_vel = csound,
+         ! L = L0 * (gradL + Y) + c0 * Af * Y_env
+         ! L = L0 * (gradL + Y) + c0 * sqrt_2_div_3 * csound * (Gamma / (1 + Gamma)) * Y
+         ! L - L0 * gradL = Y * (L0 + c0 * sqrt_2_div_3 * csound * (Gamma / (1 + Gamma)))
+         Y_face = unconvert(info%L - info%L0 * info%gradL) / (unconvert(info%L0) + unconvert(info%c0) * sqrt_2_div_3 * max_conv_vel * (info%Gamma / (1d0 + info%Gamma)))
+      end if
+
       ! Unpack output
       gradT = Y_face + gradL
       D = conv_vel*scale_height*mixing_length_alpha/3d0     ! diffusion coefficient [cm^2/sec]
@@ -225,18 +236,20 @@ end subroutine set_semiconvection
    subroutine set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
                      chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
                      gradr, grada, gradL, &
-                     Gamma, gradT, Y_face, conv_vel, D, mixing_type, ierr)
+                     Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
       use mlt
       type(auto_diff_real_star_order1), intent(in) :: chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL
       character(len=*), intent(in) :: MLT_option
-      real(dp), intent(in) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param
+      real(dp), intent(in) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel
+
       type(auto_diff_real_star_order1), intent(out) :: Gamma, gradT, Y_face, conv_vel, D
       integer, intent(out) :: mixing_type, ierr
 
       call calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, &
                      chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, &
                      gradr, grada, gradL, &
-                     Gamma, gradT, Y_face, conv_vel, D, mixing_type, ierr)
+                     Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr)
+
    end subroutine set_MLT
 
 end module turb