add tdc uq into HLLC

docs/source/changelog.rst

@@ -26,6 +26,13 @@ Removed ``file_extension`` option because it is redundant with ``file_device``.
 
 Renamed pgstar ``pause``` option to ``pause_flag``` because pause is a reserved Fortran 77 keyword.
 
+For greater consistency and clarity between ``TDC`` and ``RSP``, the controls for ``mlt_option = TDC`` have been renamed to match their respective ``RSP``
+counter part. See below :   
+
+- ``alpha_TDC_DAMP`` has been renamed to ``TDC_alpha_D`` and is analogous to ``RSP_alfad``
+- ``alpha_TDC_DAMPR`` has been renamed to ``TDC_alpha_R`` and is analogous to ``RSP_gammar``
+- ``alpha_TDC_PtdVdt`` has been renamed to ``TDC_alpha_Pt`` and is analogous to ``RSP_alfap``
+
 .. _New Features main:
 
 New Features

star/defaults/controls.defaults

@@ -2284,9 +2284,9 @@
       ! mlt_Pturb_factor
       ! ~~~~~~~~~~~~~~~~
 
-      ! include MLT turbulent pressure at face k = mlt_Pturb_factor*(1/3d0)*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3
-      ! MLT turbulent pressure for cell k = avg of values at faces. For TDC, mlt_Pturb_factor = 3d0
-      ! corresponds to alpha_TDC_Ptdvdt = 1d0, where Pt = (2/3)*rho(k)*(mlt_vc(k)/sqrt(2/3))**2
+      ! For MESA's momentume equation, include MLT turbulent pressure at face k = mlt_Pturb_factor*(1/3d0)*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3
+      ! MLT turbulent pressure for cell k = avg of values at faces. For TDC, mlt_Pturb_factor = 2d0
+      ! corresponds to TDC_alpha_Pt = 1d0, where Pt = (2/3)*rho(k)*(mlt_vc(k)/sqrt(2/3))**2
 
       ! this replaces conv_dP_term_factor.   also see extra_pressure vector and other_pressure routine
 
@@ -2328,16 +2328,16 @@
 
       ! TDC options
       ! ~~~~~~~~~~~
-      ! + ``alpha_TDC_DAMP`` : The turbulent viscous damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds.
-      ! + ``alpha_TDC_DAMPR`` : The radiative damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds.
-      ! + ``alpha_TDC_PtdVdt`` : The prefactor on the term accounting for work done against turbulent pressure (P_turb * dV/dt). Physically this should be unity.
+      ! + ``TDC_alpha_D`` : The turbulent viscous damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds. This control is analogous to ``RSP_alfad``, and previously called ``alpha_TDC_DAMP``
+      ! + ``TDC_alpha_R`` : The radiative damping parameter which determines the saturation of TDC. Increasing this decreases convection speeds. This control is analogous to ``RSP_gammar``, and previously called ``alpha_TDC_DAMPR``
+      ! + ``TDC_alpha_Pt`` : The prefactor on the term accounting for work done against turbulent pressure (P_turb * dV/dt). Physically this should be unity. This control is analogous to ``RSP_alfap``, and previously called ``alpha_TDC_PtdVdt``
       ! + ``steps_before_use_TDC`` : TDC often struggles with models on the pre-main-sequence. Set this option to pick MLT_option='Cox' for the first several steps to get past the pre-MS. Note that if this option is positive then only either TDC or Cox will be used (depending on model number). THIS OVERRIDES MLT_option!
 
       ! ::
 
-      alpha_TDC_DAMP = 1d0
-      alpha_TDC_DAMPR = 0d0
-      alpha_TDC_PtdVdt = 0d0
+      TDC_alpha_D = 1d0
+      TDC_alpha_R = 0d0
+      TDC_alpha_Pt = 0d0
       steps_before_use_TDC = 0
 
 

star/defaults/controls_dev.defaults

@@ -180,40 +180,48 @@
 ! ===
 
 
-      ! alpha_TDC_DAMPM
-      ! ~~~~~~~~~~~~~~~
-      ! alpha_TDC_C
+      ! TDC_alpha_M
+      ! ~~~~~~~~~~~
+      ! TDC_alpha_C
       ! ~~~~~~~~~~~
