changes to tdc remesh and dev cases

docs/source/changelog.rst

@@ -39,6 +39,12 @@ GYRE has been upgraded to 9.1.1, the most recent stable release. Changes since t
 
 The `MESA SDK <http://user.astro.wisc.edu/~townsend/static.php?ref=mesasdk>`__ recommended for compiling MESA has been updated to 26.6.1. Although this newer SDK is not required to successfully build MESA, it brings the benefit of restoring the cross-platform bit-for-bit compatibility that MESA once enjoyed (meaning that runs on Linux/Intel, MacOS/Intel and MacOS/ARM give identical results).
 
+Additional controls are available for TDC envelope remeshing:
+
+- ``remesh_for_TDC_pulsations_when_load`` remeshes a model after loading it from a ``.mod`` file.
+- ``TDC_hydro_nz_inner`` adds geometrically spaced zones near the inner boundary.
+- ``TDC_hydro_nz_T_gradient`` adds zones according to the variation in ``logT`` below ``TDC_hydro_T_anchor`` while retaining the mass-based mesh.
+
 .. _Bug Fixes main:
 
 Bug Fixes

star/defaults/controls.defaults

@@ -6527,37 +6527,60 @@
       ! ~~~~~~~~~~~~
       ! TDC_hydro_nz_outer
       ! ~~~~~~~~~~~~~~~~~~
+      ! TDC_hydro_nz_inner
+      ! ~~~~~~~~~~~~~~~~~~
+      ! TDC_hydro_nz_T_gradient
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~
       ! TDC_hydro_T_anchor
       ! ~~~~~~~~~~~~~~~~~~
       ! TDC_hydro_dq_1_factor
       ! ~~~~~~~~~~~~~~~~~~~~~
       ! TDC_hydro_use_mass_interp_face_values
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! remesh_for_TDC_pulsations_when_load
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       ! remesh_for_TDC_pulsations_log_core_zoning
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      ! If calling the public ``remesh_for_TDC_pulsations`` function, the model will be rezoned below ``TDC_hydro_T_anchor`` 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, similar to RSP's static model builder.
+      ! If ``remesh_for_TDC_pulsations_log_core_zoning`` is true, the model is rezoned below
+      ! ``TDC_hydro_T_anchor`` with logarithmic spacing in mass. If false, each cell increases
+      ! in mass inward following a power law, similar to the RSP static model builder.
 
       ! The total number of zones is set by ``TDC_hydro_nz``, with constant spacing in mass between the surface and
-      ! ``TDC_hydro_T_anchor``, and increasing spacing in mass between the anchor and the inner boundary condition.
+      ! ``TDC_hydro_T_anchor``. The selected core zoning is used between the anchor and the inner boundary condition.
       ! The number of zones in between the anchor and the surface is set by ``TDC_hydro_nz_outer``, and the number of
       ! zones in the interior is set by ``TDC_hydro_nz`` - ``TDC_hydro_nz_outer``.
       ! The ``TDC_hydro_dq_1_factor`` sets the dq of the surface cell to this multiple of the cell beneath it, k = 2.
 
+      ! For power-law core zoning, ``TDC_hydro_nz_inner`` reserves this many zones next to the inner boundary.
+      ! Their cell masses decrease geometrically inward to ``max_center_cell_dq``. Set this to zero to use a
+      ! single inward-increasing power law. Increase ``TDC_hydro_nz`` by the same amount to add inner zones
+      ! without removing zones from the middle of an existing mesh. This control has no effect with logarithmic
+      ! core zoning.
+
+      ! ``TDC_hydro_nz_T_gradient`` adds this many zones below ``TDC_hydro_T_anchor`` according to the variation
+      ! in ``logT``. These zones enhance the selected mass-based core zoning rather than replacing it. Increase
+      ! ``TDC_hydro_nz`` by the same amount to retain all zones in the underlying mass grid. Set this to zero to
+      ! leave the mass grid unchanged.
+
       ! ``TDC_hydro_use_mass_interp_face_values`` determines whether quantities are averaged to faces using simple average or
       ! mass weighted averaging.
 
-      ! This remeshing scheme works well for most Cepheid envelopes, but it is not designed for AGB stars with deep convective envelopes.
+      ! Use ``TDC_hydro_nz_inner`` and ``TDC_hydro_nz_T_gradient`` to resolve the inner boundary
+      ! and steep temperature gradients in AGB envelopes with deep convective regions.
+
+      ! If ``remesh_for_TDC_pulsations_when_load`` is true, a model loaded from a ``.mod`` file
+      ! is remeshed before the first step. This is not applied to photo restarts.
 
       ! ::
 
