adding eddy viscous damping to TDC

.gitignore

@@ -39,21 +39,13 @@ docs/source/_tags
 
 ## Output
 **/LOGS*/
-**/LOGS1/
-**/LOGS2/
 **/photos/
 **/png/
-**/png1/
-**/png2/
-**/png3/
+**/png*/
 **/pgstar_out/
-**/pgstar_out1/
-**/pgstar_out2/
-**/pgstar_out3/
+**/pgstar_out*/
 **/plt_out/
-**/plt_out1/
-**/plt_out2/
-**/plt_out3/
+**/plt_out*/
 
 ## Bill
 star_data/public_*

CODEOWNERS

@@ -21,7 +21,7 @@ interp_2d/ @rhdtownsend
 ionization/ @fxt44
 kap/ @fxt44 @Debraheem
 linters/
-turb/
+turb/ @Debraheem
 math/ @rhdtownsend
 mtx/ @mjoyceGR
 net/ @fxt44 @Debraheem

star/defaults/controls.defaults

@@ -2284,8 +2284,9 @@
       ! mlt_Pturb_factor
       ! ~~~~~~~~~~~~~~~~
 
-      ! include MLT turbulent pressure at face k = mlt_Pturb_factor*0.5*(rho(k) + rho(k-1))*mlt_vc(k)**2/3
-      ! MLT turbulent pressure for cell k = avg of values at faces.
+      ! 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
 
       ! this replaces conv_dP_term_factor.   also see extra_pressure vector and other_pressure routine
 
@@ -4594,6 +4595,20 @@
     Tsurf_factor = 1
 
 
+      ! use_RSP_L_eqn_outer_BC
+      ! ~~~~~~~~~~~~~~~~~~~~~~
+
+      ! use RSP2 luminosity equation
+      ! as the outer boundary condition on L rather than the
+      ! default atmospheric BC for temperature.
+      ! solves L1 - s% RSP2_Lsurf_factor*4*pi*r^2*c*a*T1^4 = 0
+      ! uses dev control s% RSP2_Lsurf_factor  = 0.5 (default),
+      ! see controls_dev.defaults.list
+
+      ! ::
+
+    use_RSP_L_eqn_outer_BC = .false.
+
       ! irradiation_flux
       ! ~~~~~~~~~~~~~~~~
       ! column_depth_for_irradiation
@@ -5880,11 +5895,15 @@
       ! max_surface_cell_dq
       ! ~~~~~~~~~~~~~~~~~~~
 
+      ! min_surface_cell_dq
+      ! ~~~~~~~~~~~~~~~~~~~
+
       ! Largest allowed dq at surface.
 
       ! ::
 
     max_surface_cell_dq = 1d-12
+    min_surface_cell_dq = 1d-99
 
 
       ! max_num_subcells
@@ -5899,12 +5918,15 @@
 
       ! max_num_merge_cells
       ! ~~~~~~~~~~~~~~~~~~~
+      ! max_num_merge_surface_cells
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
       ! Limits number of old cells to merge into 1 new one.
 
       ! ::
 
     max_num_merge_cells = 2
+    max_num_merge_surface_cells = 2
 
 
       ! mesh_adjust_get_T_from_E
@@ -6492,6 +6514,28 @@
     merge_amr_k_for_ignore_surface_cells = 2
 
 
+      ! merge_amr_ignore_core_cells
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! merge_amr_logT_for_ignore_core_cells
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! split_amr_ignore_core_cells
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! split_amr_logT_for_ignore_core_cells
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! For pulsation calculations, splitting and merging core cells can cause problems.
+      ! If merge_amr_ignore_core_cells or split_amr_ignore_surface_cells is true,
+      ! then the innermost cells above merge_amr_logT_for_ignore_core_cells or
+      ! split_amr_logT_for_ignore_core_cells  cells are ignored for split/merge.
+
+      ! ::
+
+    merge_amr_ignore_core_cells = .false.
+    merge_amr_logT_for_ignore_core_cells = 1d99
+    split_amr_ignore_core_cells = .false.
+    split_amr_logT_for_ignore_core_cells = 1d99
+
+
       ! split_merge_amr_avoid_repeated_remesh
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -7630,14 +7674,6 @@
 
     use_op_mono_alt_get_kap = .false.
 