-      ! alpha_TDC_S
+      ! TDC_alpha_S
       ! ~~~~~~~~~~~
       ! TDC_use_density_form_for_eddy_viscosity
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-
-      ! If ``alpha_TDC_DAMPM>0``, then include eddy viscous damping in TDC ``alpha_TDC_DAMPM``
+      ! If ``TDC_alpha_M>0``, then include eddy viscous damping in TDC ``TDC_alpha_M``
       ! This control is analogous to ``RSP_alfam``, where the default is ``RSP_alfam = 0.25d0``.
 
       ! If hydrostatic (``v_flag``, ``u_flag = .false.``, ``v = 0`` ) there are no velocity gradients,
       ! and thus no shear to drive turbulence. Without shear, the eddy viscosity term becomes zero.
 
-      ! ``alpha_TDC_C`` and ``alpha_TDC_S`` are pre-factors to scale the flux and source terms.
+      ! ``TDC_alpha_C`` and ``TDC_alpha_S`` are pre-factors to scale the flux and source terms.
 
       ! If ``u_flag = .true.`` or ``TDC_use_density_form_for_eddy_viscosity = .true.``, use density
       ! derivative from newton solver to form d(v/r)/dr, used to compute Eq and Uq.
 
+      ! if ``TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = .true.``, ``mlt_vc_old`` (the start of step value)
+      ! is used in place of the current solver iterate ``mlt_vc_ad`` inside the hydro momentum equation. This option can
+      ! provide greater numerical stability in extremely dynamic situations, with very little loss in accuracy.
+
       ! ::
 
-    alpha_TDC_DAMPM = 0d0
-    alpha_TDC_C = 1.0d0
-    alpha_TDC_S = 1.0d0
+    TDC_alpha_M = 0d0
+    TDC_alpha_C = 1.0d0
+    TDC_alpha_S = 1.0d0
     TDC_use_density_form_for_eddy_viscosity = .false.
-
+    TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = .false.
 
       ! TDC_include_eturb_in_energy_equation
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      ! calculates and incorporates det/dt from TDC as a negative source term in
-      ! the total energy equation, along with the eddy viscous heating Eq
+      ! While ``TDC`` is already implicitly coupled to the total energy equation, this 
+      ! controls feedsback the turbulent energy into the total energy equation by
+      ! calculating and incorporating det/dt from TDC as a negative source term in
+      ! the total energy equation, along with the eddy viscous heating Eq. This comes
+      ! at the cost of numerical stability on long timescales, be warned. 
 
       ! ::
 
@@ -247,14 +255,18 @@
 
       ! TDC_num_innermost_cells_forced_nonturbulent
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! TDC_num_outermost_cells_forced_nonturbulent
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      ! Forces innermost TDC_num_innermost_cells_forced_nonturbulent cells
+      ! Forces innermost TDC_num_innermost_cells_forced_nonturbulent or
+      ! TDC_num_outermost_cells_forced_nonturbulent cells
       ! to be nonturbulent, and sets gradT = gradr in these cells.
       ! Useful for pulsation models in TDC.
 
       ! ::
 
     TDC_num_innermost_cells_forced_nonturbulent = 0
+    TDC_num_outermost_cells_forced_nonturbulent = 0
 
 
       ! include_mlt_corr_to_TDC
@@ -268,6 +280,16 @@
 
     include_mlt_corr_to_TDC = .true.
 
+      ! use_TDC_enthalpy_flux_limiter
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! If use_TDC_enthalpy_flux_limiter = .true. , apply enthalpy flux limiter to TDC 
+      ! similar to the form presented in Wuchterl & Feuchtinger 1998, and Smolec 2008.
+
+      ! ::
+
+    use_TDC_enthalpy_flux_limiter = .false.
+
 
       ! include_mlt_Pturb_in_thermodynamic_gradients
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -293,6 +315,9 @@
     include_mlt_in_velocity_time_centering = .false.
 
 
+! Remeshing 
+! =========
+
       ! remesh_for_TDC_pulsations_log_core_zoning
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -302,13 +327,22 @@
 
       ! ::
 
-      ! Remeshing parameters
-      ! ~~~~~~~~~~~~~~~~~~~~
-
     remesh_for_TDC_pulsations_log_core_zoning = .false.
     TDC_hydro_use_mass_interp_face_values = .false.
 
     TDC_hydro_nz = 150
     TDC_hydro_nz_outer = 40
     TDC_hydro_T_anchor = 11d3
     TDC_hydro_dq_1_factor = 2d0
+
+
+      ! use_hydro_merge_limits_in_mesh_plan
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Allows use of ``merge_amr_max_abs_du_div_cs``, ``merge_amr_du_div_cs_limit_only_for_compression``,
+      ! and ``merge_amr_inhibit_at_jumps`` in default mesh_plan when ``use_split_merge_amr = .false.``.
+      ! Only supports ``v_flag``.
+
+      ! ::
+
+    use_hydro_merge_limits_in_mesh_plan = .false.