+    remesh_for_TDC_pulsations_when_load = .false.
     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_nz_inner = 0
+    TDC_hydro_nz_T_gradient = 0
     TDC_hydro_T_anchor = 11d3
     TDC_hydro_dq_1_factor = 2d0
 

star/defaults/controls_dev.defaults

@@ -201,15 +201,6 @@
 ! ===
 
 
-      ! TDC_use_density_form_for_eddy_viscosity
-      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-      ! 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.
-
-    TDC_use_density_form_for_eddy_viscosity = .false.
-
-
       ! TDC_adjust_mass_fallback_to_mlt
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pgstar

@@ -193,7 +193,7 @@ Profile_Panels1_other_yaxis_log(2) = .false.
 !Profile_Panels1_ymax(2) = 210
 !Profile_Panels1_ymin(2) = 160
 
-Profile_Panels1_yaxis_name(3) = 'Lc_div_L'!'gradT'
+Profile_Panels1_yaxis_name(3) = 'logdq'!'gradT'
 Profile_Panels1_other_yaxis_name(3) = 'conv_vel'!'lum_div_Ledd'!'grada'
 Profile_Panels1_same_yaxis_range(3) = .false.
 Profile_Panels1_other_dymin(3) = 0.08
@@ -209,7 +209,7 @@ Profile_Panels1_other_dymin(4) = 0.14
 ! x-axis limits and properties
 Profile_Panels1_xaxis_name = 'logtau'!'zone'
 Profile_Panels1_xmin = -101d0!-101d0
-Profile_Panels1_xmax = 5.5d0!-101d0!9d0 !-101d0 !8.1
+Profile_Panels1_xmax = -101d0!9d0 !-101d0 !8.1
 Profile_Panels1_xaxis_reversed = .true.
 
 !Profile_Panels1_xaxis_name = 'zone'

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/inlist_pulses

@@ -27,7 +27,7 @@ cubic_interpolation_in_Z = .false.
 
    relax_to_this_tau_factor = 1d-3
    dlogtau_factor = 0.1d0
-   relax_tau_factor = .true.
+   relax_tau_factor = .false.
 
    set_initial_dt = .true.
    years_for_initial_dt = -1