-
-      ! min_kap_for_dPrad_dm_eqn
-      ! ~~~~~~~~~~~~~~~~~~~~~~~~
-
-      ! ::
-
-    min_kap_for_dPrad_dm_eqn = 1d-4
-
 ! asteroseismology controls
 ! =========================
 
@@ -8246,6 +8282,8 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       ! use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! min_kap_for_dPrad_dm_eqn
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~
 
 
       ! These are for alternatives ways to determine the T gradient.
@@ -8284,6 +8322,7 @@
       ! ::
 
     use_dPrad_dm_form_of_T_gradient_eqn = .false.
+    min_kap_for_dPrad_dm_eqn = 1d-4
     use_flux_limiting_with_dPrad_dm_form = .false.
     use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.
 

star/defaults/controls_dev.defaults

@@ -162,7 +162,7 @@
 
 
 ! split burn
-! ================
+! ==========
 
 
       ! op_split_burn_min_T_for_variable_T_solver
@@ -174,3 +174,140 @@
       ! ::
 
     op_split_burn_min_T_for_variable_T_solver = 1d99
+
+
+! TDC
+! ===
+
+
+      ! alpha_TDC_DAMPM
+      ! ~~~~~~~~~~~~~~~
+      ! alpha_TDC_C
+      ! ~~~~~~~~~~~
+      ! alpha_TDC_S
+      ! ~~~~~~~~~~~
+      ! TDC_use_density_form_for_eddy_viscosity
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+      ! If ``alpha_TDC_DAMPM>0``, then include eddy viscous damping in TDC ``alpha_TDC_DAMPM``
+      ! 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.
+
+      ! 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.
+
+      ! ::
+
+    alpha_TDC_DAMPM = 0.25d0
+    alpha_TDC_C = 1.0d0
+    alpha_TDC_S = 1.0d0
+    TDC_use_density_form_for_eddy_viscosity = .false.
+
+
+      ! TDC_include_eturb_in_energy_equation
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! calculates and incorporates det/dt from TDC as a negative source term in
+      ! the total energy equation
+
+      ! ::
+
+    TDC_include_eturb_in_energy_equation = .false.
+
+      ! use_rsp_form_of_scale_height
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! calculate scale height by averaging P/rho onto faces together instead
+      ! of of wrapping P and rho separately by default. Both forms assume HSE.
+
+      ! ::
+
+      use_rsp_form_of_scale_height = .false.
+
+      ! TDC_num_innermost_cells_forced_nonturbulent
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Forces innermost TDC_num_innermost_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
+
+
+      ! include_mlt_corr_to_TDC
+      ! ~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! If include_mlt_corr_to_TDC = .true. , mlt correction to TDC follows MESA VI (Jermyn et a. 2023)
+      ! in which we correct TDC on long timescales with Y_face ~ Y*Gamma/(1+Gamma). If .false., abandon
+      ! mlt limit of tdc and adopt original form of (Kuhfuß 1986) convection model in the local limit.
+
+      ! ::
+
+    include_mlt_corr_to_TDC = .true.
+
+
+      ! include_mlt_Pturb_in_thermodynamic_gradients
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! If mlt_Pturb_factor > 0, includes the effect of turbulent pressure Pturb
+      ! on thermodynamic gradients gradr, grada, and hence gradT. The correction
+      ! is approximated as correction = Prad + Pgas / Ptotal, where Ptotal = Prad + Pgas + Pturb.
+      ! Not supported by other_alpha_mlt or starspots
+
+      ! ::
+
+    include_mlt_Pturb_in_thermodynamic_gradients = .false.
+
+
+      ! include_mlt_in_velocity_time_centering
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! If doing velocity time_centering, include time centering in mlt. Depending on the time
+      ! centering flags adopted, include_mlt_in_velocity_time_centering includes time centering
+      ! for P, L, and r in geff, gradr, grada, and mlt.
+
+      ! ::
+
+    include_mlt_in_velocity_time_centering = .false.
+
+      ! make_mlt_hydrodynamic
+      ! ~~~~~~~~~~~~~~~~~~~~~
+
+
+      ! If make_mlt_hydrodynamic = .true., adopt geff = g - dv/dt from newton solver.
+      ! gradr and scale_height are modified by this flag to hydro form.
+      ! If .false. default to QHSE form, where g = GM/r^2
+      ! semiconvection, and use_mass_corrections ignore this flag, see star_utils.
+      ! not supported by RSP or RSP2.
+
+      ! ::
+
+    make_mlt_hydrodynamic = .false.
+
+
+      ! remesh_for_TDC_pulsations_log_core_zoning
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! If remesh_for_TDC_pulsations, .true. will rezone below RSP2_Tanchor with logarithmic
+      ! spacing in mass. If .false., remeshing scheme will do a bisection root find such
+      ! that each zone increases in mass inward following a power law.
+
+      ! ::
+
+      ! Remeshing parameters
+      ! ~~~~~~~~~~~~~~~~~~~~
+
+    remesh_for_TDC_pulsations_log_core_zoning = .false.
+    TDC_pulse_use_mass_interp_face_values = .false.
+
+    TDC_pulse_nz = 150
+    TDC_pulse_nz_outer = 40
+    TDC_pulse_T_anchor = 11d3
+    TDC_pulse_dq_1_factor = 2d0

