From 748422000c86342b0be342e4fe36254590bfcf55 Mon Sep 17 00:00:00 2001 From: Ebraheem Farag <63124736+Debraheem@users.noreply.github.com> Date: Wed, 26 Aug 2026 20:45:43 -0400 Subject: [PATCH] 2: TDC consistency and energy conservation --- docs/source/changelog.rst | 24 ++++ star/defaults/controls.defaults | 31 ++++- star/private/ctrls_io.f90 | 5 +- star/private/hydro_energy.f90 | 60 ++++++++- star/private/mix_info.f90 | 26 ++-- star/private/reconstructed_face_support.f90 | 62 +++++++++ star/private/star_utils.f90 | 25 ++++ star/private/tdc_hydro.f90 | 132 ++++++++++---------- star/private/turb_info.f90 | 52 +++++--- star/private/turb_support.f90 | 65 +++++++--- star_data/private/star_controls.inc | 2 +- star_data/public/star_data_step_work.inc | 2 +- turb/private/mlt.f90 | 5 +- turb/private/tdc.f90 | 2 +- turb/private/tdc_support.f90 | 14 ++- turb/public/turb.f90 | 18 +-- turb/test/src/test_turb.f90 | 19 +-- 17 files changed, 405 insertions(+), 139 deletions(-) diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index ca38cb36fa..91d44b989f 100644 --- a/docs/source/changelog.rst +++ b/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 ` and `gh-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 ` and `gh-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. diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 7f051d8511..2ba166820c 100644 --- a/star/defaults/controls.defaults +++ b/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. ! :: diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 5769321d92..afed2120ae 100644 --- a/star/private/ctrls_io.f90 +++ b/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 diff --git a/star/private/hydro_energy.f90 b/star/private/hydro_energy.f90 index 1dc228af38..5257682a4a 100644 --- a/star/private/hydro_energy.f90 +++ b/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 diff --git a/star/private/mix_info.f90 b/star/private/mix_info.f90 index f155e2eb35..e7e351ae2f 100644 --- a/star/private/mix_info.f90 +++ b/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 diff --git a/star/private/reconstructed_face_support.f90 b/star/private/reconstructed_face_support.f90 index ac73e3c55b..7a224da04c 100644 --- a/star/private/reconstructed_face_support.f90 +++ b/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 diff --git a/star/private/star_utils.f90 b/star/private/star_utils.f90 index 4ad8f8ef57..f6d29cc142 100644 --- a/star/private/star_utils.f90 +++ b/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 diff --git a/star/private/tdc_hydro.f90 b/star/private/tdc_hydro.f90 index 5bb50f0dac..04d4c98c0f 100644 --- a/star/private/tdc_hydro.f90 +++ b/star/private/tdc_hydro.f90 @@ -25,7 +25,8 @@ module tdc_hydro use auto_diff use auto_diff_support use star_utils - use reconstructed_face_support, only: get_reconstructed_scale_height_ad + use reconstructed_face_support, only: & + get_reconstructed_scale_height_ad, get_reconstructed_hse_scale_height_ad implicit none @@ -135,53 +136,42 @@ function get_TDC_Hp_face(s, k, ierr) result(Hp_face) end function get_TDC_Hp_face - function wrap_Hp_cell(s, k, ierr) result(Hp_cell) ! cm , different than rsp2 + function get_TDC_mixing_length_face(s, k, ierr) result(Lambda_face) type(star_info), pointer :: s integer, intent(in) :: k integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: Hp1, Hp0, Hp_cell + type(auto_diff_real_star_order1) :: Lambda_face, Hp_face + ierr = 0 - Hp0 = get_TDC_Hp_face(s, k, ierr) - if (ierr /= 0) return - Hp1 = 0d0 - if (k+1 < s%nz) then - Hp1 = shift_p1(get_TDC_Hp_face(s, k+1, ierr)) + if (s% harmonic_dissipation_length_beta > 0d0) then + call get_reconstructed_hse_scale_height_ad(s, k, Hp_face, ierr) + if (ierr /= 0) return + else + Hp_face = get_TDC_Hp_face(s, k, ierr) if (ierr /= 0) return end if - Hp_cell = 0.5d0*(Hp0 + Hp1) - !0.5d0*(wrap_Hp_00(s, k) + wrap_Hp_p1(s, k)) - end function wrap_Hp_cell + Lambda_face = get_mlt_mixing_length( & + s, Hp_face, wrap_r_00(s,k), s%mixing_length_alpha) + end function get_TDC_mixing_length_face + - function Hp_cell_for_Chi(s, k, ierr) result(Hp_cell) ! cm + function get_TDC_mixing_length_cell(s, k, ierr) result(Lambda_cell) type(star_info), pointer :: s integer, intent(in) :: k integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: Hp_cell - type(auto_diff_real_star_order1) :: d_00, Peos_00, rmid - real(dp) :: mmid, cgrav_mid - include 'formats' - ierr = 0 + type(auto_diff_real_star_order1) :: Lambda0, Lambda1, Lambda_cell - Hp_cell = wrap_Hp_cell(s, k, ierr) + ierr = 0 + Lambda0 = get_TDC_mixing_length_face(s, k, ierr) if (ierr /= 0) return - return ! below is skipped, for now. - - d_00 = wrap_d_00(s, k) - Peos_00 = wrap_Peos_00(s, k) - if (k < s%nz) then - rmid = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k)) - mmid = 0.5d0*(s%m(k) + s%m(k + 1)) - cgrav_mid = 0.5d0*(s%cgrav(k) + s%cgrav(k + 1)) - else - rmid = 0.5d0*(wrap_r_00(s, k) + s%r_center) - mmid = 0.5d0*(s%m(k) + s%m_center) - cgrav_mid = s%cgrav(k) - end if - Hp_cell = pow2(rmid)*Peos_00/(d_00*cgrav_mid*mmid) - if (s%alt_scale_height_flag) then - call mesa_error(__FILE__, __LINE__, 'Hp_cell_for_Chi: cannot use alt_scale_height_flag') + Lambda1 = 0d0 + if (k+1 < s%nz) then + Lambda1 = shift_p1(get_TDC_mixing_length_face(s, k+1, ierr)) + if (ierr /= 0) return end if - end function Hp_cell_for_Chi + Lambda_cell = 0.5d0*(Lambda0 + Lambda1) + end function get_TDC_mixing_length_cell + ! this function is only called internally in TDC_Uq_face, and for v_flag only. function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! does not update s% Chi or Chi_ad @@ -191,8 +181,8 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! does not update s% Chi type(auto_diff_real_star_order1) :: Chi_cell integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: & - rho2, r6_cell, d_v_div_r, Hp_cell, w_00, d_00, r_00, r_p1 - real(dp) :: f, ALFAM_ALFA + rho2, r6_cell, d_v_div_r, Lambda_cell, w_00, d_00, r_00, r_p1 + real(dp) :: f, ALFAM logical :: dbg include 'formats' ierr = 0 @@ -200,17 +190,17 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! does not update s% Chi ! check where we are getting alfam from. if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then - ALFAM_ALFA = s%TDC_alpha_M*s%mixing_length_alpha + ALFAM = s%TDC_alpha_M else ! this is for safety, but probably is never called. - ALFAM_ALFA = 0d0 + ALFAM = 0d0 end if - if (ALFAM_ALFA == 0d0 .or. & + if (ALFAM == 0d0 .or. & k <= s% TDC_num_outermost_cells_forced_nonturbulent .or. & k > s% nz - s% TDC_num_innermost_cells_forced_nonturbulent) then Chi_cell = 0d0 else - Hp_cell = Hp_cell_for_Chi(s, k, ierr) + Lambda_cell = get_TDC_mixing_length_cell(s, k, ierr) if (ierr /= 0) return d_v_div_r = compute_d_v_div_r(s, k, ierr) if (ierr /= 0) return @@ -232,12 +222,12 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! does not update s% Chi end if end if d_00 = wrap_d_00(s, k) - f = (16d0/3d0)*pi*ALFAM_ALFA/s%dm(k) + f = (16d0/3d0)*pi*ALFAM/s%dm(k) rho2 = pow2(d_00) r_00 = wrap_r_00(s, k) r_p1 = wrap_r_p1(s, k) r6_cell = 0.5d0*(pow6(r_00) + pow6(r_p1)) - Chi_cell = f*rho2*r6_cell*d_v_div_r*Hp_cell*w_00 + Chi_cell = f*rho2*r6_cell*d_v_div_r*Lambda_cell*w_00 ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm ! = g cm^2 s^-2 ! = erg @@ -248,8 +238,8 @@ function compute_Chi_cell(s, k, ierr) result(Chi_cell) ! does not update s% Chi !s%Chi_ad(k) = Chi_cell if (dbg .and. k == -100) then - write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA - write (*, *) 'Hp_cell', Hp_cell%val + write (*, *) ' s% ALFAM', ALFAM + write (*, *) 'Lambda_cell', Lambda_cell%val write (*, *) 'd_v_div_r', d_v_div_r%val write (*, *) ' f', f write (*, *) 'w_00', w_00%val @@ -269,8 +259,8 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) type(auto_diff_real_star_order1) :: Chi_face integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: & - rho2, r6_face, d_v_div_r, Hp_face, w_00, d_00, r_00, r_p1 - real(dp) :: f, ALFAM_ALFA, dmbar + rho2, r6_face, d_v_div_r, Lambda_face, w_00, d_00, r_00, r_p1 + real(dp) :: f, ALFAM, dmbar logical :: dbg include 'formats' ierr = 0 @@ -278,16 +268,17 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) ! check where we are getting alfam from. if (s%MLT_option == 'TDC' .and. .not. s%RSP2_flag) then - ALFAM_ALFA = s%TDC_alpha_M*s%mixing_length_alpha + ALFAM = s%TDC_alpha_M else ! this is for safety, but probably is never called. - ALFAM_ALFA = 0d0 + ALFAM = 0d0 end if - if (ALFAM_ALFA == 0d0 .or. & + if (ALFAM == 0d0 .or. (s% u_flag .and. k == 1) .or. & k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then + ! The outer boundary is stress free for cell-centered u hydro. Chi_face = 0d0 else - Hp_face = get_TDC_Hp_face(s, k, ierr) + Lambda_face = get_TDC_mixing_length_face(s, k, ierr) if (ierr /= 0) return d_v_div_r = compute_d_v_div_r_face(s, k, ierr) if (ierr /= 0) return @@ -298,12 +289,12 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) dmbar = 0.5d0*s% dm(k) end if d_00 = get_rho_face(s, k) - f = (16d0/3d0)*pi*ALFAM_ALFA/dmbar + f = (16d0/3d0)*pi*ALFAM/dmbar rho2 = pow2(d_00) r_00 = wrap_r_00(s, k) !r_p1 = wrap_r_p1(s, k) r6_face = pow6(r_00) !0.5d0*(pow6(r_00) + pow6(r_p1)) - Chi_face = f*rho2*r6_face*d_v_div_r*Hp_face!*w_00 + Chi_face = f*rho2*r6_face*d_v_div_r*Lambda_face!*w_00 ! units = g^-1 cm s^-1 g^2 cm^-6 cm^6 s^-1 cm * [s/cm] ! [1/w_00] = [s/cm] ! = g cm^2 s^-2 * [s/cm] ! = erg ! * [s / cm] - > [erg] * [s/cm] @@ -321,8 +312,8 @@ function compute_Chi_div_w_face(s, k, ierr) result(Chi_face) s%Chi_ad(k) = Chi_face*w_00 if (dbg .and. k == -100) then - write (*, *) ' s% ALFAM_ALFA', ALFAM_ALFA - write (*, *) 'Hp_face', Hp_face%val + write (*, *) ' s% ALFAM', ALFAM + write (*, *) 'Lambda_face', Lambda_face%val write (*, *) 'd_v_div_r', d_v_div_r%val write (*, *) ' f', f write (*, *) 'w_00', w_00%val @@ -417,7 +408,7 @@ function compute_tdc_Uq_dm_cell(s, k, ierr) result(Uq_cell) ! cm s^-2, accelera type(star_info), pointer :: s integer, intent(in) :: k integer, intent(out) :: ierr - type(auto_diff_real_star_order1) :: Chi_00, Chi_p1, r_00, r_p1, w_00, w_p1, r_cell, Uq_cell + type(auto_diff_real_star_order1) :: Chi_00, Chi_p1, w_00, w_p1, r_cell, Uq_cell include 'formats' ierr = 0 if (s%mixing_length_alpha == 0d0 .or. & @@ -425,9 +416,8 @@ function compute_tdc_Uq_dm_cell(s, k, ierr) result(Uq_cell) ! cm s^-2, accelera k > s%nz - s% TDC_num_innermost_cells_forced_nonturbulent) then Uq_cell = 0d0 else - r_00 = wrap_opt_time_center_r_00(s, k) - r_p1 = wrap_opt_time_center_r_p1(s, k) - r_cell = 0.5d0*(r_00+r_p1) ! not staggered unlike terms inside chi_div_w_face + ! Use the same lagged cell radius as the u/r strain. + r_cell = s% rmid_start(k) if (s% okay_to_set_mlt_vc .and. & s% TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation) then @@ -458,7 +448,7 @@ function compute_tdc_Uq_dm_cell(s, k, ierr) result(Uq_cell) ! cm s^-2, accelera if (k == -56) then write (*, 3) 'TDC Uq chi_m1 chi_00 r', k, s%solver_iter, & - Uq_cell%val, Chi_p1%val, Chi_00%val, r_00%val + Uq_cell%val, Chi_p1%val, Chi_00%val, r_cell%val end if end if @@ -531,9 +521,14 @@ function compute_d_v_div_r_face(s, k, ierr) result(d_v_div_r) ! s^-1 if (s% v_flag) then r_00 = 0.5d0*(wrap_r_00(s, k) + wrap_r_p1(s, k)) r_m1 = 0.5d0*(wrap_r_00(s, k) + wrap_r_m1(s, k)) - else if(s% u_flag) then ! stagger r for u_flag to retain tridiagonality. - r_00 = wrap_r_00(s, k) - r_m1 = wrap_r_m1(s, k) + else if(s% u_flag) then + ! Lag the cell-centered radius to retain tridiagonality. + r_00 = s% rmid_start(k) + if (k > 1) then + r_m1 = s% rmid_start(k-1) + else + r_m1 = 1d0 + end if end if if (r_00%val == 0d0) r_00 = 1d0 @@ -571,9 +566,14 @@ function compute_d_v_div_r_opt_time_center_face(s, k, ierr) result(d_v_div_r) ! if (s% v_flag) then r_00 = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_p1(s, k)) r_m1 = 0.5d0*(wrap_opt_time_center_r_00(s, k) + wrap_opt_time_center_r_m1(s, k)) - else if(s% u_flag) then ! stagger r for u_flag to retain tridiagonality. - r_00 = wrap_opt_time_center_r_00(s, k) - r_m1 = wrap_opt_time_center_r_m1(s, k) + else if(s% u_flag) then + ! Lag the cell-centered radius to retain tridiagonality. + r_00 = s% rmid_start(k) + if (k > 1) then + r_m1 = s% rmid_start(k-1) + else + r_m1 = 1d0 + end if end if if (r_00%val == 0d0) r_00 = 1d0 diff --git a/star/private/turb_info.f90 b/star/private/turb_info.f90 index 4c98d2684f..27a21a5ef5 100644 --- a/star/private/turb_info.f90 +++ b/star/private/turb_info.f90 @@ -22,7 +22,8 @@ module turb_info use star_private_def use const_def, only: dp, i8, ln10, pi4, no_mixing, convective_mixing, crystallized, phase_separation_mixing - use reconstructed_face_support, only: get_reconstructed_face_state_ad + use reconstructed_face_support, only: get_reconstructed_face_state_ad, get_reconstructed_hse_scale_height_ad, & + get_effective_gradr_factor_ad, get_Lrad_per_gradT_face_ad use num_lib use utils_lib use auto_diff_support @@ -111,7 +112,7 @@ subroutine do1_mlt_2(s, k, & type(auto_diff_real_star_order1) :: & T_face_ad, P_face_ad, energy_face_ad, opacity_face_ad, rho_face_ad, chiRho_face_ad, chiT_face_ad, Cp_face_ad, & grada_face_ad, scale_height_ad, gradr_ad, & - gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad + gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, mixing_length_ad, hse_scale_height_ad, L0_ad, Lrad_ad include 'formats' ierr = 0 @@ -164,11 +165,16 @@ subroutine do1_mlt_2(s, k, & s% tau_conv_start(k) = tau_conv end if - 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) + if (s% harmonic_dissipation_length_beta > 0d0) then + call get_reconstructed_hse_scale_height_ad(s, k, hse_scale_height_ad, ierr) + if (ierr /= 0) return + mixing_length_ad = get_mlt_mixing_length( & + s, hse_scale_height_ad, wrap_r_00(s,k), mixing_length_alpha) else - gradr_factor = s% gradr_factor(k) + mixing_length_ad = mixing_length_alpha*scale_height_ad end if + + gradr_factor = get_effective_gradr_factor_ad(s, k) if (is_bad_num(gradr_factor% val)) then ierr = -1 return @@ -267,7 +273,7 @@ subroutine do1_mlt_2(s, k, & call do1_mlt_eval(s, k, s% MLT_option, gradL_composition_term, & T_face_ad, P_face_ad, energy_face_ad, opacity_face_ad, rho_face_ad, chiRho_face_ad, chiT_face_ad, Cp_face_ad, & gradr_ad, grada_face_ad, scale_height_ad, mixing_length_alpha, & - mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, ierr) + mixing_type, gradT_ad, Y_face_ad, mlt_vc_ad, D_ad, Gamma_ad, mixing_length_ad, ierr) if (ierr /= 0) then if (s% report_ierr) then write(*,*) 'ierr in do1_mlt_eval for k', k @@ -284,10 +290,13 @@ subroutine do1_mlt_2(s, k, & end if call adjust_gradT_fraction(s, k, f) - if (s% mlt_mixing_type(k) == no_mixing .or. abs(s% gradr(k)) < 1d-20) then + if (s% mlt_mixing_type(k) == no_mixing .or. abs(gradr_factor%val) <= 1d-20) then s% L_conv(k) = 0d0 else - s% L_conv(k) = s% L(k) * (1d0 - s% gradT(k)/s% gradr(k)) ! C&G 14.109 + L0_ad = get_Lrad_per_gradT_face_ad( & + s, k, T_face_ad, P_face_ad, opacity_face_ad, gradr_factor) + Lrad_ad = L0_ad*s% gradT_ad(k) + s% L_conv(k) = s% L(k) - Lrad_ad%val ! C&G 14.109 end if contains @@ -325,7 +334,7 @@ subroutine store_results s% scale_height_ad(k) = scale_height_ad s% scale_height(k) = scale_height_ad%val - s% Lambda_ad(k) = mixing_length_alpha*scale_height_ad + s% Lambda_ad(k) = mixing_length_ad s% mlt_mixing_length(k) = s% Lambda_ad(k)%val end subroutine store_results @@ -366,7 +375,7 @@ subroutine set_no_mixing(str) s% scale_height_ad(k) = scale_height_ad s% scale_height(k) = scale_height_ad%val - s% Lambda_ad(k) = mixing_length_alpha*scale_height_ad + s% Lambda_ad(k) = mixing_length_ad s% mlt_mixing_length(k) = s% Lambda_ad(k)%val s% L_conv(k) = 0d0 @@ -607,7 +616,7 @@ subroutine check_for_redo_MLT(s, nzlo, nzhi, ierr) subroutine end_of_convective_region() integer :: kk, op_err - real(dp) :: Hp + real(dp) :: mixing_length logical :: end_dbg 9 format(a40, 3i7, 99(1pd26.16)) include 'formats' @@ -616,8 +625,12 @@ subroutine end_of_convective_region() top_r = s% r(k) top_Hp = s% scale_height(k) dr = top_r - bot_r - Hp = (bot_Hp + top_Hp)/2 - if (dr < s% alpha_mlt(k)*min(top_Hp, bot_Hp) .and. & + if (s% harmonic_dissipation_length_beta > 0d0) then + mixing_length = min(s% mlt_mixing_length(k), s% mlt_mixing_length(k_bot)) + else + mixing_length = s% alpha_mlt(k)*min(top_Hp, bot_Hp) + end if + if (dr < mixing_length .and. & s% redo_conv_for_dr_lt_mixing_length) then !$OMP PARALLEL DO PRIVATE(kk,op_err) SCHEDULE(dynamic,2) do kk = k, k_bot @@ -635,13 +648,24 @@ subroutine redo1_mlt(s, k, dr, ierr) real(dp), intent(in) :: dr integer, intent(out) :: ierr logical :: make_gradr_sticky_in_solver_iters + real(dp) :: mixing_length_alpha, radial_length + type(auto_diff_real_star_order1) :: hse_scale_height include 'formats' ierr = 0 if (dr >= s% mlt_mixing_length(k)) return ! if convection zone is smaller than mixing length ! redo MLT with reduced alpha so mixing_length = dr + if (s% harmonic_dissipation_length_beta > 0d0) then + call get_reconstructed_hse_scale_height_ad(s, k, hse_scale_height, ierr) + if (ierr /= 0) return + radial_length = s% harmonic_dissipation_length_beta*s% r(k) + if (hse_scale_height%val <= 0d0 .or. radial_length <= dr) return + mixing_length_alpha = dr*radial_length/(hse_scale_height%val*(radial_length - dr)) + else + mixing_length_alpha = dr/s% scale_height(k) + end if call do1_mlt_2(s, k, make_gradr_sticky_in_solver_iters, ierr, & - mixing_length_alpha_in = dr/s% scale_height(k)) + mixing_length_alpha_in = mixing_length_alpha) end subroutine redo1_mlt end subroutine check_for_redo_MLT diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90 index 6f0c3aab92..0c67fca4fa 100644 --- a/star/private/turb_support.f90 +++ b/star/private/turb_support.f90 @@ -113,9 +113,9 @@ subroutine do1_mlt_eval( & s, k, MLT_option, gradL_composition_term, & T_in, P_in, energy_in, opacity_in, rho_in, chiRho_in, chiT_in, Cp_in, & gradr_in, grada, scale_height, mixing_length_alpha, & - mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr) + mixing_type, gradT, Y_face, mlt_vc, D, Gamma, mixing_length, ierr) use chem_def, only: ih1 - use const_def, only: ln10 + use const_def, only: convective_mixing use starspots, only: starspot_tweak_gradr type (star_info), pointer :: s integer, intent(in) :: k @@ -126,13 +126,13 @@ subroutine do1_mlt_eval( & real(dp), intent(in) :: gradL_composition_term, mixing_length_alpha integer, intent(out) :: mixing_type type(auto_diff_real_star_order1), intent(out) :: & - gradT, Y_face, mlt_vc, D, Gamma + gradT, Y_face, mlt_vc, D, Gamma, mixing_length integer, intent(out) :: ierr real(dp) :: cgrav, m, XH1, P_theta, L_theta integer :: iso type(auto_diff_real_star_order1) :: gradr, r, L, T, P, opacity, rho, dV, & - chiRho, chiT, Cp, rho_start, energy + chiRho, chiT, Cp, rho_start, energy, conv_vel include 'formats' ierr = 0 @@ -182,11 +182,12 @@ subroutine do1_mlt_eval( & XH1 = s% xa(s% net_iso(ih1),k) if (s% use_other_mlt_results) then + mixing_length = mixing_length_alpha*scale_height call s% other_mlt_results(s% id, k, MLT_option, & r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & s% alpha_semiconvection, s% thermohaline_coeff, & - mixing_type, gradT, Y_face, mlt_vc, D, Gamma, energy, ierr) + mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr) else ! starspot YREC routine if (s% do_starspots) then @@ -197,7 +198,15 @@ subroutine do1_mlt_eval( & r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & s% alpha_semiconvection, s% thermohaline_coeff, & - mixing_type, gradT, Y_face, mlt_vc, D, Gamma, energy, ierr) + mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr, mixing_length) + end if + if (ierr /= 0) return + + ! Do not use diffusion velocities as MLT/TDC turbulent velocity. + if (mixing_type == convective_mixing) then + mlt_vc = conv_vel + else + mlt_vc = 0d0 end if end subroutine do1_mlt_eval @@ -207,7 +216,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, & iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, & alpha_semiconvection, thermohaline_coeff, & - mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr) + mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr, mixing_length) use star_utils use tdc_hydro, only: compute_tdc_Eq_div_w_face type (star_info), pointer :: s @@ -222,8 +231,10 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg integer, intent(out) :: mixing_type type(auto_diff_real_star_order1), intent(out) :: gradT, Y_face, conv_vel, D, Gamma integer, intent(out) :: ierr + type(auto_diff_real_star_order1), intent(out), optional :: mixing_length - type(auto_diff_real_star_order1) :: Pr, Pg, grav, Lambda, gradL, beta + type(auto_diff_real_star_order1) :: & + Pr, Pg, grav, Hp_for_mlt, Lambda, legacy_Lambda, mixing_length_alpha_ad, gradL, beta real(dp) :: conv_vel_start, scale, max_conv_vel, Y_face_guess ! these are used by use_superad_reduction @@ -262,7 +273,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg Pr = crad*pow4(T)/3d0 Pg = Ptot - Pr beta = Pg / Ptot - Lambda = mixing_length_alpha*scale_height + legacy_Lambda = mixing_length_alpha*scale_height if (k == 0) then grav = cgrav*m/pow2(r) @@ -270,6 +281,23 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg grav = cgrav*m/pow2(r) !try replacing with wrap_geff_face(s,k) end if + if (s% harmonic_dissipation_length_beta > 0d0) then + if (P > 0d0 .and. grav > 0d0 .and. rho > 0d0) then + Hp_for_mlt = P/(rho*grav) + Lambda = get_mlt_mixing_length(s, Hp_for_mlt, r, mixing_length_alpha) + mixing_length_alpha_ad = Lambda/Hp_for_mlt + else + Hp_for_mlt = 0d0 + Lambda = 0d0 + mixing_length_alpha_ad = 0d0 + end if + else + Hp_for_mlt = scale_height + Lambda = legacy_Lambda + mixing_length_alpha_ad = mixing_length_alpha + end if + if (present(mixing_length)) mixing_length = Lambda + if (s% use_Ledoux_criterion) then gradL = grada + gradL_composition_term ! Ledoux temperature gradient else @@ -356,10 +384,10 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg end if call set_TDC(& - conv_vel_start, mixing_length_alpha, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, & + conv_vel_start, mixing_length_alpha_ad, Lambda, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, & s%dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr, r, Ptot, T, rho, dV, Cp, opacity, & - scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, & + Hp_for_mlt, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, & Eq_div_w, grav, & s% include_mlt_corr_to_TDC, s% TDC_alpha_C, s% TDC_alpha_S, s% use_TDC_enthalpy_flux_limiter, energy, & Y_face_guess, ierr) @@ -378,10 +406,10 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg call set_superad_reduction if (Gamma_factor > 1d0) then call set_TDC(& - conv_vel_start, mixing_length_alpha, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, & + conv_vel_start, mixing_length_alpha_ad, Lambda, s%TDC_alpha_D, s%TDC_alpha_R, s%TDC_alpha_Pt, & s%dt, cgrav, m, report, & mixing_type, scale, chiT, chiRho, gradr_scaled, r, Ptot, T, rho, dV, Cp, opacity, & - scale_height, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, & + Hp_for_mlt, gradL, grada, conv_vel, D, Y_face, gradT, s%tdc_num_iters(k), max_conv_vel, & Eq_div_w, grav, & s% include_mlt_corr_to_TDC, s% TDC_alpha_C, s% TDC_alpha_S, s% use_TDC_enthalpy_flux_limiter, energy, & Y_face_guess, ierr) @@ -395,7 +423,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg else if (gradr > gradL) then if (report) write(*,3) 'call set_MLT', k, s% solver_iter - call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, & + call set_MLT(MLT_option, mixing_length_alpha_ad, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, & chiT, chiRho, Cp, grav, Lambda, rho, Ptot, T, opacity, & gradr, grada, gradL, & Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr) @@ -413,7 +441,7 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg if (s% use_superad_reduction) then call set_superad_reduction if (Gamma_factor > 1d0) then - call set_MLT(MLT_option, mixing_length_alpha, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, & + call set_MLT(MLT_option, mixing_length_alpha_ad, s% Henyey_MLT_nu_param, s% Henyey_MLT_y_param, & chiT, chiRho, Cp, grav, Lambda, rho, Ptot, T, opacity, & gradr_scaled, grada, gradL, & Gamma, gradT, Y_face, conv_vel, D, mixing_type, max_conv_vel, ierr) @@ -430,7 +458,9 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg if (mixing_type == no_mixing) then if (gradL_composition_term < 0) then if (report) write(*,3) 'call set_thermohaline', k, s% solver_iter - call set_thermohaline(s%thermohaline_option, Lambda, grada, gradr, T, opacity, rho, Cp, gradL_composition_term, & + if (present(mixing_length)) mixing_length = legacy_Lambda + call set_thermohaline(s%thermohaline_option, legacy_Lambda, grada, gradr, T, opacity, rho, Cp, & + gradL_composition_term, & iso, XH1, thermohaline_coeff, & D, gradT, Y_face, conv_vel, mixing_type, ierr) if (ierr /= 0) then @@ -439,7 +469,8 @@ subroutine Get_results(s, k, MLT_option, & ! NOTE: k=0 is a valid arg end if else if (gradr > grada) then if (report) write(*,3) 'call set_semiconvection', k, s% solver_iter - call set_semiconvection(L, Lambda, m, T, Ptot, Pr, beta, opacity, rho, alpha_semiconvection, & + if (present(mixing_length)) mixing_length = legacy_Lambda + call set_semiconvection(L, legacy_Lambda, m, T, Ptot, Pr, beta, opacity, rho, alpha_semiconvection, & s% semiconvection_option, cgrav, Cp, gradr, grada, gradL, & gradL_composition_term, & gradT, Y_face, conv_vel, D, mixing_type, ierr) diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index de0ded144a..a3f62b5b84 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -251,7 +251,7 @@ real(dp) :: thermohaline_coeff character (len=32) :: & thermohaline_option, MLT_option, semiconvection_option - real(dp) :: mixing_length_alpha, remove_small_D_limit + real(dp) :: mixing_length_alpha, harmonic_dissipation_length_beta, remove_small_D_limit logical :: use_Ledoux_criterion integer :: num_cells_for_smooth_gradL_composition_term integer :: steps_before_use_TDC diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index 4f9f176c03..59a1c7f2ec 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -361,7 +361,7 @@ ! NOTE: the following values are from the most recent call on mlt ! values for these at start of step have "_start" at end of name - real(dp), pointer :: mlt_mixing_length(:) ! mixing_length_alpha*scale_height as computed by mlt + real(dp), pointer :: mlt_mixing_length(:) ! mixing length from mlt (cm) integer, pointer :: mlt_mixing_type(:) ! as defined in mesa/const. before overshooting has been added. real(dp), pointer :: mlt_D(:) ! mlt diffusion coeff. (cm^2/sec) real(dp), pointer :: mlt_Gamma(:) ! local mlt convection efficiency. diff --git a/turb/private/mlt.f90 b/turb/private/mlt.f90 index fd5e0e9ef8..3ef7584d9e 100644 --- a/turb/private/mlt.f90 +++ b/turb/private/mlt.f90 @@ -58,9 +58,10 @@ subroutine calc_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, & 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 + type(auto_diff_real_star_order1), intent(in) :: & + mixing_length_alpha, 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, max_conv_vel + real(dp), intent(in) :: 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 diff --git a/turb/private/tdc.f90 b/turb/private/tdc.f90 index f4bd3c534f..d3b9f225c8 100644 --- a/turb/private/tdc.f90 +++ b/turb/private/tdc.f90 @@ -64,7 +64,7 @@ subroutine get_TDC_solution(info, scale, Zlb, Zub, conv_vel, Y_face, tdc_num_ite include 'formats' ierr = 0 - if (info%mixing_length_alpha == 0d0 .or. info%dt <= 0d0) then + if (info%Lambda <= 0d0 .or. info%dt <= 0d0) then call mesa_error(__FILE__,__LINE__,'bad call to TDC get_TDC_solution') end if diff --git a/turb/private/tdc_support.f90 b/turb/private/tdc_support.f90 index 99d02e5b76..6c4ca26598 100644 --- a/turb/private/tdc_support.f90 +++ b/turb/private/tdc_support.f90 @@ -43,7 +43,8 @@ module tdc_support !! do not depend on Y. !! !! @param report Write debug output if true, not if false. - !! @param mixing_length_alpha Mixing length parameter + !! @param mixing_length_alpha Effective mixing length parameter Lambda/Hp + !! @param Lambda Mixing length !! @param TDC_alpha_D TDC turbulent damping parameter !! @param TDC_alpha_R TDC radiative damping parameter !! @param TDC_alpha_Pt TDC coefficient on P_turb*dV/dt. Physically should probably be 1. @@ -63,9 +64,9 @@ module tdc_support !! @param Gamma Gamma is the MLT Gamma efficiency parameter, which we evaluate in steady state from MLT. type tdc_info logical :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter - real(dp) :: mixing_length_alpha, TDC_alpha_C, TDC_alpha_S, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, e + real(dp) :: TDC_alpha_C, TDC_alpha_S, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, e type(auto_diff_real_tdc) :: A0, c0, L, L0, gradL, grada - type(auto_diff_real_star_order1) :: T, rho, dV, Cp, kap, Hp, Gamma, Eq_div_w, P, h + type(auto_diff_real_star_order1) :: mixing_length_alpha, Lambda, T, rho, dV, Cp, kap, Hp, Gamma, Eq_div_w, P, h end type tdc_info contains @@ -534,7 +535,8 @@ subroutine eval_xis(info, Y, xi0, xi1, xi2) scale = 1d0 if (Y > 0d0 .and. info%use_TDC_enthalpy_flux_limiter) then ! X = G/F - X = convert(info%Cp*info%T/info%h)*info%mixing_length_alpha * info%TDC_alpha_S * x_ALFAS * Y / sqrt_2_div_3 + X = convert(info%Cp*info%T/info%h*info%mixing_length_alpha) * & + info%TDC_alpha_S * x_ALFAS * Y / sqrt_2_div_3 FL = flux_limiter_function(X) ! Avoid 0/0 or tiny/tiny; for X ≈ 0, FL ≈ X so scale ~ 1 anyway. if (abs(X%val) >= 0.95d0) then @@ -551,8 +553,8 @@ subroutine eval_xis(info, Y, xi0, xi1, xi2) S0 = S0*Y + convert(info%Eq_div_w) end if - D0 = convert(info%TDC_alpha_D*x_CEDE/(info%mixing_length_alpha*info%Hp)) - gammar_div_alfa = info%TDC_alpha_R*x_GAMMAR/(info%mixing_length_alpha*info%Hp) + D0 = convert(info%TDC_alpha_D*x_CEDE/info%Lambda) + gammar_div_alfa = info%TDC_alpha_R*x_GAMMAR/info%Lambda DR0 = convert(4d0*boltz_sigma*pow2(gammar_div_alfa)*pow3(info%T)/(pow2(info%rho)*info%Cp*info%kap)) Pt0 = info%TDC_alpha_Pt*x_ALFAP*info%rho dVdt = info%dV/info%dt diff --git a/turb/public/turb.f90 b/turb/public/turb.f90 index 20789e5552..a2d80697ec 100644 --- a/turb/public/turb.f90 +++ b/turb/public/turb.f90 @@ -113,23 +113,24 @@ end subroutine set_thermohaline !! @param ierr Tracks errors (output). !! @param Y_face_guess Candidate superadiabaticity for the local solve. Non-positive values disable seeding. subroutine set_TDC( & - conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, report, & + conv_vel_start, mixing_length_alpha, Lambda, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, 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, & max_conv_vel, Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, & TDC_alpha_S, use_TDC_enthalpy_flux_limiter, energy, Y_face_guess, ierr) use tdc use tdc_support - real(dp), intent(in) :: conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt + real(dp), intent(in) :: conv_vel_start, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt real(dp), intent(in) :: dt, cgrav, m, scale, max_conv_vel, TDC_alpha_C, TDC_alpha_S type(auto_diff_real_star_order1), intent(in) :: & - chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, scale_height, gradL, grada, Eq_div_w, grav, energy + mixing_length_alpha, Lambda, chiT, chiRho, gradr, r, P, T, rho, dV, Cp, opacity, & + scale_height, gradL, grada, Eq_div_w, grav, energy logical, intent(in) :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter real(dp), intent(in) :: Y_face_guess type(auto_diff_real_star_order1),intent(out) :: conv_vel, Y_face, gradT, D integer, intent(out) :: tdc_num_iters, mixing_type, ierr type(tdc_info) :: info - type(auto_diff_real_star_order1) :: L, Lambda, Gamma, h + type(auto_diff_real_star_order1) :: L, Gamma, h 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 @@ -140,7 +141,6 @@ subroutine set_TDC( & !grav = cgrav * m / pow2(r) L = 64d0 * pi * boltz_sigma * pow4(T) * grav * pow2(r) * gradr / (3d0 * P * opacity) if (include_mlt_corr_to_TDC) then - Lambda = mixing_length_alpha * scale_height call set_MLT('Cox', mixing_length_alpha, 0d0, 0d0, & chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, & gradr, grada, gradL, & @@ -161,6 +161,7 @@ subroutine set_TDC( & info%include_mlt_corr_to_TDC = include_mlt_corr_to_TDC info%use_TDC_enthalpy_flux_limiter = use_TDC_enthalpy_flux_limiter info%mixing_length_alpha = mixing_length_alpha + info%Lambda = Lambda info%TDC_alpha_D = TDC_alpha_D info%TDC_alpha_R = TDC_alpha_R info%TDC_alpha_Pt = TDC_alpha_Pt @@ -206,7 +207,7 @@ subroutine set_TDC( & ! Unpack output gradT = Y_face + gradL - D = conv_vel*scale_height*mixing_length_alpha/3d0 ! diffusion coefficient [cm^2/sec] + D = conv_vel*Lambda/3d0 ! diffusion coefficient [cm^2/sec] if (conv_vel > 0d0) then mixing_type = convective_mixing else @@ -289,9 +290,10 @@ subroutine set_MLT(MLT_option, mixing_length_alpha, Henyey_MLT_nu_param, Henyey_ gradr, grada, gradL, & 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 + type(auto_diff_real_star_order1), intent(in) :: & + mixing_length_alpha, 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, max_conv_vel + real(dp), intent(in) :: 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 diff --git a/turb/test/src/test_turb.f90 b/turb/test/src/test_turb.f90 index 00e500715b..51b5566b0a 100644 --- a/turb/test/src/test_turb.f90 +++ b/turb/test/src/test_turb.f90 @@ -25,9 +25,10 @@ subroutine header(text) end subroutine header subroutine check_efficient_MLT_scaling() - type(auto_diff_real_star_order1) :: chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL + type(auto_diff_real_star_order1) :: & + mixing_length_alpha, chiT, chiRho, Cp, grav, Lambda, rho, P, T, opacity, gradr, grada, gradL character(len=3) :: MLT_option - real(dp) :: mixing_length_alpha, Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel + real(dp) :: Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel type(auto_diff_real_star_order1) :: Gamma, gradT, Y_face, conv_vel, conv_vel2, D, r, L integer :: mixing_type, ierr @@ -80,10 +81,11 @@ subroutine check_efficient_MLT_scaling() end subroutine check_efficient_MLT_scaling subroutine compare_TDC_and_Cox_MLT() - real(dp) :: mixing_length_alpha, conv_vel_start, & + real(dp) :: conv_vel_start, & TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, scale, L_start, TDC_alpha_C, TDC_alpha_S type(auto_diff_real_star_order1) :: & - r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL, grav, Lambda + mixing_length_alpha, r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, & + scale_height, gradL, grav, Lambda type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Gamma, Eq_div_w, energy real(dp) :: Henyey_MLT_nu_param, Henyey_MLT_y_param, max_conv_vel @@ -147,7 +149,7 @@ subroutine compare_TDC_and_Cox_MLT() write (*, 1) 'gradR - gradA', gradr%val - grada%val call set_TDC( & - conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, report, & + conv_vel_start, mixing_length_alpha, Lambda, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, 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, max_conv_vel, & Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, TDC_alpha_S, use_TDC_enthalpy_flux_limiter, & @@ -166,10 +168,10 @@ subroutine compare_TDC_and_Cox_MLT() end subroutine compare_TDC_and_Cox_MLT subroutine check_TDC() - real(dp) :: mixing_length_alpha, conv_vel_start + real(dp) :: conv_vel_start real(dp) :: TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, scale, max_conv_vel, L_start, TDC_alpha_C, TDC_alpha_S type(auto_diff_real_star_order1) :: & - r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL + mixing_length_alpha, Lambda, r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, gradL type(auto_diff_real_star_order1) :: gradT, Y_face, conv_vel, D, Eq_div_w, grav, energy integer :: mixing_type, ierr, tdc_num_iters logical :: report, include_mlt_corr_to_TDC, use_TDC_enthalpy_flux_limiter @@ -202,6 +204,7 @@ subroutine check_TDC() Cp = 6628075118.4606590d0 opacity = 9.0750171231469945d-2 scale_height = 2638686602.0063782d0 + Lambda = mixing_length_alpha*scale_height gradL = 0.25207587267343501d0 grada = 0.25204697256872738d0 report = .false. @@ -221,7 +224,7 @@ subroutine check_TDC() do j = 0, 30 dt = 500d0*pow(1.02d0, j) call set_TDC( & - conv_vel_start, mixing_length_alpha, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, dt, cgrav, m, report, & + conv_vel_start, mixing_length_alpha, Lambda, TDC_alpha_D, TDC_alpha_R, TDC_alpha_Pt, 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, max_conv_vel, & Eq_div_w, grav, include_mlt_corr_to_TDC, TDC_alpha_C, TDC_alpha_S, use_TDC_enthalpy_flux_limiter, &