@@ -59,17 +59,19 @@ cubic_interpolation_in_Z = .false.
 
    x_logical_ctrl(7) = .true. ! doing pulses
    x_integer_ctrl(7) = -1 ! 3 ! which period to check (<= 0 means don't check any)
-   x_ctrl(7) = 12d0 ! expected period (in days) 
+   x_ctrl(7) = -1d0 ! expected period in days; <= 0 uses the dynamical time until a period is measured
    x_ctrl(8) = -1 ! min_deltaR_for_periods (Rsun)
    x_ctrl(9) = 1d0 ! KE_growth_avg_abs_frac_new ! for averaging growth rates across steps.
-   x_ctrl(10) = 0.3 ! min_period_div_target
+   x_ctrl(10) = 0.3d0 ! min_period_div_target
 
 
-   ! new rsp style meshing, thanks to Bill P.
-   x_logical_ctrl(23) = .true. ! .true. = Remesh for TDC
+   ! TDC pulsation mesh
+   remesh_for_TDC_pulsations_when_load = .true.
    TDC_hydro_use_mass_interp_face_values = .true.
-   TDC_hydro_nz = 150 
+   TDC_hydro_nz = 190
    TDC_hydro_nz_outer = 40 
+   TDC_hydro_nz_inner = 20
+   TDC_hydro_nz_T_gradient = 20
    TDC_hydro_T_anchor = 11d3
    TDC_hydro_dq_1_factor = 2d0
    remesh_for_TDC_pulsations_log_core_zoning = .false. ! .false means do rsp style core
@@ -79,8 +81,8 @@ cubic_interpolation_in_Z = .false.
    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
+   !max_surface_cell_dq= 1d0
+   !min_surface_cell_dq = 1d-8
 
 ! GYRE set starting velocities, kick!
 ! kick when true and not restarting.
@@ -101,10 +103,6 @@ cubic_interpolation_in_Z = .false.
    ! turn of wind during Pulsation
    use_other_wind = .false.
 
-   ! can optionally use QHSE form of dp/dm or dprad/dm form below.
-   use_dPrad_dm_form_of_T_gradient_eqn = .false. ! if false use QHSE form.
-   min_kap_for_dPrad_dm_eqn = 1d-5
-
 ! timesteps for saturation
 
   ! TDC Pulsation timestepping.
@@ -145,13 +143,12 @@ cubic_interpolation_in_Z = .false.
    RSP2_Lsurf_factor = 0.5d0
    use_zero_Pgas_outer_BC = .true.
    use_compression_outer_BC = .false.
-   use_momentum_outer_BC = .false. ! adopt this if not using rsp and zero P bc.
+   use_momentum_outer_BC = .true. ! adopt this if not using rsp and zero P bc.
 
 ! Convection model
    MLT_option = 'TDC'
    include_mlt_corr_to_TDC = .false. ! true K. 1986 model, no mlt limiting
    TDC_include_eturb_in_energy_equation = .true.
-   use_rsp_form_of_scale_height = .true.
    include_mlt_in_velocity_time_centering = .false.
    include_mlt_Pturb_in_thermodynamic_gradients = .false.
    mixing_length_alpha = 1.5d0
@@ -163,11 +160,16 @@ cubic_interpolation_in_Z = .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_num_innermost_cells_forced_nonturbulent = 0 ! for envelope models only.
 
-   TDC_use_density_form_for_eddy_viscosity = .false. ! always true for u_flag.
    mlt_make_surface_no_mixing = .false.
 
+
+   ! hydro solver
+   use_TDC_Y_face_seeded_newton = .true.
+   use_face_reconstruction = .true. ! accurate face reconstruction
+   hydro_matrix_solver = 'banded'
+
 ! controls for shock capturing, relaxed for pulsations
    ! main purpose is to force radiative in shock face
     max_abs_du_div_cs_for_convection = 1d0 !1d-1 !0.03d0

star/dev_cases_TDC_Pulsation/dev_TDC_Cepheid_6M/src/run_star_extras.f90

@@ -44,7 +44,7 @@ module run_star_extras
    ! values specified on inlist_common, inlist_pulses
       !!!!!!!!!!!!!!!!!!!!!!!!!
 
-   logical :: in_inlist_pulses, remesh_for_envelope_model, turn_off_remesh
+   logical :: in_inlist_pulses, turn_off_remesh
    integer :: kick_model_number, timestep_drop_model_number, turn_off_remesh_model_number
    integer :: initial_model_number
    real(dp) :: max_dt_before_pulse, max_dt_during_pulse
@@ -87,11 +87,10 @@ subroutine extras_controls(id, ierr)
       in_inlist_pulses = s%x_logical_ctrl(22)
       max_dt_before_pulse = s%x_ctrl(17)
       max_dt_during_pulse = s%x_ctrl(18)
-      remesh_for_envelope_model = s%x_logical_ctrl(23)
       turn_off_remesh = s%x_logical_ctrl(24)
-      kick_model_number = s%x_ctrl(11)
-      timestep_drop_model_number = s%x_ctrl(13)
-      turn_off_remesh_model_number = s%x_ctrl(12)
+      kick_model_number = int(s%x_ctrl(11))
+      timestep_drop_model_number = int(s%x_ctrl(13))
+      turn_off_remesh_model_number = int(s%x_ctrl(12))
    end subroutine extras_controls
 
    subroutine brott_wind(id, Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z, w, ierr)
@@ -245,7 +244,6 @@ subroutine my_adjust_mdot(id, ierr)
       integer, intent(in) :: id
       integer, intent(out) :: ierr
       type(star_info), pointer :: s
-      real(dp) :: Lrad_div_Ledd
       ierr = 0
       call star_ptr(id, s, ierr)
       if (ierr /= 0) return
@@ -299,9 +297,7 @@ subroutine my_other_kap_get( &
 
       type(star_info), pointer :: s
       real(dp) :: velocity
-      real(dp) :: radius, logR
-      real(dp) :: logT_alt, inv_diff
-      real(dp) :: log_kap, alpha
+      real(dp) :: radius
 
       ierr = 0
       call star_ptr(id, s, ierr)
@@ -378,18 +374,12 @@ subroutine extras_startup(id, restart, ierr)
       end if
       !initial_model_number = 0 ! since we are setting model # to 0 in inlist_pulses
 
-      ! for rsp style mesh
-      if (.not. restart .and. in_inlist_pulses .and. remesh_for_envelope_model) then
-         call remesh_for_TDC_pulsation(id, ierr)
-      end if
    end subroutine extras_startup
 
    subroutine extras_after_evolve(id, ierr)
       integer, intent(in) :: id
       integer, intent(out) :: ierr
       type(star_info), pointer :: s
-      real(dp) :: dt
-      character(len=strlen) :: test
       ierr = 0
       call star_ptr(id, s, ierr)
       if (ierr /= 0) return
@@ -400,10 +390,9 @@ subroutine extras_after_evolve(id, ierr)
    end subroutine extras_after_evolve
 
    ! returns either keep_going, retry, or terminate.
-   integer function extras_check_model(id)
+      integer function extras_check_model(id)
       integer, intent(in) :: id
-      integer :: ierr, k
-      real(dp) :: max_v
+      integer :: ierr
       type(star_info), pointer :: s
       include 'formats'
       ierr = 0
@@ -426,10 +415,9 @@ end function how_many_extra_history_columns
    subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
       integer, intent(in) :: id, n
       character(len=maxlen_history_column_name) :: names(n)
-      real(dp) :: vals(n), v_esc
+      real(dp) :: vals(n)
       integer, intent(out) :: ierr
       type(star_info), pointer :: s
-      integer :: k, k0
       ierr = 0
       call star_ptr(id, s, ierr)
       if (ierr /= 0) return
@@ -456,7 +444,6 @@ subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr)
       real(dp) :: vals(nz, n)
       integer, intent(out) :: ierr
       type(star_info), pointer :: s
-      integer :: k
       ierr = 0
       call star_ptr(id, s, ierr)
       if (ierr /= 0) return
@@ -517,11 +504,13 @@ end function extras_start_step
    subroutine my_before_struct_burn_mix(id, dt, res)
       use const_def, only: dp
       use star_def
+      use utils_lib, only: is_bad
       integer, intent(in) :: id
       real(dp), intent(in) :: dt
       integer, intent(out) :: res ! keep_going, redo, retry, terminate
-      real(dp) :: power_photo, v_esc
-      integer :: ierr, k
+      real(dp) :: dt_limit
+      integer :: ierr
+      logical :: have_dt_limit
       type(star_info), pointer :: s
       include 'formats'
       ierr = 0
@@ -530,15 +519,32 @@ subroutine my_before_struct_burn_mix(id, dt, res)
 
       if (in_inlist_pulses) then
          if (s%model_number > timestep_drop_model_number) then
-            s%max_timestep = max_dt_during_pulse
+            if (max_dt_during_pulse > 0d0) s%max_timestep = max_dt_during_pulse
          else
-            s%max_timestep = max_dt_before_pulse
+            if (max_dt_before_pulse > 0d0) s%max_timestep = max_dt_before_pulse
          end if
 
          ! time step control on pulsations
-         if (period > 0d0 .and. period/s%max_timestep < 600 .and. &
-             s%model_number > timestep_drop_model_number) then
-            s%max_timestep = period/600d0
+         have_dt_limit = .false.
+         dt_limit = 0d0
+         if (s%model_number > timestep_drop_model_number) then
+            if (num_periods < 1) then
+               if (.not. is_bad(s%dynamic_timescale) .and. s%dynamic_timescale > 0d0) then
+                  dt_limit = s%dynamic_timescale/600d0
+                  have_dt_limit = .true.
+               end if
+            else if (period > 0d0) then
+               dt_limit = period/600d0
+               have_dt_limit = .true.
+            end if
+         end if
+
+         if (have_dt_limit) then
+            if (s%max_timestep <= 0d0) then
+               s%max_timestep = dt_limit
+            else
+               s%max_timestep = min(s%max_timestep, dt_limit)
+            end if
          end if
 
          if (s%model_number > turn_off_remesh_model_number .and. turn_off_remesh) then
@@ -563,10 +569,10 @@ integer function extras_finish_step(id)
       use run_star_support
       use math_lib
       integer, intent(in) :: id
-      integer :: ierr, k
-      real(dp) :: max_vel_inside, vesc_for_cell, vesc_surf !check_avg_v_div_vesc
+      integer :: ierr
       type(star_info), pointer :: s
       include 'formats'
+      extras_finish_step = terminate
       ierr = 0
       call star_ptr(id, s, ierr)
       if (ierr /= 0) return
@@ -594,4 +600,3 @@ subroutine photo_read(id, iounit, ierr)
    end subroutine photo_read
 
 end module run_star_extras
-