From 126df74eef189748821e677c8197a7ef39d8a584 Mon Sep 17 00:00:00 2001 From: Debraheem Date: Mon, 1 Jun 2026 22:45:43 -0400 Subject: [PATCH 1/2] compute dynamical time and limit drag when dt is small to prevent blow up --- star/defaults/controls.defaults | 5 ++++- star/private/alloc.f90 | 2 ++ star/private/hydro_energy.f90 | 14 +++++++++----- star/private/hydro_momentum.f90 | 9 +++++++-- star/private/star_utils.f90 | 1 + star/private/struct_burn_mix.f90 | 10 ++++++++-- star_data/public/star_data_step_work.inc | 1 + 7 files changed, 32 insertions(+), 10 deletions(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 3f15424e3..e0296c9db 100644 --- a/star/defaults/controls.defaults +++ b/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. diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index 675345b21..735805865 100644 --- a/star/private/alloc.f90 +++ b/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) diff --git a/star/private/hydro_energy.f90 b/star/private/hydro_energy.f90 index b22d3c801..6690b341d 100644 --- a/star/private/hydro_energy.f90 +++ b/star/private/hydro_energy.f90 @@ -217,7 +217,7 @@ 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 - real(dp) :: alfa, beta + real(dp) :: alfa, beta, drag_rate 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, & @@ -283,18 +283,22 @@ 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 + drag_rate = s% drag_coefficient/max(s% dt, s% dynamic_timescale_start(k)) v_00 = wrap_v_00(s,k) - drag_force = s% drag_coefficient*v_00/s% dt + drag_force = drag_rate*v_00 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 + drag_rate = s% drag_coefficient/max(s% dt, s% dynamic_timescale_start(k+1)) v_p1 = wrap_v_p1(s,k) - drag_force = s% drag_coefficient*v_p1/s% dt + drag_force = drag_rate*v_p1 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 diff --git a/star/private/hydro_momentum.f90 b/star/private/hydro_momentum.f90 index ef8e620a3..23fe4a509 100644 --- a/star/private/hydro_momentum.f90 +++ b/star/private/hydro_momentum.f90 @@ -382,7 +382,7 @@ subroutine expected_non_HSE_term( & integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: extra_ad, v_00, & drag - real(dp) :: accel, d_accel_dv + real(dp) :: accel, d_accel_dv, drag_rate logical :: test_partials, local_v_flag include 'formats' @@ -423,7 +423,12 @@ subroutine expected_non_HSE_term( & if (s% q(k) > s% min_q_for_drag .and. s% drag_coefficient > 0) then v_00 = wrap_v_00(s,k) - drag = -s% drag_coefficient*v_00/s% dt + if (s% dynamic_timescale_start(k) > 0d0) then + drag_rate = s% drag_coefficient/max(s% dt, s% dynamic_timescale_start(k)) + else + drag_rate = 0d0 + end if + drag = -drag_rate*v_00 s% dvdt_drag(k) = drag%val end if diff --git a/star/private/star_utils.f90 b/star/private/star_utils.f90 index 051678cdf..0fda0a386 100644 --- a/star/private/star_utils.f90 +++ b/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 diff --git a/star/private/struct_burn_mix.f90 b/star/private/struct_burn_mix.f90 index 584f758ac..8dddabef8 100644 --- a/star/private/struct_burn_mix.f90 +++ b/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 - diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index 1991b965d..0c279912a 100644 --- a/star_data/public/star_data_step_work.inc +++ b/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) From f4268141ddb00f81386ca9b9bac1412f4b5cd214 Mon Sep 17 00:00:00 2001 From: Debraheem Date: Mon, 1 Jun 2026 23:10:39 -0400 Subject: [PATCH 2/2] clean up --- star/private/hydro_energy.f90 | 8 +++----- star/private/hydro_momentum.f90 | 13 +++++-------- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/star/private/hydro_energy.f90 b/star/private/hydro_energy.f90 index 6690b341d..2c1815048 100644 --- a/star/private/hydro_energy.f90 +++ b/star/private/hydro_energy.f90 @@ -217,7 +217,7 @@ 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 - real(dp) :: alfa, beta, drag_rate + 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, & @@ -285,9 +285,8 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad (s% drag_coefficient > 0) .and. & s% use_drag_energy .and. & s% dynamic_timescale_start(k) > 0d0) then - drag_rate = s% drag_coefficient/max(s% dt, s% dynamic_timescale_start(k)) v_00 = wrap_v_00(s,k) - drag_force = drag_rate*v_00 + 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 @@ -296,9 +295,8 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad (s% drag_coefficient > 0) .and. & s% use_drag_energy .and. & s% dynamic_timescale_start(k+1) > 0d0) then - drag_rate = s% drag_coefficient/max(s% dt, s% dynamic_timescale_start(k+1)) v_p1 = wrap_v_p1(s,k) - drag_force = drag_rate*v_p1 + 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 diff --git a/star/private/hydro_momentum.f90 b/star/private/hydro_momentum.f90 index 23fe4a509..cc913c752 100644 --- a/star/private/hydro_momentum.f90 +++ b/star/private/hydro_momentum.f90 @@ -382,7 +382,7 @@ subroutine expected_non_HSE_term( & integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: extra_ad, v_00, & drag - real(dp) :: accel, d_accel_dv, drag_rate + real(dp) :: accel, d_accel_dv logical :: test_partials, local_v_flag include 'formats' @@ -421,14 +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) - if (s% dynamic_timescale_start(k) > 0d0) then - drag_rate = s% drag_coefficient/max(s% dt, s% dynamic_timescale_start(k)) - else - drag_rate = 0d0 - end if - drag = -drag_rate*v_00 + drag = -s% drag_coefficient*v_00/max(s% dt, s% dynamic_timescale_start(k)) s% dvdt_drag(k) = drag%val end if