2: TDC consistency and energy conservation

docs/source/changelog.rst

@@ -66,6 +66,30 @@ AMR does not currently support RSP2.
 Bug Fixes
 ---------
 
+Fixed the TDC eddy-viscosity boundary condition for cell-centered Riemann
+hydrodynamics. The outer boundary now has zero turbulent stress while the
+surface cell retains the force from its inner turbulent face.
+
+Fixed a radius-collocation inconsistency in TDC eddy viscosity for
+cell-centered Riemann hydrodynamics. The strain, viscous heating, and momentum
+force now use the same lagged cell-midpoint radius, preserving the tridiagonal
+velocity coupling and nonpositive discrete viscous momentum work.
+
+Fixed an ambiguity between the generic transport velocity ``conv_vel`` and
+the MLT/TDC turbulent velocity ``mlt_vc``. ``do1_mlt_eval`` now copies
+``conv_vel`` into ``mlt_vc`` only for ``convective_mixing``. This preserves
+semiconvective and thermohaline diffusion coefficients while preventing their
+diffusion-equivalent velocities from entering turbulent energy, turbulent
+pressure, or TDC momentum terms.
+
+Fixed local eddy-viscous energy accounting when the ``dedt`` energy equation
+explicitly includes radial kinetic energy. TDC now includes the midpoint
+mechanical work from its cell-centered ``u_flag`` acceleration or its
+half-cell ``v_flag`` accelerations in addition to viscous heating. RSP2 uses
+the same ``v_flag`` correction. Energy equations that omit ``dKE/dt``,
+including the time-centered ``P d(1/rho)`` form used by pulsation models,
+continue to include viscous heating alone.
+
 Important bug fix for ``r26.4.1`` identified by Emily Sandford and Louis Siebenaler: the ``lowT_Freedman11`` opacity option used ``[M/H]`` labels as the metal mass fraction when interpolating in ``Z``, resulting in incorrect opacities. We recommend users who use these low-temperature opacities, such as in planet models, update to the latest MESA version or employ the fixes in :ref:`the known bugs entry <freedman_lowt_z_bug>` and `gh-993 <https://github.com/MESAHub/mesa/pull/993>`_.
 
 The plasmon neutrino cooling rate used a hardcoded prefactor calculated with a Weinberg angle of 0.2319, while all other neutrino cooling processes used calculated prefactors taking the Weinberg angle as input, with default value 0.22290. Thus, modifying the value of the Weinberg angle resulted in changes to neutrino cooling processes except for the plasmon neutrinos. This affects all previous MESA versions, and was found and fixed by user Garv Chauhan, see :ref:`the known bugs entry <plasmon_weinberg_angle_bug>` and `gh-998 <https://github.com/MESAHub/mesa/pull/998>`_. Plasmon neutrinos now use the same Weinberg angle as all other processes and changing its value will affect the corresponding cooling rate. Changes to the plasmon neutrino prefactor for MESA's default Weinberg angle result in small numerical differences for stars where plasmon neutrino cooling is significant.

star/defaults/controls.defaults

@@ -2153,6 +2153,33 @@
     mixing_length_alpha = 2
 
 
+      ! harmonic_dissipation_length_beta
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! If positive, limit the MLT and TDC mixing length using
+
+      ! .. math::
+
+      !    \frac{1}{\Lambda}
+      !    =
+      !    \frac{1}{\alpha_{\rm MLT} H_P}
+      !    +
+      !    \frac{1}{\mathtt{harmonic\_dissipation\_length\_beta}\,r}.
+
+      ! The control is the dimensionless :math:`\beta_r` in Eq. (17) of
+      ! `Wuchterl & Feuchtinger (1998)`_; :math:`\beta_r r` is the radial
+      ! limiting length. Here :math:`H_P=P/(\rho g)` is the usual pressure
+      ! scale height. This option ignores ``alt_scale_height_flag`` when
+      ! constructing the MLT and TDC mixing length. Non-positive values
+      ! preserve the existing behavior.
+
+      ! .. _Wuchterl & Feuchtinger (1998): https://ui.adsabs.harvard.edu/abs/1998A%26A...340..419W
+
+      ! ::
+
+    harmonic_dissipation_length_beta = 0d0
+
+
       ! remove_small_D_limit
       ! ~~~~~~~~~~~~~~~~~~~~
 