star/defaults/profile_columns.list

@@ -883,7 +883,8 @@
    !log_du_kick_div_du
    !du
    !dvdt_dPdm
-   !dvdt_grav
+   !dvdt_grav ! gravity = -grav = -Gm/r^2
+   !grav_eff ! effective gravity = -(grav - dvdt)
 
    !tau_conv
    !tau_cool

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_15M/README.rst

@@ -0,0 +1,21 @@
+.. _ppisn:
+
+*****
+ppisn
+*****
+
+This test case evolves a very massive helium star from the He-ZAMS
+up to the ocurrence of a pulsational pair-instability event (see |Marchant2019|).
+
+.. |Marchant2019| replace:: `Marchant et al. 2019 <https://ui.adsabs.harvard.edu/abs/2019ApJ...882...36M/abstract>`__
+
+Initialization of the model
+===========================
+The initial mass of the helium star is set in ``inlist_extra``
+
+.. literalinclude:: ../../../star/test_suite/ppisn/inlist_extra
+
+In this case we use a :math:`72 M_\odot`
+
+Last-Updated: 2019-11-12 (mesa r12413) by Pablo Marchant
+

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_15M/ck

@@ -0,0 +1,7 @@
+#!/bin/bash
+
+# this provides the definition of check_one
+# check_one
+source "${MESA_DIR}/star/test_suite/test_suite_helpers"
+
+check_one

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_15M/clean

@@ -0,0 +1,4 @@
+#!/bin/bash
+
+cd make
+make clean

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_15M/gyre.in

@@ -0,0 +1,34 @@
+&model
+  add_center = .FALSE.
+/
+
+&mode
+   l = 0
+/
+
+&osc
+   inner_bound = 'ZERO_R' ! for envelopes
+   nonadiabatic = .TRUE.
+/
+
+&rot
+/
+
+&num
+   diff_scheme = 'MAGNUS_GL2'!'COLLOC_GL2' 
+/
+
+&scan
+   grid_type = 'LINEAR'
+   freq_min = 0.5!1d-4
+   freq_max = 5.0!10.0
+   freq_min_units = 'ACOUSTIC_DELTA'
+   freq_max_units = 'ACOUSTIC_DELTA'
+   n_freq = 50
+/
+
+&grid
+   w_osc = 10
+   w_exp = 2
+   w_ctr = 10
+/