Skip to content

Fix accurate AD summation in residual equations #1035

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Debraheem wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix accurate AD summation in residual equations #1035

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 35 additions & 10 deletions star/private/hydro_energy.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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'
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
15 changes: 11 additions & 4 deletions star/private/hydro_momentum.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
8 changes: 6 additions & 2 deletions star/private/hydro_riemann.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
18 changes: 15 additions & 3 deletions star/private/hydro_rsp2.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
1 change: 0 additions & 1 deletion star/private/hydro_rsp2_support.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 0 additions & 1 deletion star/private/tdc_hydro.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 0 additions & 1 deletion star/private/tdc_hydro_support.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading