rates/net/solver robustness tweaks/fixes

net/private/net_eval.f90

@@ -607,7 +607,7 @@ subroutine get_rates_with_screening(n, ierr)
          if (dbg) write(*,*) 'call eval_using_rate_tables'
          call eval_using_rate_tables( &
             g% num_reactions, g% reaction_id, g% rate_table, g% rattab_f1, nrattab,  &
-            n% ye, n% logT, n% temp, n% rho, n% rate_factors, g% logttab, &
+            n% ye, n% logT, n% temp, n% rho, n% rate_factors, g% ttab, g% logttab, &
             n% rate_raw, n% rate_raw_dT, n% rate_raw_dRho, ierr)
          if (ierr /= 0) then
             if (dbg) write(*,*) 'ierr from eval_using_rate_tables'

net/test/src/test_net_do_one.f90

@@ -136,6 +136,10 @@ subroutine do1_net(handle, symbolic)
          call mesa_error(__FILE__, __LINE__)
       end if
 
+      if (net_file == 'approx21.net') then
+         call test_approx21_reduced_flow_derivatives(n)
+      end if
+
       if (symbolic .and. .not. qt) then
          write (*, *) 'nonzero d_dxdt_dx entries'
          k = 0
@@ -180,6 +184,99 @@ subroutine do1_net(handle, symbolic)
 
    end subroutine do1_net
 
+   subroutine test_approx21_reduced_flow_derivatives(n)
+      use net_approx21, only: eval_approx21_reduced_flows
+      type(Net_Info), intent(in) :: n
+
+      integer :: i, j, i_max, j_max
+      real(qp), parameter :: temp = 2.0e9_qp, rel_step = 1.0e-5_qp, &
+         tolerance = 1.0e-7_qp
+      real(qp) :: deriv, error, max_error, scale, step
+      real(qp), allocatable :: y(:), y_plus(:), y_minus(:)
+      real(qp), allocatable :: rate(:), rate_plus(:), rate_minus(:), drate(:)
+      real(qp), allocatable :: dydt(:), dydt_plus(:), dydt_minus(:), drate_dydt(:)
+      real(qp), allocatable :: dfdy(:,:)
+
+      allocate( &
+         y(size(n% y)), y_plus(size(n% y)), y_minus(size(n% y)), &
+         rate(size(n% rate_screened)), rate_plus(size(n% rate_screened)), &
+         rate_minus(size(n% rate_screened)), drate(size(n% rate_screened)), &
+         dydt(size(n% y)), dydt_plus(size(n% y)), dydt_minus(size(n% y)), &
+         drate_dydt(size(n% y)), dfdy(size(n% y),size(n% y)))
+
+      do i = 1, size(y)
+         y(i) = 0.2_qp + 0.01_qp*real(mod(11*i,17),kind=qp)
+      end do
+      do i = 1, size(rate)
+         rate(i) = 0.75_qp + 0.01_qp*real(mod(19*i,43),kind=qp)
+      end do
+
+      call eval_approx21_reduced_flows(y, rate, temp, dydt, dfdy)
+      max_error = 0.0_qp
+      i_max = 0
+      j_max = 0
+      do j = 1, size(y)
+         step = rel_step*max(1.0_qp,abs(y(j)))
+         y_plus = y
+         y_minus = y
+         y_plus(j) = y_plus(j) + step
+         y_minus(j) = y_minus(j) - step
+         call eval_approx21_reduced_flows(y_plus, rate, temp, dydt_plus)
+         call eval_approx21_reduced_flows(y_minus, rate, temp, dydt_minus)
+         do i = 1, size(y)
+            deriv = (dydt_plus(i) - dydt_minus(i))/(2.0_qp*step)
+            scale = max(1.0_qp,abs(deriv),abs(dfdy(i,j)))
+            error = abs(deriv - dfdy(i,j))/scale
+            if (error > max_error) then
+               max_error = error
+               i_max = i
+               j_max = j
+            end if
+         end do
+      end do
+      if (max_error > tolerance) then
+         write (*, '(a,2i6,1x,es14.6)') &
+            'bad approx21 reduced-flow composition derivative', i_max, j_max, max_error
+         call mesa_error(__FILE__, __LINE__)
+      end if
+
+      max_error = 0.0_qp
+      i_max = 0
+      j_max = 0
+      do j = 1, size(rate)
+         drate = 0.0_qp
+         drate(j) = 1.0_qp
+         call eval_approx21_reduced_flows( &
+            y, rate, temp, dydt, drate=drate, reduced_drate=drate_dydt)
+         step = rel_step*max(1.0_qp,abs(rate(j)))
+         rate_plus = rate
+         rate_minus = rate
+         rate_plus(j) = rate_plus(j) + step
+         rate_minus(j) = rate_minus(j) - step
+         call eval_approx21_reduced_flows(y, rate_plus, temp, dydt_plus)
+         call eval_approx21_reduced_flows(y, rate_minus, temp, dydt_minus)
+         do i = 1, size(y)
+            deriv = (dydt_plus(i) - dydt_minus(i))/(2.0_qp*step)
+            scale = max(1.0_qp,abs(deriv),abs(drate_dydt(i)))
+            error = abs(deriv - drate_dydt(i))/scale
+            if (error > max_error) then
+               max_error = error
+               i_max = i
+               j_max = j
+            end if
+         end do
+      end do
+      if (max_error > tolerance) then
+         write (*, '(a,2i6,1x,es14.6)') &
+            'bad approx21 reduced-flow rate derivative', i_max, j_max, max_error
+         call mesa_error(__FILE__, __LINE__)
+      end if
+
+      deallocate( &
+         y, y_plus, y_minus, rate, rate_plus, rate_minus, drate, &
+         dydt, dydt_plus, dydt_minus, drate_dydt, dfdy)
+   end subroutine test_approx21_reduced_flow_derivatives
+
    subroutine show_results( &
       g, n, logT, logRho, species, num_reactions, xin, &
       eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, &

net/test/test_output

@@ -194,19 +194,19 @@
  test_one_zone_burn_small_net
                                       number of species          21
                                  large final abundances    1.0000000000000000D-02
-                                                    c12           1    9.8714342254068388D-01
-                                                    he4           2    1.2802186836379657D-02
+                                                    c12           1    9.8714342211786166D-01
+                                                    he4           2    1.2802187260317500D-02
 
-                                                   xsum    1.0000000000006330D+00
+                                                   xsum    1.0000000000006328D+00
 
 
  test_one_zone_burn_const_P
                                       number of species          21
                                  large final abundances    1.0000000000000000D-02
-                                                    o16           1    6.9559810129898159D-01
-                                                   si28           2    2.0468825030127175D-01
-                                                    s32           3    4.4368437194011109D-02
-                                                   mg24           4    3.7858761619376645D-02
+                                                    o16           1    6.9559810129898514D-01
+                                                   si28           2    2.0468825030127566D-01
+                                                    s32           3    4.4368437194014287D-02
+                                                   mg24           4    3.7858761619371552D-02
 
                                                    xsum    1.0000000000000000D+00
 

rates/private/rates_support.f90

@@ -18,7 +18,7 @@
 ! ***********************************************************************
 
       module rates_support
-      use const_def, only: dp, use_mesa_temp_cache, missing_value, ln10
+      use const_def, only: dp, qp, use_mesa_temp_cache, missing_value, ln10
       use math_lib
       use rates_def
       use utils_lib, only: mv, switch_str, mesa_error
@@ -32,12 +32,12 @@ module rates_support
 
       subroutine do_get_raw_rates( &
             num_reactions, reaction_id, rattab, rattab_f1, nT8s, &
-            ye, logtemp_in, btemp, bden, raw_rate_factor, logttab, &
+            ye, logtemp_in, btemp, bden, raw_rate_factor, ttab, logttab, &
             rate_raw, rate_raw_dT, rate_raw_dRho, ierr)
          integer, intent(in) :: num_reactions, reaction_id(:), nT8s
          real(dp), intent(in) ::  &
             ye, logtemp_in, btemp, bden, raw_rate_factor(:),  &
-            rattab(:,:), logttab(:)
+            rattab(:,:), ttab(:), logttab(:)
          real(dp), pointer, intent(in) :: rattab_f1(:)
          real(dp), intent(inout), dimension(:) :: rate_raw, rate_raw_dT, rate_raw_dRho
          integer, intent(out) :: ierr
@@ -155,7 +155,13 @@ subroutine get_rates_from_table(r1, r2)
             do while (logtemp > logttab(k+1) .and. k+1 < nrattab)
                k = k+1
             end do
-            dt = logtemp - logttab(k)
+            if (logtemp_in >= max_safe_logT_for_rates) then
+               dt = logtemp - logttab(k)
+            else
+               ! Preserve local changes lost by subtracting absolute logarithms.
+               dt = real(log(real(btemp,kind=qp)/real(ttab(k),kind=qp))/ &
+                  real(ln10,kind=qp),kind=dp)
+            end if
 
             do i = r1,r2
 

rates/public/rates_lib.f90

@@ -769,18 +769,18 @@ end subroutine eval_weak_reaction_info
 
       subroutine eval_using_rate_tables( &
             num_reactions, reaction_id, rattab, rattab_f1, nT8s, &
-            ye, logtemp, btemp, bden, raw_rate_factor, logttab, &
+            ye, logtemp, btemp, bden, raw_rate_factor, ttab, logttab, &
             rate_raw, rate_raw_dT, rate_raw_dRho, ierr)
          use rates_support, only : do_get_raw_rates
          integer, intent(in) :: num_reactions, reaction_id(:), nT8s
          real(dp), intent(in) ::  &
             ye, logtemp, btemp, bden, raw_rate_factor(:),  &
-            rattab(:,:), logttab(:)
+            rattab(:,:), ttab(:), logttab(:)
          real(dp), pointer :: rattab_f1(:)
          real(dp), intent(out), dimension(:) :: rate_raw, rate_raw_dT, rate_raw_dRho
          integer, intent(out) :: ierr
          call do_get_raw_rates(num_reactions, reaction_id, rattab, rattab_f1, nT8s, &
-               ye, logtemp, btemp, bden, raw_rate_factor, logttab, &
+               ye, logtemp, btemp, bden, raw_rate_factor, ttab, logttab, &
                rate_raw, rate_raw_dT, rate_raw_dRho, ierr)
       end subroutine eval_using_rate_tables
 

star/defaults/controls.defaults

@@ -9149,6 +9149,8 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       ! ignore_species_in_max_correction
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      ! solver_use_T_reference
+      ! ~~~~~~~~~~~~~~~~~~~~~~
       ! num_times_solver_reuse_mtx
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -9165,6 +9167,13 @@
     tiny_corr_factor = 2
     ignore_min_corr_coeff_for_scale_max_correction = .false.
     ignore_species_in_max_correction = .false.
+
+      ! During solver iterations, calculate the trial temperature from the
+      ! temperature at the start of the solver call and the logarithmic
+      ! temperature correction. This avoids losing small corrections when
+      ! adding them to the absolute value of lnT.
+
+    solver_use_T_reference = .true.
     num_times_solver_reuse_mtx = 0
 
 

star/defaults/controls_dev.defaults

@@ -24,6 +24,16 @@
 
       use_TDC_Y_face_seeded_newton = .false.
 
+      ! solver_damping_on_retry
+      ! ~~~~~~~~~~~~~~~~~~~~~~~
+      ! If true, enable the global correction coefficient limit immediately
+      ! on retries. If false, retries outside relaxation use the same correction
+      ! coefficients as the original timestep. Relaxation retries remain damped.
+      !
+      ! ::
+
+      solver_damping_on_retry = .true.
+
       ! Hydro matrix solver options
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
       ! hydro_matrix_solver selects the mesa-star hydro matrix solve.

star/private/alloc.f90

@@ -1188,6 +1188,9 @@ subroutine star_info_arrays(s, c_in, action_in, ierr)
             call do2(s% xa_sub_xa_start, c% xa_sub_xa_start, species, 'xa_sub_xa_start')
             if (failed('xa_sub_xa_start')) exit
 
+            call do1(s% solver_T_reference, c% solver_T_reference)
+            if (failed('solver_T_reference')) exit
+
             call do1(s% lnd_start, c% lnd_start)
             if (failed('lnd_start')) exit
             call do1(s% lnPgas_start, c% lnPgas_start)

star/private/ctrls_io.f90

@@ -361,7 +361,7 @@ module ctrls_io
     P_theta_for_velocity_time_centering, L_theta_for_velocity_time_centering, &
     max_logT_for_include_P_and_L_in_velocity_time_centering, &
     steps_before_use_TDC, use_P_d_1_div_rho_form_of_work_when_time_centering_velocity, compare_TDC_to_MLT, &
-    use_TDC_Y_face_seeded_newton, &
+    use_TDC_Y_face_seeded_newton, solver_damping_on_retry, &
     hydro_matrix_solver, &
     remesh_for_TDC_pulsations_log_core_zoning, velocity_logT_lower_bound, &
     max_dt_yrs_for_velocity_logT_lower_bound, velocity_tau_lower_bound, velocity_q_upper_bound, &
@@ -403,6 +403,7 @@ module ctrls_io
     tiny_corr_coeff_limit, scale_correction_norm, corr_param_factor, num_times_solver_reuse_mtx, &
     scale_max_correction, ignore_min_corr_coeff_for_scale_max_correction, &
     ignore_too_large_correction, ignore_species_in_max_correction, &
+    solver_use_T_reference, &
     corr_norm_jump_limit, max_corr_jump_limit, resid_norm_jump_limit, max_resid_jump_limit, RSP2_use_mass_interp_face_values, &
     corr_coeff_limit, tiny_corr_factor, solver_test_partials_call_number, solver_test_partials_iter_number, &
     max_tries1, solver_max_tries_before_reject, max_tries_for_retry, max_tries_after_5_retries, solver_test_partials_sink_name, &
@@ -2017,6 +2018,7 @@ subroutine store_controls(s, ierr)
  s% ignore_min_corr_coeff_for_scale_max_correction = ignore_min_corr_coeff_for_scale_max_correction
  s% ignore_too_large_correction = ignore_too_large_correction
  s% ignore_species_in_max_correction = ignore_species_in_max_correction
+ s% solver_use_T_reference = solver_use_T_reference
 
  s% corr_norm_jump_limit = corr_norm_jump_limit
  s% max_corr_jump_limit = max_corr_jump_limit
@@ -2105,6 +2107,7 @@ subroutine store_controls(s, ierr)
  s% include_mlt_in_velocity_time_centering = include_mlt_in_velocity_time_centering
  s% compare_TDC_to_MLT = compare_TDC_to_MLT
  s% use_TDC_Y_face_seeded_newton = use_TDC_Y_face_seeded_newton
+ s% solver_damping_on_retry = solver_damping_on_retry
  s% hydro_matrix_solver = hydro_matrix_solver
  s% TDC_hydro_use_mass_interp_face_values = TDC_hydro_use_mass_interp_face_values
  s% TDC_hydro_nz = TDC_hydro_nz
@@ -3731,6 +3734,7 @@ subroutine set_controls_for_writing(s, ierr)
  ignore_min_corr_coeff_for_scale_max_correction = s% ignore_min_corr_coeff_for_scale_max_correction
  ignore_too_large_correction = s% ignore_too_large_correction
  ignore_species_in_max_correction = s% ignore_species_in_max_correction
+ solver_use_T_reference = s% solver_use_T_reference
 
  corr_norm_jump_limit = s% corr_norm_jump_limit
  max_corr_jump_limit = s% max_corr_jump_limit
@@ -3819,6 +3823,7 @@ subroutine set_controls_for_writing(s, ierr)
  include_mlt_in_velocity_time_centering = s% include_mlt_in_velocity_time_centering
  compare_TDC_to_MLT = s% compare_TDC_to_MLT
  use_TDC_Y_face_seeded_newton = s% use_TDC_Y_face_seeded_newton
+ solver_damping_on_retry = s% solver_damping_on_retry
  hydro_matrix_solver = s% hydro_matrix_solver
  TDC_hydro_use_mass_interp_face_values = s% TDC_hydro_use_mass_interp_face_values
  TDC_hydro_nz = s% TDC_hydro_nz

star/private/evolve.f90

@@ -672,6 +672,8 @@ integer function do_step_part2(id, first_try)
             end if
 
             skip_global_corr_coeff_limit = (first_try .or. &
+                (.not. s% doing_relax .and. .not. s% solver_damping_on_retry .and. &
+                 s% retry_cnt > 0) .or. &
                 s% model_number_for_last_retry /= s% model_number)  ! last alternative is for redo's
 
             s% doing_struct_burn_mix = .true.