From 746d93a57c5937e9273b82ed0e5779fb8436d21e Mon Sep 17 00:00:00 2001 From: Ebraheem Farag <63124736+Debraheem@users.noreply.github.com> Date: Wed, 8 Jul 2026 23:10:52 -0400 Subject: [PATCH] Fix accurate AD summation in residual equations --- star/private/hydro_energy.f90 | 45 ++++++++++++++++++++++------- star/private/hydro_momentum.f90 | 15 +++++++--- star/private/hydro_riemann.f90 | 8 +++-- star/private/hydro_rsp2.f90 | 18 ++++++++++-- star/private/hydro_rsp2_support.f90 | 1 - star/private/tdc_hydro.f90 | 1 - star/private/tdc_hydro_support.f90 | 1 - 7 files changed, 67 insertions(+), 22 deletions(-) diff --git a/star/private/hydro_energy.f90 b/star/private/hydro_energy.f90 index b22d3c8013..1dc228af38 100644 --- a/star/private/hydro_energy.f90 +++ b/star/private/hydro_energy.f90 @@ -100,15 +100,29 @@ subroutine get1_energy_eqn( & ! eps_WD_sedimentation, eps_diffusion, eps_pre_mix, eps_phase_separation ! sum terms in esum_ad using accurate_auto_diff_real_star_order1 if (eps_grav_form) then ! for this case, dwork_dm doesn't include work by P since that is in eps_grav - esum_ad = - dL_dm_ad + sources_ad + & - others_ad - d_turbulent_energy_dt_ad - dwork_dm_ad + eps_grav_ad + esum_ad = -dL_dm_ad + esum_ad = esum_ad + sources_ad + esum_ad = esum_ad + others_ad + esum_ad = esum_ad - d_turbulent_energy_dt_ad + esum_ad = esum_ad - dwork_dm_ad + esum_ad = esum_ad + eps_grav_ad else if (s% using_velocity_time_centering .and. & s% use_P_d_1_div_rho_form_of_work_when_time_centering_velocity) then - esum_ad = - dL_dm_ad + sources_ad + & - others_ad - d_turbulent_energy_dt_ad - dwork_dm_ad - de_dt_ad + esum_ad = -dL_dm_ad + esum_ad = esum_ad + sources_ad + esum_ad = esum_ad + others_ad + esum_ad = esum_ad - d_turbulent_energy_dt_ad + esum_ad = esum_ad - dwork_dm_ad + esum_ad = esum_ad - de_dt_ad else - esum_ad = - dL_dm_ad + sources_ad + & - others_ad - d_turbulent_energy_dt_ad - dwork_dm_ad - dke_dt_ad - dpe_dt_ad - de_dt_ad + esum_ad = -dL_dm_ad + esum_ad = esum_ad + sources_ad + esum_ad = esum_ad + others_ad + esum_ad = esum_ad - d_turbulent_energy_dt_ad + esum_ad = esum_ad - dwork_dm_ad + esum_ad = esum_ad - dke_dt_ad + esum_ad = esum_ad - dpe_dt_ad + esum_ad = esum_ad - de_dt_ad end if resid_ad = esum_ad ! convert back to auto_diff_real_star_order1 s% ergs_error(k) = -dm*dt*resid_ad%val ! save ergs_error before scaling @@ -222,6 +236,7 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad type(auto_diff_real_star_order1) :: & eps_nuc_ad, non_nuc_neu_ad, extra_heat_ad, Eq_ad, RTI_diffusion_ad, & v_00, v_p1, drag_force, drag_energy + type(accurate_auto_diff_real_star_order1) :: sources_sum_ad include 'formats' ierr = 0 @@ -301,7 +316,13 @@ subroutine setup_sources_and_others(ierr) ! sources_ad, others_ad end if end if - sources_ad = eps_nuc_ad - non_nuc_neu_ad + extra_heat_ad + Eq_ad + RTI_diffusion_ad + drag_energy + sources_sum_ad = eps_nuc_ad + sources_sum_ad = sources_sum_ad - non_nuc_neu_ad + sources_sum_ad = sources_sum_ad + extra_heat_ad + sources_sum_ad = sources_sum_ad + Eq_ad + sources_sum_ad = sources_sum_ad + RTI_diffusion_ad + sources_sum_ad = sources_sum_ad + drag_energy + sources_ad = sources_sum_ad sources_ad%val = sources_ad%val + s% irradiation_heat(k) @@ -560,7 +581,6 @@ end subroutine get1_energy_eqn subroutine eval_dwork(s, k, skip_P, dwork_ad, dwork, & d_dwork_dxam1, d_dwork_dxa00, d_dwork_dxap1, ierr) - use accurate_sum_auto_diff_star_order1 use auto_diff_support use star_utils, only: calc_Ptot_ad_tw type (star_info), pointer :: s @@ -630,6 +650,7 @@ subroutine eval1_work(s, k, skip_Peos, & P_face_ad, A_times_v_face_ad, mlt_Pturb_ad, & PtrbR_ad, PtrbL_ad, PvscL_ad, PvscR_ad, Ptrb_div_etrb, PL_ad, PR_ad, & Peos_ad, Ptrb_ad, Pvsc_ad, extra_P + type(accurate_auto_diff_real_star_order1) :: P_face_sum_ad logical :: test_partials integer :: j include 'formats' @@ -759,7 +780,12 @@ subroutine eval1_work(s, k, skip_Peos, & if (s% mlt_Pturb_factor > 0d0 .and. s% mlt_vc_old(k) > 0d0) & mlt_Pturb_ad = s% mlt_Pturb_factor*pow2(s% mlt_vc_old(k))*get_rho_face(s,k)/3d0 - P_face_ad = Peos_ad + Pvsc_ad + Ptrb_ad + mlt_Pturb_ad + extra_P + P_face_sum_ad = Peos_ad + P_face_sum_ad = P_face_sum_ad + Pvsc_ad + P_face_sum_ad = P_face_sum_ad + Ptrb_ad + P_face_sum_ad = P_face_sum_ad + mlt_Pturb_ad + P_face_sum_ad = P_face_sum_ad + extra_P + P_face_ad = P_face_sum_ad end if @@ -823,7 +849,6 @@ end subroutine eval1_A_times_v_face_ad subroutine eval_simple_PdV_work( & s, k, skip_P, dwork_ad, dwork, d_dwork_dxa00, ierr) - use accurate_sum_auto_diff_star_order1 use auto_diff_support use star_utils, only: calc_Ptot_ad_tw type (star_info), pointer :: s diff --git a/star/private/hydro_momentum.f90 b/star/private/hydro_momentum.f90 index ef8e620a3e..c156b26be6 100644 --- a/star/private/hydro_momentum.f90 +++ b/star/private/hydro_momentum.f90 @@ -138,8 +138,11 @@ subroutine get1_momentum_eqn( & RTI_terms_dm_div_A_ad = RTI_terms_ad*dm_div_A_ad ! sum terms in residual_sum_ad using accurate_auto_diff_real_star_order1 - residual_sum_ad = & - other_dm_div_A_ad + grav_dm_div_A_ad - dPtot_ad - d_mlt_Pturb_ad + RTI_terms_dm_div_A_ad + residual_sum_ad = other_dm_div_A_ad + residual_sum_ad = residual_sum_ad + grav_dm_div_A_ad + residual_sum_ad = residual_sum_ad - dPtot_ad + residual_sum_ad = residual_sum_ad - d_mlt_Pturb_ad + residual_sum_ad = residual_sum_ad + RTI_terms_dm_div_A_ad resid1_ad = residual_sum_ad ! convert back to auto_diff_real_star_order1 resid_ad = resid1_ad*iPtotavg_ad ! scaling @@ -382,6 +385,7 @@ subroutine expected_non_HSE_term( & integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: extra_ad, v_00, & drag + type(accurate_auto_diff_real_star_order1) :: other_sum_ad real(dp) :: accel, d_accel_dv logical :: test_partials, local_v_flag @@ -437,7 +441,11 @@ subroutine expected_non_HSE_term( & end if ! v_flag - other_ad = extra_ad - accel_ad + drag + Uq_ad + other_sum_ad = extra_ad + other_sum_ad = other_sum_ad - accel_ad + other_sum_ad = other_sum_ad + drag + other_sum_ad = other_sum_ad + Uq_ad + other_ad = other_sum_ad other = other_ad%val !test_partials = (k == s% solver_test_partials_k) @@ -458,7 +466,6 @@ subroutine get_dPtot_face_info(s, k, P_surf_ad, & dPtot_ad, dPtot, d_dPtot_dxam1, d_dPtot_dxa00, & iPtotavg_ad, iPtotavg, d_iPtotavg_dxam1, d_iPtotavg_dxa00, ierr) use star_utils, only: calc_Ptot_ad_tw - use accurate_sum_auto_diff_star_order1 use auto_diff_support type (star_info), pointer :: s integer, intent(in) :: k diff --git a/star/private/hydro_riemann.f90 b/star/private/hydro_riemann.f90 index ca7e83578c..3516f20cef 100644 --- a/star/private/hydro_riemann.f90 +++ b/star/private/hydro_riemann.f90 @@ -131,8 +131,12 @@ subroutine do1_dudt_eqn( & if (ierr /= 0) return end if - sum_ad = flux_in_ad - flux_out_ad + & - geometry_source_ad + gravity_source_ad + diffusion_source_ad + Uq_cell + sum_ad = flux_in_ad + sum_ad = sum_ad - flux_out_ad + sum_ad = sum_ad + geometry_source_ad + sum_ad = sum_ad + gravity_source_ad + sum_ad = sum_ad + diffusion_source_ad + sum_ad = sum_ad + Uq_cell dudt_expected_ad = sum_ad dudt_expected_ad = dudt_expected_ad/dm diff --git a/star/private/hydro_rsp2.f90 b/star/private/hydro_rsp2.f90 index 8bd134a00d..08cbe00c59 100644 --- a/star/private/hydro_rsp2.f90 +++ b/star/private/hydro_rsp2.f90 @@ -118,6 +118,7 @@ subroutine do1_rsp2_L_eqn(s, k, nvar, ierr) integer, intent(out) :: ierr type(auto_diff_real_star_order1) :: & L_expected, L_actual,resid + type(accurate_auto_diff_real_star_order1) :: L_sum real(dp) :: scale, residual, L_start_max logical :: test_partials include 'formats' @@ -132,7 +133,10 @@ subroutine do1_rsp2_L_eqn(s, k, nvar, ierr) ierr = 0 !L_expected = compute_L_face(s, k, ierr) !if (ierr /= 0) return - L_expected = s% Lr_ad(k) + s% Lc_ad(k) + s% Lt_ad(k) + L_sum = s% Lr_ad(k) + L_sum = L_sum + s% Lc_ad(k) + L_sum = L_sum + s% Lt_ad(k) + L_expected = L_sum L_actual = wrap_L_00(s, k) L_start_max = maxval(s% L_start(1:s% nz)) scale = 1d0/L_start_max @@ -317,7 +321,11 @@ subroutine do1_turbulent_energy_eqn(s, k, nvar, ierr) call setup_dt_Eq_ad(ierr); if (ierr /= 0) return ! erg g^-1 call set_energy_eqn_scal(s, k, scal, ierr); if (ierr /= 0) return ! 1/(erg g^-1 s^-1) ! sum terms in esum_ad using accurate_auto_diff_real_star_order1 - esum_ad = d_turbulent_energy_ad + Ptrb_dV_ad + dt_dLt_dm_ad - dt_C_ad - dt_Eq_ad ! erg g^-1 + esum_ad = d_turbulent_energy_ad + esum_ad = esum_ad + Ptrb_dV_ad + esum_ad = esum_ad + dt_dLt_dm_ad + esum_ad = esum_ad - dt_C_ad + esum_ad = esum_ad - dt_Eq_ad ! erg g^-1 resid_ad = esum_ad if (k==-35 .and. s% solver_iter == 1) then @@ -971,6 +979,7 @@ subroutine compute_L_terms(s, k, L, Lr, Lc, Lt, ierr) type (star_info), pointer, intent(in) :: s integer, intent(in) :: k type(auto_diff_real_star_order1), intent(out) :: L, Lr, Lc, Lt + type(accurate_auto_diff_real_star_order1) :: L_sum integer, intent(out) :: ierr include 'formats' ierr = 0 @@ -993,7 +1002,10 @@ subroutine compute_L_terms(s, k, L, Lr, Lc, Lt, ierr) Lt = compute_Lt(s, k, ierr) if (ierr /= 0) return end if - L = Lr + Lc + Lt + L_sum = Lr + L_sum = L_sum + Lc + L_sum = L_sum + Lt + L = L_sum s% Lr_ad(k) = Lr s% Lc_ad(k) = Lc s% Lt_ad(k) = Lt diff --git a/star/private/hydro_rsp2_support.f90 b/star/private/hydro_rsp2_support.f90 index 0d6aaabcdf..ccb0df188a 100644 --- a/star/private/hydro_rsp2_support.f90 +++ b/star/private/hydro_rsp2_support.f90 @@ -24,7 +24,6 @@ module hydro_rsp2_support use utils_lib, only: is_bad use auto_diff use auto_diff_support - use accurate_sum_auto_diff_star_order1 use star_utils implicit none diff --git a/star/private/tdc_hydro.f90 b/star/private/tdc_hydro.f90 index 3b10b394c2..3491aca94a 100644 --- a/star/private/tdc_hydro.f90 +++ b/star/private/tdc_hydro.f90 @@ -24,7 +24,6 @@ module tdc_hydro use utils_lib, only: is_bad use auto_diff use auto_diff_support - use accurate_sum_auto_diff_star_order1 use star_utils implicit none diff --git a/star/private/tdc_hydro_support.f90 b/star/private/tdc_hydro_support.f90 index 161d1f6791..820b53e514 100644 --- a/star/private/tdc_hydro_support.f90 +++ b/star/private/tdc_hydro_support.f90 @@ -24,7 +24,6 @@ module tdc_hydro_support use utils_lib, only: is_bad use auto_diff use auto_diff_support - use accurate_sum_auto_diff_star_order1 use star_utils implicit none