@@ -2264,8 +2291,8 @@
 
       ! If false, then stick to the usual definition -- P/(g*rho).
       ! If true, use min of the usual and sound speed * hydro time scale, sqrt(P/G)/rho.
-      ! Note that the 'TDC' ``MLT_option`` does not respect the ``alt_scale_height`` option, and continues to use ``h = P / rho g``
-      ! even if that flag is set.
+      ! A positive ``harmonic_dissipation_length_beta`` overrides this option
+      ! when constructing the MLT and TDC mixing length.
 
       ! ::
 

star/private/ctrls_io.f90

@@ -144,7 +144,8 @@ module ctrls_io
     superad_reduction_gamma_inv_scale, superad_reduction_diff_grads_limit, superad_reduction_limit, &
     superad_reduction_use_turnover_limit, superad_reduction_turnover_limit_function, &
     make_gradr_sticky_in_solver_iters, min_logT_for_make_gradr_sticky_in_solver_iters, &
-    max_logT_for_mlt, thermohaline_coeff, thermohaline_option, mixing_length_alpha, remove_small_D_limit, &
+    max_logT_for_mlt, thermohaline_coeff, thermohaline_option, mixing_length_alpha, &
+    harmonic_dissipation_length_beta, remove_small_D_limit, &
     alt_scale_height_flag, Henyey_MLT_y_param, Henyey_MLT_nu_param, no_MLT_below_shock, mlt_make_surface_no_mixing, &
     MLT_option, mlt_use_rotation_correction, mlt_Pturb_factor, do_normalize_dqs_as_part_of_set_qs, &
     max_Y_for_burn_z_mix_region, max_X_for_burn_he_mix_region, &
@@ -1083,6 +1084,7 @@ subroutine store_controls(s)
  s% thermohaline_coeff = thermohaline_coeff
  s% thermohaline_option = thermohaline_option
  s% mixing_length_alpha = mixing_length_alpha
+ s% harmonic_dissipation_length_beta = harmonic_dissipation_length_beta
  s% remove_small_D_limit = remove_small_D_limit
  s% alt_scale_height_flag = alt_scale_height_flag
  s% Henyey_MLT_y_param = Henyey_MLT_y_param
@@ -2820,6 +2822,7 @@ subroutine set_controls_for_writing(s, ierr)
  thermohaline_coeff = s% thermohaline_coeff
  thermohaline_option = s% thermohaline_option
  mixing_length_alpha = s% mixing_length_alpha
+ harmonic_dissipation_length_beta = s% harmonic_dissipation_length_beta
  remove_small_D_limit = s% remove_small_D_limit
  alt_scale_height_flag = s% alt_scale_height_flag
  Henyey_MLT_y_param = s% Henyey_MLT_y_param

star/private/hydro_energy.f90