star/defaults/profile_columns.list

@@ -420,6 +420,7 @@
    !cell_internal_energy_fraction
    !cell_internal_energy_fraction_start
    !cell_specific_PE
+   !dwork_dm ! cell specific work per unit time, Work = dwork_dm*dm*dt
 
    !log_cell_ie_div_star_ie
    !log_cell_specific_IE

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_common

@@ -124,11 +124,14 @@ use_Type2_opacities = .true.
 
    ! convection controls
     MLT_option = 'TDC'
-    alpha_TDC_DAMP = 1d0
-    alpha_TDC_DAMPR = 0d0
-    alpha_TDC_PtdVdt = 0d0
-    alpha_TDC_DAMPM = 0.25d0
+    TDC_alpha_D = 1d0
+    TDC_alpha_R = 0d0
+    TDC_alpha_Pt = 0d0
+    TDC_alpha_M = 0.25d0
+    TDC_alpha_M_use_explicit_mlt_vc_in_momentum_equation = .false.
     steps_before_use_TDC = 0
+    include_mlt_corr_to_TDC = .false. 
+    TDC_include_eturb_in_energy_equation = .false. ! unstable for evol
 
     use_ledoux_criterion = .true.
     alpha_semiconvection = 1d0

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pgstar

@@ -3,7 +3,7 @@
 
 !pause = .true.
 
-pgstar_interval = 50 ! making this too small slows the model down.
+pgstar_interval = 200 ! making this too small slows the model down.
 pgstar_show_age_in_years = .true.
 pgstar_show_age_in_seconds = .false.
 pgstar_sleep = 0.0
@@ -181,9 +181,9 @@ Profile_Panels1_title = ''
 Profile_Panels1_txt_scale = 0.7
 Profile_Panels1_num_panels = 4
 
-Profile_Panels1_yaxis_name(1) = 'Lc_div_L'
+Profile_Panels1_yaxis_name(1) = 'Eq'
 !Profile_Panels1_dymin(1) = 0.14
-Profile_Panels1_other_yaxis_name(1) = 'Frad_div_cUrad'!'lum_conv'
+Profile_Panels1_other_yaxis_name(1) = 'Uq'!'lum_conv'
 !Profile_Panels1_other_dymin(1) = 0.14
 
 Profile_Panels1_yaxis_name(2) = 'logT'!'logdq'!'radius'!'entropy'

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pulses

@@ -52,7 +52,7 @@ cubic_interpolation_in_Z = .false.
 &controls
 
    ! limit max model number as part of test suite.
-   max_model_number = 3000
+   !max_model_number = 3000
 
    ! controls for analyzing pulsations using GYRE in MESA
    x_integer_ctrl(1) = 1000 ! gyre interval to check
@@ -66,7 +66,7 @@ cubic_interpolation_in_Z = .false.
 
 
    ! new rsp style meshing, thanks to Bill P.
-   x_logical_ctrl(23) = .true. ! .true. = Remesh for TDC 
+   x_logical_ctrl(23) = .true. ! .true. = Remesh for TDC
    TDC_hydro_use_mass_interp_face_values = .true.
    TDC_hydro_nz = 150 
    TDC_hydro_nz_outer = 40 
@@ -75,9 +75,13 @@ cubic_interpolation_in_Z = .false.
    remesh_for_TDC_pulsations_log_core_zoning = .false. ! .false means do rsp style core
 
    okay_to_remesh = .false. ! freeze mesh after initial remesh.
-   x_logical_ctrl(24) = .true. ! if true turn off remesh at the following model number
+   x_logical_ctrl(24) = .false. ! if true turn off remesh at the following model number
    x_ctrl(12) = 200 ! model number to turn off remesh ( only if if okay_to_remesh = .true.)
 
+   !max_num_merge_surface_cells = 5
+   !max_surface_cell_dq= 1d-4
+   !min_surface_cell_dq = 1d-6
+
 ! GYRE set starting velocities, kick!
 ! kick when true and not restarting.
    x_logical_ctrl(5) = .true. ! to turn on gyre kick
@@ -151,16 +155,16 @@ cubic_interpolation_in_Z = .false.
    include_mlt_in_velocity_time_centering = .false.
    include_mlt_Pturb_in_thermodynamic_gradients = .false.
    mixing_length_alpha = 1.5d0
-   alpha_TDC_DAMP = 1d0
-   alpha_TDC_DAMPR = 0d0
-   alpha_TDC_PtdVdt = 0d0
-   alpha_TDC_DAMPM = 0.25d0
+   TDC_alpha_D = 1d0
+   TDC_alpha_R = 0d0
+   TDC_alpha_Pt = 0d0
+   TDC_alpha_M = 0.25d0
+   use_TDC_enthalpy_flux_limiter = .false.
    steps_before_use_TDC = 0
    mlt_Pturb_factor = 0d0 ! use 2d0 = 1d0 in TDC. time lagged.
    alt_scale_height_flag = .false. ! ignore eggleton in the core.
    TDC_num_innermost_cells_forced_nonturbulent = 2 ! for envelope models only.
 
-
    TDC_use_density_form_for_eddy_viscosity = .false. ! always true for u_flag.
    mlt_make_surface_no_mixing = .false.
 

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/profile_columns.list

@@ -418,6 +418,7 @@
    !cell_internal_energy_fraction
    !cell_internal_energy_fraction_start
    !cell_specific_PE
+   dwork_dm
 
    !log_cell_ie_div_star_ie
    !log_cell_specific_IE

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras.f90

@@ -490,6 +490,7 @@ integer function extras_start_step(id)
          ! if v= 0, turn on v so we can kick
          if (.not. s%v_flag .and. .not. s%u_flag) then
             call star_set_v_flag(id, .true., ierr)
+            write (*,*) 'turning on v_flag hydro for kick'
          end if
 
          call gyre_in_mesa_extras_set_velocities(s, .false., ierr)

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras_TDC_pulsation.inc

@@ -158,8 +158,8 @@
                v_surf = s% v(1)
                v_surf_start = s% v_start(1)
             else if (s% u_flag) then
-               v_surf = s% u_face_val(1)
-               v_surf_start = s% u_face_start(1)
+               v_surf = s% u(1) ! s% u_face_val(1)
+               v_surf_start = s% u_start(1) ! s% u_face_start(1)
             else !
                v_surf = 0d0
                v_surf_start = 0d0

star/other/other_mlt_results.f90

@@ -31,15 +31,15 @@ subroutine null_other_mlt_results(id, k, MLT_option, &  ! NOTE: k=0 is a valid a
                                      r, L, T, P, opacity, rho, 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, ierr)
+                                     mixing_type, gradT, Y_face, conv_vel, D, Gamma, energy, ierr)
       use const_def, only: dp
       use auto_diff
       use star_def
       integer, intent(in) :: id
       integer, intent(in) :: k
       character(len=*), intent(in) :: MLT_option
       type(auto_diff_real_star_order1), intent(in) :: &
-         r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height
+         r, L, T, P, opacity, rho, chiRho, chiT, Cp, gradr, grada, scale_height, energy
       integer, intent(in) :: iso
       real(dp), intent(in) :: &
          XH1, cgrav, m, gradL_composition_term, &

star/private/auto_diff_support.f90

@@ -1235,6 +1235,38 @@ function wrap_u_p1(s, k) result(v_p1)
          end if
       end function wrap_u_p1
 
+      function wrap_opt_time_center_u_m1(s, k) result(v_tc)
+         type (star_info), pointer :: s
+         type(auto_diff_real_star_order1) :: v_tc
+         integer, intent(in) :: k
+         v_tc = 0d0
+         if (k == 1) return
+         v_tc = wrap_u_m1(s,k)
+         if (s% using_velocity_time_centering) &
+            v_tc = 0.5d0*(v_tc + s% u_start(k-1))
+      end function wrap_opt_time_center_u_m1
+
+      function wrap_opt_time_center_u_00(s, k) result(v_tc)
+         type (star_info), pointer :: s
+         type(auto_diff_real_star_order1) :: v_tc
+         integer, intent(in) :: k
+         v_tc = wrap_u_00(s,k)
+         if (s% using_velocity_time_centering) &
+            v_tc = 0.5d0*(v_tc + s% u_start(k))
+      end function wrap_opt_time_center_u_00
+
+      function wrap_opt_time_center_u_p1(s, k) result(v_tc)
+         type (star_info), pointer :: s
+         type(auto_diff_real_star_order1) :: v_tc
+         integer, intent(in) :: k
+         v_tc = 0d0
+         if (k == s%nz) return
+         v_tc = wrap_u_p1(s,k)
+         if (s% using_velocity_time_centering) &
+            v_tc = 0.5d0*(v_tc + s% u_start(k+1))
+      end function wrap_opt_time_center_u_p1
+
+
       function wrap_Hp_m1(s, k) result(Hp_m1)
          type (star_info), pointer :: s
          type(auto_diff_real_star_order1) :: Hp_m1