remove/modify fe_core_mass limit for tracking velocity

docs/source/changelog.rst

@@ -53,6 +53,10 @@ By user request, and motivated by the underestimation of line opacities from exp
 see also section 2.2 in `Morozova et al. (2015) <https://ui.adsabs.harvard.edu/abs/2015ApJ...814...63M/abstract>`_. An optional control for an opacity floor, ``opacity_min``, has
 been introduced.
 
+The Fe core-collapse infall condition ``fe_core_infall_limit`` has been adjusted, and the old control remains an optional infall condition.
+Users can switch between either choice with the new logical control ``report_max_infall_inside_fe_core``.
+See the ``&controls`` for further details.
+
 .. _Bug Fixes main:
 
 Bug Fixes

star/defaults/controls.defaults

@@ -1780,6 +1780,15 @@
     center_R_lower_limit = -1
 
 
+      ! report_max_infall_inside_fe_core
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! use ``fe_core_infall_mass`` to check amount of mass that is collapsing.
+      ! Then, If ``.true.``, fe_core_infall is determined by finding the maximum velocity that is inside
+      ! ``fe_core_mass`` in units of Msun. if ``.false.``, report the maximum infall velocity at any location.
+
+    report_max_infall_inside_fe_core = .true.
+
       ! fe_core_infall_limit
       ! ~~~~~~~~~~~~~~~~~~~~
 
@@ -1792,11 +1801,11 @@
       ! fe_core_infall_mass
       ! ~~~~~~~~~~~~~~~~~~~
 
-      ! Amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to ``fe_core_infall`` but there will be a greater chance of a
+      ! if ``check_mass_sum_for_infall = .true.``, ``fe_core_infall_mass`` is the amount of mass to check if collapsing, the smaller this is the closer the velocity minima will be to ``fe_core_infall`` but there will be a greater chance of a
       ! transistent velocity spike causing the model to prematurely exit.  In solar masses
 
 
-      fe_core_infall_mass = 0.1d0
+    fe_core_infall_mass = 0.1d0
 
 
       ! non_fe_core_infall_limit

star/private/ctrls_io.f90

@@ -355,8 +355,8 @@ module ctrls_io
     include_P_in_velocity_time_centering, include_L_in_velocity_time_centering, &
     P_theta_for_velocity_time_centering, L_theta_for_velocity_time_centering, &
     steps_before_use_TDC, use_P_d_1_div_rho_form_of_work_when_time_centering_velocity, compare_TDC_to_MLT, &
-    velocity_logT_lower_bound, max_dt_yrs_for_velocity_logT_lower_bound, velocity_tau_lower_bound, velocity_q_upper_bound, &
-    use_drag_energy, drag_coefficient, min_q_for_drag, &
+    velocity_logT_lower_bound, max_dt_yrs_for_velocity_logT_lower_bound, velocity_tau_lower_bound, velocity_q_upper_bound,  &
+    report_max_infall_inside_fe_core, use_drag_energy, drag_coefficient, min_q_for_drag, &
     v_drag_factor, v_drag, q_for_v_drag_full_off, q_for_v_drag_full_on, &
     retry_for_v_above_clight, &
 
@@ -1889,7 +1889,7 @@ subroutine store_controls(s, ierr)
  s% max_dt_yrs_for_velocity_logT_lower_bound = max_dt_yrs_for_velocity_logT_lower_bound
  s% velocity_tau_lower_bound = velocity_tau_lower_bound
  s% velocity_q_upper_bound = velocity_q_upper_bound
-
+ s% report_max_infall_inside_fe_core = report_max_infall_inside_fe_core
  s% retry_for_v_above_clight = retry_for_v_above_clight
 
  ! solvers
@@ -3571,7 +3571,7 @@ subroutine set_controls_for_writing(s, ierr)
  max_dt_yrs_for_velocity_logT_lower_bound = s% max_dt_yrs_for_velocity_logT_lower_bound
  velocity_tau_lower_bound = s% velocity_tau_lower_bound
  velocity_q_upper_bound = s% velocity_q_upper_bound