@@ -74,7 +74,7 @@ subroutine get1_energy_eqn( &
          real(dp), dimension(s% species) :: &
             d_dwork_dxam1, d_dwork_dxa00, d_dwork_dxap1
          integer :: nz, i_dlnE_dt, i_lum, i_v
-         logical :: test_partials, doing_op_split_burn, eps_grav_form
+         logical :: test_partials, doing_op_split_burn, eps_grav_form, include_dke_dt
 
          include 'formats'
 
@@ -85,6 +85,10 @@ subroutine get1_energy_eqn( &
          call init
 
          call setup_eps_grav(ierr); if (ierr /= 0) return  ! do this first - it sets eps_grav_form
+         ! Only total-energy forms need work from momentum sources.
+         include_dke_dt = .not. eps_grav_form .and. &
+            .not. (s% using_velocity_time_centering .and. &
+               s% use_P_d_1_div_rho_form_of_work_when_time_centering_velocity)
          call setup_de_dt_and_friends(ierr); if (ierr /= 0) return
          call setup_dwork_dm(ierr); if (ierr /= 0) return
          call setup_dL_dm(ierr); if (ierr /= 0) return
@@ -229,14 +233,18 @@ end subroutine setup_dL_dm
 
 
          subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
-            use hydro_rsp2, only: compute_Eq_cell
-            use tdc_hydro, only: compute_tdc_Eq_div_w_face
+            use hydro_rsp2, only: compute_Eq_cell, compute_Uq_face
+            use tdc_hydro, only: &
+               compute_tdc_Eq_div_w_face, compute_tdc_Uq_face, compute_tdc_Uq_dm_cell
             real(dp) :: alfa, beta
             integer, intent(out) :: ierr
             type(auto_diff_real_star_order1) :: &
-               eps_nuc_ad, non_nuc_neu_ad, extra_heat_ad, Eq_ad, RTI_diffusion_ad, &
+               eps_nuc_ad, non_nuc_neu_ad, extra_heat_ad, Eq_ad, viscous_work_ad, &
+               Uq_00, Uq_p1, RTI_diffusion_ad, &
                v_00, v_p1, drag_force, drag_energy
             type(accurate_auto_diff_real_star_order1) :: sources_sum_ad
+            real(dp) :: kinetic_mass_factor
+            logical :: have_v_viscous_work
             include 'formats'
             ierr = 0
 
@@ -277,9 +285,22 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
                others_ad%val = others_ad%val + s% eps_phase_separation(k)
 
             Eq_ad = 0d0
+            viscous_work_ad = 0d0
+            Uq_00 = 0d0
+            Uq_p1 = 0d0
+            have_v_viscous_work = .false.
             if (s% RSP2_flag) then
                Eq_ad = s% Eq_ad(k)  ! compute_Eq_cell(s, k, ierr)
                if (ierr /= 0) return
+               if (include_dke_dt .and. s% v_flag) then
+                  Uq_00 = compute_Uq_face(s, k, ierr)
+                  if (ierr /= 0) return
+                  if (k < s% nz) then
+                     Uq_p1 = shift_p1(compute_Uq_face(s, k+1, ierr))
+                     if (ierr /= 0) return
+                  end if
+                  have_v_viscous_work = .true.
+               end if
             else if (s% TDC_alpha_M >0d0 .and. s% MLT_option == 'TDC' .and. &
                s% TDC_include_eturb_in_energy_equation .and. (s% v_flag .or. s% u_flag)) then
                 if (k < s% nz) then
@@ -289,10 +310,38 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
                      Eq_ad = 0.5d0*compute_tdc_Eq_div_w_face(s, k, ierr)*s% mlt_vc_ad(k)/sqrt_2_div_3
                 end if
                 if (ierr /= 0) return
+                if (include_dke_dt) then
+                   if (s% u_flag) then
+                      Uq_00 = compute_tdc_Uq_dm_cell(s, k, ierr)/s% dm(k)
+                      if (ierr /= 0) return
+                      v_00 = 0.5d0*(wrap_u_00(s,k) + s% u_start(k))
+                      viscous_work_ad = v_00*Uq_00
+                   else if (s% v_flag) then
+                      Uq_00 = compute_tdc_Uq_face(s, k, ierr)
+                      if (ierr /= 0) return
+                      if (k < s% nz) then
+                         Uq_p1 = shift_p1(compute_tdc_Uq_face(s, k+1, ierr))
+                         if (ierr /= 0) return
+                      end if
+                      have_v_viscous_work = .true.
+                   end if
+                end if
             end if
 
-            call setup_RTI_diffusion(RTI_diffusion_ad)
+            if (have_v_viscous_work) then
+               v_00 = 0.5d0*(wrap_v_00(s,k) + s% v_start(k))
+               v_p1 = 0d0
+               if (k < s% nz) &
+                  v_p1 = 0.5d0*(wrap_v_p1(s,k) + s% v_start(k+1))
+               kinetic_mass_factor = 1d0
+               if (s% use_mass_corrections) &
+                  kinetic_mass_factor = s% mass_correction(k)
+               ! Match the half-cell kinetic-energy quadrature.
+               viscous_work_ad = 0.5d0*kinetic_mass_factor* &
+                  (v_00*Uq_00 + v_p1*Uq_p1)
+            end if
 
+            call setup_RTI_diffusion(RTI_diffusion_ad)
             drag_energy = 0d0
             s% FdotV_drag_energy(k) = 0
             if (k /= s% nz) then
@@ -320,6 +369,7 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
             sources_sum_ad = sources_sum_ad - non_nuc_neu_ad
             sources_sum_ad = sources_sum_ad + extra_heat_ad
             sources_sum_ad = sources_sum_ad + Eq_ad
+            sources_sum_ad = sources_sum_ad + viscous_work_ad
             sources_sum_ad = sources_sum_ad + RTI_diffusion_ad
             sources_sum_ad = sources_sum_ad + drag_energy
             sources_ad = sources_sum_ad

star/private/mix_info.f90

@@ -60,7 +60,7 @@ subroutine set_mixing_info(s, skip_set_cz_bdy_mass, ierr)
 
          integer :: nz, k, max_conv_bdy, max_mix_bdy, k_Tmax, i_h1, i_he4, i_c12
          real(dp) :: rho_face, f, Tmax, min_conv_vel_for_convective_mixing_type, &
-            region_bottom_q, region_top_q, L_val
+            region_bottom_q, region_top_q, L_val, mixing_length
          real(dp), allocatable, dimension(:) :: eps_h, eps_he, eps_z, cdc_factor
 
          logical :: RSP2_or_RSP
@@ -362,9 +362,13 @@ subroutine set_mixing_info(s, skip_set_cz_bdy_mass, ierr)
          ! as last thing, update conv_vel from D_mix and mixing length.
          ! this updates the effective conv vel for rotation and overshooting effects
          do k=2,nz
-            if (s% alpha_mlt(k)*s% scale_height(k) > 0) then
-               s% conv_vel(k) = &
-                  3d0*s% D_mix(k)/(s% alpha_mlt(k)*s% scale_height(k))
+            if (s% harmonic_dissipation_length_beta > 0d0 .and. .not. RSP2_or_RSP) then
+               mixing_length = s% mlt_mixing_length(k)
+            else
+               mixing_length = s% alpha_mlt(k)*s% scale_height(k)
+            end if
+            if (mixing_length > 0d0) then
+               s% conv_vel(k) = 3d0*s% D_mix(k)/mixing_length
             else
                s% conv_vel(k) = 0
             end if
@@ -1017,9 +1021,17 @@ subroutine remove_mixing_singletons(s, ierr)
                if (s% cdc(k-1) /= 0 .and. s% cdc(k+1) /= 0) then
                   s% cdc(k) = (s% cdc(k-1) + s% cdc(k+1))/2
                   s% D_mix(k) = s% cdc(k)/pow2(pi4*s% r(k)*s% r(k)*s% rho(k))
-                  lambda = s% alpha_mlt(k)* &
-                     (s% scale_height(k-1) + s% scale_height(k+1))/2
-                  s% conv_vel(k) = 3*s% D_mix(k)/lambda
+                  if (s% harmonic_dissipation_length_beta > 0d0) then
+                     lambda = (s% mlt_mixing_length(k-1) + s% mlt_mixing_length(k+1))/2
+                  else
+                     lambda = s% alpha_mlt(k)* &
+                        (s% scale_height(k-1) + s% scale_height(k+1))/2
+                  end if
+                  if (lambda > 0d0) then
+                     s% conv_vel(k) = 3d0*s% D_mix(k)/lambda
+                  else
+                     s% conv_vel(k) = 0d0
+                  end if
                   s% mixing_type(k) = max(s% mixing_type(k-1), s% mixing_type(k+1))
                if (dbg) write(*,3) 'remove radiative singleton', k, nz
                end if

star/private/reconstructed_face_support.f90

@@ -30,9 +30,38 @@ module reconstructed_face_support
   public :: get_reconstructed_face_state_ad
   public :: get_reconstructed_face_eos_kap_ad
   public :: get_reconstructed_scale_height_ad
+  public :: get_reconstructed_hse_scale_height_ad
+  public :: get_effective_gradr_factor_ad
+  public :: get_Lrad_per_gradT_face_ad
 
 contains
 
+  function get_effective_gradr_factor_ad(s, k) result(gradr_factor)
+    type(star_info), pointer :: s
+    integer, intent(in) :: k
+    type(auto_diff_real_star_order1) :: gradr_factor
+
+    if (s%rotation_flag .and. s%mlt_use_rotation_correction) then
+       gradr_factor = s%ft_rot(k)/s%fp_rot(k)*s%gradr_factor(k)
+    else
+       gradr_factor = s%gradr_factor(k)
+    end if
+  end function get_effective_gradr_factor_ad
+
+
+  function get_Lrad_per_gradT_face_ad( &
+       s, k, T_face, P_face, opacity_face, gradr_factor) result(L0)
+    type(star_info), pointer :: s
+    integer, intent(in) :: k
+    type(auto_diff_real_star_order1), intent(in) :: &
+       T_face, P_face, opacity_face, gradr_factor
+    type(auto_diff_real_star_order1) :: L0, Pr_face
+
+    Pr_face = crad*pow4(T_face)/3d0
+    L0 = 4d0*pi4*clight*s%m_grav(k)*s%cgrav(k)*Pr_face/ &
+       (P_face*opacity_face*gradr_factor)
+  end function get_Lrad_per_gradT_face_ad
+
   ! Returns the MLT/TDC face thermodynamic state as
   ! auto_diff_real_star_order1 quantities, either from recomputed face
   ! EOS/opacity data or from the stored face quantities.
@@ -390,6 +419,39 @@ subroutine get_reconstructed_scale_height_ad(s, k, scale_height_face, ierr)
   end subroutine get_reconstructed_scale_height_ad
 
 
+  subroutine get_reconstructed_hse_scale_height_ad(s, k, scale_height_face, ierr)
+    use auto_diff_support, only: wrap_r_00
+    use star_utils, only: get_Peos_face, get_rho_face
+
+    type(star_info), pointer :: s
+    integer, intent(in) :: k
+    type(auto_diff_real_star_order1), intent(out) :: scale_height_face
+    integer, intent(out) :: ierr
+
+    real(dp) :: G
+    type(auto_diff_real_star_order1) :: P_face, rho_face, grav
+
+    ierr = 0
+    if (s%use_face_reconstruction) then
+       call ensure_reconstructed_face_state_ad(s, k, ierr)
+       if (ierr /= 0) return
+       P_face = s%reconstructed_P_face_ad(k)
+       rho_face = s%reconstructed_rho_face_ad(k)
+    else
+       P_face = get_Peos_face(s, k)
+       rho_face = get_rho_face(s, k)
+    end if
+
+    G = s%cgrav(k)
+    if (G <= 0d0 .or. s%m_grav(k) <= 0d0 .or. rho_face <= 0d0) then
+       scale_height_face = 0d0
+       return
+    end if
+    grav = G*s%m_grav(k)/pow2(wrap_r_00(s,k))
+    scale_height_face = P_face/(grav*rho_face)
+  end subroutine get_reconstructed_hse_scale_height_ad
+
+
   subroutine set_scale_height_from_face_state(s, k, P_face, rho_face, scale_height_face)
     use auto_diff_support, only: wrap_r_00
 

star/private/star_utils.f90

@@ -96,6 +96,7 @@ module star_utils
       public :: get_grada_face
       public :: get_gradr_face
       public :: get_scale_height_face
+      public :: get_mlt_mixing_length
       public :: get_tau
       public :: after_c_burn
       public :: get_shock_info
@@ -4032,6 +4033,30 @@ function get_scale_height_face(s,k) result(scale_height)
       end function get_scale_height_face
 
 
+      function get_mlt_mixing_length(s, pressure_scale_height, r, mixing_length_alpha) result(Lambda)
+         type (star_info), pointer :: s
+         type(auto_diff_real_star_order1), intent(in) :: pressure_scale_height, r
+         real(dp), intent(in) :: mixing_length_alpha
+         type(auto_diff_real_star_order1) :: Lambda, Lambda0, radial_length
+
+         Lambda0 = mixing_length_alpha*pressure_scale_height
+         if (s% harmonic_dissipation_length_beta <= 0d0) then
+            Lambda = Lambda0
+            return
+         end if
+
+         ! Wuchterl & Feuchtinger (1998), Eq. (17).
+         radial_length = s% harmonic_dissipation_length_beta*r
+         if (Lambda0 <= 0d0 .or. radial_length <= 0d0) then
+            Lambda = 0d0
+         else if (Lambda0 <= radial_length) then
+            Lambda = Lambda0/(1d0 + Lambda0/radial_length)
+         else
+            Lambda = radial_length/(1d0 + radial_length/Lambda0)
+         end if
+      end function get_mlt_mixing_length
+
+
       real(dp) function get_scale_height_face_val(s,k) result(scale_height)
          type (star_info), pointer :: s
          integer, intent(in) :: k