limit drag when dt is small to prevent blow up

star/defaults/controls.defaults

@@ -8496,7 +8496,10 @@
 
       ! only when v_flag = .true.. Adjusts both v and energy transfer from kinetic to thermal.
       ! only for v(k) when q(k) > min_q_for_drag.
-      ! kill off fraction of v = drag_coefficient (i.e. set to 1 to keep v near 0)
+      ! drag_coefficient sets the drag rate relative to the local face
+      ! dynamical timescale at the start of the step. The effective implicit
+      ! per-step damping strength is
+      ! ``drag_coefficient*dt/max(dt,t_dyn_face)``.
       ! useful for preventing the development radial pulsations during advanced
       ! burning in massive stars and AGB stars.
 

star/private/alloc.f90

@@ -1220,6 +1220,8 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             if (failed('L_start')) exit
             call do1(s% r_start, c% r_start)
             if (failed('r_start')) exit
+            call do1(s% dynamic_timescale_start, c% dynamic_timescale_start)
+            if (failed('dynamic_timescale_start')) exit
             call do1(s% rmid_start, c% rmid_start)
             if (failed('rmid_start')) exit
             call do1(s% omega_start, c% omega_start)

star/private/hydro_energy.f90

@@ -283,18 +283,20 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad
             if (k /= s% nz) then
                if ((s% q(k) > s% min_q_for_drag) .and. &
                     (s% drag_coefficient > 0) .and. &
-                    s% use_drag_energy) then
+                    s% use_drag_energy .and. &
+                    s% dynamic_timescale_start(k) > 0d0) then
                   v_00 = wrap_v_00(s,k)
-                  drag_force = s% drag_coefficient*v_00/s% dt
+                  drag_force = s% drag_coefficient*v_00/max(s% dt, s% dynamic_timescale_start(k))
                   drag_energy = 0.5d0*v_00*drag_force
                   s% FdotV_drag_energy(k) = drag_energy%val
                ! drag energy for outer half-cell.   the 0.5d0 is for dm/2
                end if
                if ((s% q(k+1) > s% min_q_for_drag) .and. &
                     (s% drag_coefficient > 0) .and. &
-                    s% use_drag_energy) then
+                    s% use_drag_energy .and. &
+                    s% dynamic_timescale_start(k+1) > 0d0) then
                   v_p1 = wrap_v_p1(s,k)
-                  drag_force = s% drag_coefficient*v_p1/s% dt
+                  drag_force = s% drag_coefficient*v_p1/max(s% dt, s% dynamic_timescale_start(k+1))
                   drag_energy = drag_energy + 0.5d0*v_p1*drag_force
                   s% FdotV_drag_energy(k) = drag_energy%val
                ! drag energy for inner half-cell.   the 0.5d0 is for dm/2

star/private/hydro_momentum.f90

@@ -421,9 +421,11 @@ subroutine expected_non_HSE_term( &
             accel_ad%val = accel
             accel_ad%d1Array(i_v_00) = d_accel_dv
 
-            if (s% q(k) > s% min_q_for_drag .and. s% drag_coefficient > 0) then
+            if (s% q(k) > s% min_q_for_drag .and. &
+                  s% drag_coefficient > 0 .and. &
+                  s% dynamic_timescale_start(k) > 0d0) then
                v_00 = wrap_v_00(s,k)
-               drag = -s% drag_coefficient*v_00/s% dt
+               drag = -s% drag_coefficient*v_00/max(s% dt, s% dynamic_timescale_start(k))
                s% dvdt_drag(k) = drag%val
             end if
 

star/private/star_utils.f90

@@ -1567,6 +1567,7 @@ subroutine reset_starting_vectors(s)
          do k=1,s% nz
             s% T_start(k) = -1d99
             s% r_start(k) = -1d99
+            s% dynamic_timescale_start(k) = -1d99
             s% rmid_start(k) = -1d99
             s% v_start(k) = -1d99
             s% u_start(k) = -1d99

star/private/struct_burn_mix.f90

@@ -20,7 +20,8 @@
       module struct_burn_mix
 
       use star_private_def
-      use const_def, only: dp, i8, ln10, secyer, lsun
+      use const_def, only: dp, i8, ln10, pi, secyer, lsun
+      use math_lib, only: pow3
       use utils_lib, only: is_bad
 
       implicit none
@@ -328,6 +329,12 @@ subroutine save_start_values(s, ierr)
             s% eps_nuc_start(k) = s% eps_nuc(k)
             s% opacity_start(k) = s% opacity(k)
             s% m_grav_start(k) = s% m_grav(k)
+            if (s% r_start(k) > 0d0 .and. s% m_grav_start(k) > 0d0 .and. s% cgrav(k) > 0d0) then
+               s% dynamic_timescale_start(k) = &
+                  2d0*pi*sqrt(pow3(s% r_start(k))/(s% cgrav(k)*s% m_grav_start(k)))
+            else
+               s% dynamic_timescale_start(k) = 0d0
+            end if
          end do
 
          if (s% RSP2_flag) then
@@ -1211,4 +1218,3 @@ end subroutine burn1_zone
 
       end module struct_burn_mix
 
-

star_data/public/star_data_step_work.inc

@@ -828,6 +828,7 @@
          real(dp), pointer :: u_start(:) ! (nz)
          real(dp), pointer :: L_start(:) ! (nz)
          real(dp), pointer :: r_start(:) ! (nz)
+         real(dp), pointer :: dynamic_timescale_start(:) ! local face dynamical timescale at start of step (nz)
          real(dp), pointer :: rmid_start(:) ! (nz)
          real(dp), pointer :: rho_start(:) ! (nz)
          real(dp), pointer :: omega_start(:) ! (nz)