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..2c1815048 100644 --- a/star/private/hydro_energy.f90 +++ b/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 diff --git a/star/private/hydro_momentum.f90 b/star/private/hydro_momentum.f90 index ef8e620a3..cc913c752 100644 --- a/star/private/hydro_momentum.f90 +++ b/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 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)