-
+ report_max_infall_inside_fe_core = s% report_max_infall_inside_fe_core
  retry_for_v_above_clight = s% retry_for_v_above_clight
 
  ! solvers

star/private/report.f90

@@ -150,7 +150,7 @@ subroutine do_report(s, ierr)
          type (star_info), pointer :: s
          integer, intent(out) :: ierr
 
-         integer :: k, nz, h1, h2, he3, he4, c12, n14, o16, ne20, si28, co56, ni56, k_min
+         integer :: k, nz, h1, h2, he3, he4, c12, n14, o16, ne20, si28, co56, ni56, k_min, k_fe
          real(dp) :: radius, dr, non_fe_core_mass, nu_for_delta_Pg, v, mstar, luminosity, mass_sum
          integer, pointer :: net_iso(:)
          real(dp), pointer :: velocity(:) => null()
@@ -381,10 +381,12 @@ subroutine do_report(s, ierr)
          if (failed('get_burn_zone_info')) return
 
 
-         s% fe_core_infall = 0
-         s% non_fe_core_infall = 0
-         s% non_fe_core_rebound = 0
-         s% max_infall_speed_mass = 0
+         s% fe_core_infall = 0d0
+         s% non_fe_core_infall = 0d0
+         s% non_fe_core_rebound = 0d0
+         s% max_infall_speed_mass = 0d0
+         k_fe = -1 ! initialize, in case u_flag & v_flag = .false.
+         k_min = -1 ! initialize, in case u_flag & v_flag = .false.
 
          if(s% u_flag .or. s% v_flag) then
 
@@ -399,14 +401,26 @@ subroutine do_report(s, ierr)
 
             mass_sum = 0d0
             if (s% fe_core_mass > 0) then
-               do k=1, nz
-                  if (s% m(k) > Msun*s% fe_core_mass) cycle
-                  if(-velocity(k) > s% fe_core_infall) mass_sum = mass_sum + s% dm(k)
+                ! check if [> fe_core_infall_mass] of core is infalling
+                ! instead of a for loop, we move inside out and only check regions inside the fe_core
+                k = nz
+                do while (k > 1)
+                   if (s% m(k) > Msun * s% fe_core_mass) exit  ! exit when outside Fe core
+                   if (-velocity(k) > s% fe_core_infall) then
+                      mass_sum = mass_sum + s% dm(k)
+                   end if
+                   k = k-1 ! loop outwards
                end do
+               k_fe = k ! mark fe_core_mass boundary cell location, k_fe = nz if no fe_core
 
-               if ((mass_sum > s% fe_core_infall_mass*msun) .and. &
-                   (s%m(k_min) <= s%fe_core_mass*msun)) then
-                  s% fe_core_infall = -velocity(k_min)
+               if (s% report_max_infall_inside_fe_core) then  ! report peak infall velocity inside fe_core (not necessarily the maximum, since infall tends to start outside in)
+                  if (mass_sum > s% fe_core_infall_mass*msun) then ! prevents toggling when k == nz
+                     s% fe_core_infall = - minval(s%v(k_fe:nz))
+                  end if
+               else !(default in r24.03.1 prior) report max infall velocity anywhere
+                  if(mass_sum > s% fe_core_infall_mass*msun) then
+                     s% fe_core_infall = -velocity(k_min)
+                  end if
                end if
             end if
 

star_data/private/star_controls.inc

@@ -836,7 +836,7 @@
             L_theta_for_velocity_time_centering
          integer :: steps_before_use_velocity_time_centering
 
-         logical :: use_drag_energy
+         logical :: use_drag_energy, report_max_infall_inside_fe_core
          real(dp) :: &
             drag_coefficient, &
             min_q_for_drag, &