From b935998f25b2383ee720e46c7f6cbec1f7eba44a Mon Sep 17 00:00:00 2001 From: Debraheem Date: Sat, 11 Jan 2025 19:07:13 -0500 Subject: [PATCH 01/10] [ci optional] add and test use_eos_partials_for_Brunt --- star/defaults/controls.defaults | 9 ++ star/private/brunt.f90 | 176 ++++++++++++++++++++++++++++ star/private/ctrls_io.f90 | 6 +- star_data/private/star_controls.inc | 2 +- 4 files changed, 190 insertions(+), 3 deletions(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 3d257a93c..92013b656 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -7713,6 +7713,15 @@ calculate_Brunt_B = .true. calculate_Brunt_N2 = .true. + ! calculate_Brunt_N2 + ! ~~~~~~~~~~~~~~~~~~ + + ! If true calculate B using equation 6 in MESA VI, otherwise use equation 8 + + + ! :: + + use_eos_partials_for_Brunt = .true. ! brunt_N2_coefficient ! ~~~~~~~~~~~~~~~~~~~~ diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index 2af4acdd1..c0154e017 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -69,6 +69,17 @@ subroutine do_brunt_B(s,nzlo,nzhi,ierr) s% retry_message = 'failed in other_brunt' if (s% report_ierr) write(*, *) s% retry_message end if + + else if (s% use_eos_partials_for_Brunt) then + ! ------------------------------------------------------------ + ! chain-rule approach for calculating B from eos partials, see MESA II equation 6. + ! ------------------------------------------------------------ + call do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) + if (ierr /= 0) then + s% retry_message = 'failed in do_brunt_B_eos_partials_form' + if (s% report_ierr) write(*, *) s% retry_message + end if + else call do_brunt_B_MHM_form(s,nzlo,nzhi,ierr) if (ierr /= 0) then @@ -340,4 +351,169 @@ subroutine get_brunt_B(s, species, nz, k, T_face, rho_face, chiT_face, chiRho_fa end subroutine get_brunt_B + subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) + ! Brassard from Mike Montgomery (MHM) + use star_utils, only: get_face_values + use interp_1d_def + + type (star_info), pointer :: s + integer, intent(in) :: nzlo, nzhi + integer, intent(out) :: ierr + + + real(dp), allocatable, dimension(:) :: T_face, rho_face, chiT_face, chiRho_face + ! real(dp) :: mass_corr_factor, delta_lnP, delta_lnMbar, B_cell_centered + integer :: nz, species, k, op_err + logical, parameter :: dbg = .false. + + include 'formats' + + ierr = 0 + + nz = s% nz + species = s% species + + allocate(T_face(nz), rho_face(nz), chiT_face(nz), chiRho_face(nz)) ! ,B_cell_centered(nz) + + call get_face_values(s, s% chiT, chiT_face, ierr) + if (ierr /= 0) return + + call get_face_values(s, s% chiRho, chiRho_face, ierr) + if (ierr /= 0) return + + call get_face_values(s, s% T, T_face, ierr) + if (ierr /= 0) return + + call get_face_values(s, s% rho, rho_face, ierr) + if (ierr /= 0) return + + !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) + do k=nzlo,nzhi + op_err = 0 + call get_brunt_B(& + s, species, nz, k, T_face(k), rho_face(k), chiT_face(k), chiRho_face(k), op_err) + if (op_err /= 0) ierr = op_err + end do + + ! Probably not needed. + ! Because I initially thought B was returned on the cell not face, the items below + ! were intended to apply corrections and move back to face centered. + !!!$OMP END PARALLEL DO + ! + !! store center cell info for B into new variable for cell centered_B + !B_cell_centered = s% brunt_B + !! Average cell_centered B onto faces and return to s% Brunt_B + !call get_face_values(s, B_cell_centered, s% brunt_B, ierr) + !if (ierr /= 0) return + ! + !! Add mass correction term if enabled + !!$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) + !do k=nzlo,nzhi + !if (s% use_mass_corrections) then + ! delta_lnMbar = log(s% mass_correction(k-1)) - log(s% mass_correction(k)) + ! delta_lnP = s% lnPeos(k-1) - s% lnPeos(k) + ! mass_corr_factor = chiRho_face(k) * delta_lnMbar / (delta_lnP *chiT_face(k)) + ! s% brunt_B(k) = s% brunt_B(k) - mass_corr_factor + !end if + !end do + + !!!$OMP END PARALLEL DO + + + end subroutine do_brunt_B_eos_partials_form + + subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, chiT_face, chiRho_face, ierr) + use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnPgas + use eos_support, only: get_eos + + type (star_info), pointer :: s + integer, intent(in) :: species, nz, k + real(dp), intent(in) :: T_face, rho_face, chiT_face, chiRho_face + integer, intent(out) :: ierr + + real(dp) :: logRho_face, logT_face, Prad_face + real(dp), dimension(num_eos_basic_results) :: res, d_eos_dlnd, d_eos_dlnT + real(dp), dimension(num_eos_d_dxa_results, species) :: d_eos_dxa + real(dp) :: mass_corr_factor, delta_lnP, delta_lnMbar, Ppoint, dlnP_dm, alfa + real(dp) :: chiT, B_term, total_derivative_P_dX, total_derivative_X_dP + integer :: i + + logical, parameter :: dbg = .false. + + include 'formats' + + ierr = 0 + s% brunt_B(k) = 0 + if (k <= 1) return + + logT_face = log10(T_face) + logRho_face = log10(rho_face) + Prad_face = crad * T_face**4 / 3 + + if (is_bad_num(logT_face) .or. is_bad_num(logRho_face)) then + ierr = -1 + return + end if + + ! Call the EOS to get the required partial derivatives + call get_eos( & + s, 0, s% xa(:, k), & + rho_face, logRho_face, T_face, logT_face, & + res, d_eos_dlnd, d_eos_dlnT, & + d_eos_dxa, ierr) + if (ierr /= 0) return + + ! Extract chiT + chiT = d_eos_dlnT(i_lnPgas) + + ! Initialize B_term to accumulate the contribution from each species + B_term = 0.0 + + ! Compute the Brunt B composition term + do i = 1, species + ! Total derivative of P with respect to X_i + total_derivative_P_dX = d_eos_dxa(i_lnPgas, i) & + + d_eos_dlnd(i_lnPgas) * d_eos_dxa(i, i) & + + d_eos_dlnT(i_lnPgas) * d_eos_dxa(i_lnPgas, i) + + ! Total derivative of X_i with respect to P + total_derivative_X_dP = 1.0 / total_derivative_P_dX + + ! Accumulate the contribution to B_term + B_term = B_term - (d_eos_dxa(i_lnPgas, i) * total_derivative_X_dP) + end do + + ! Final calculation of B using chiT + s% brunt_B(k) = B_term / chiT + + + delta_lnP = s% lnPeos(k-1) - s% lnPeos(k) + if (delta_lnP > -1d-50) then + alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k)) + Ppoint = alfa*s% Peos(k) + (1-alfa)*s% Peos(k-1) + dlnP_dm = -s% cgrav(k)*s% m(k)/(pi4*pow4(s% r(k))*Ppoint) + delta_lnP = dlnP_dm*s% dm_bar(k) + end if + + ! add term accounting for the composition-related gradient in gravitational mass + if (s% use_mass_corrections) then + delta_lnMbar = log(s% mass_correction(k-1)) - log(s% mass_correction(k)) + s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP/chiT_face + end if + + ! Check for bad numbers + if (is_bad_num(s% brunt_B(k))) then + ierr = -1 + s% retry_message = 'bad num for brunt_B' + if (s% report_ierr) then + write(*,2) 's% brunt_B(k)', k, s% brunt_B(k) + write(*,2) 'chiT', k, chiT + write(*,2) 'B_term', k, B_term + call mesa_error(__FILE__, __LINE__, 'get_brunt_B_from_eos_partials') + end if + end if + end subroutine get_brunt_B_from_eos_partials + + + end module brunt diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 526da6265..adbdee981 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -151,8 +151,8 @@ module ctrls_io max_v_for_convection, max_q_for_convection_with_hydro_on, alpha_RTI_src_max_q, & max_v_div_cs_for_convection, max_abs_du_div_cs_for_convection, RSP_max_dt, RSP_relax_dm_tolerance, & calculate_Brunt_B, calculate_Brunt_N2, brunt_N2_coefficient, num_cells_for_smooth_brunt_B, & - threshold_for_smooth_brunt_B, min_magnitude_brunt_B, RSP_max_dt_times_min_rad_diff_time, & - min_overshoot_q, overshoot_alpha, RSP_target_steps_per_cycle, & + threshold_for_smooth_brunt_B, min_magnitude_brunt_B, use_eos_partials_for_Brunt, & + RSP_max_dt_times_min_rad_diff_time, min_overshoot_q, overshoot_alpha, RSP_target_steps_per_cycle, & RSP_max_num_periods, RSP_min_max_R_for_periods, RSP_min_deltaR_for_periods, & RSP_min_PERIOD_div_PERIODLIN, RSP_report_limit_dt, RSP_mode_for_setting_PERIODLIN, RSP_initial_dt_factor, & RSP_v_div_cs_threshold_for_dt_limit, RSP_max_dt_times_min_dr_div_cs, RSP_thetae, & @@ -1129,6 +1129,7 @@ subroutine store_controls(s, ierr) s% num_cells_for_smooth_brunt_B = num_cells_for_smooth_brunt_B s% threshold_for_smooth_brunt_B = threshold_for_smooth_brunt_B s% min_magnitude_brunt_B = min_magnitude_brunt_B + s% use_eos_partials_for_Brunt = use_eos_partials_for_Brunt s% min_overshoot_q = min_overshoot_q s% overshoot_alpha = overshoot_alpha @@ -2815,6 +2816,7 @@ subroutine set_controls_for_writing(s, ierr) brunt_N2_coefficient = s% brunt_N2_coefficient threshold_for_smooth_brunt_B = s% threshold_for_smooth_brunt_B min_magnitude_brunt_B = s% min_magnitude_brunt_B + use_eos_partials_for_Brunt = s% use_eos_partials_for_Brunt min_overshoot_q = s% min_overshoot_q overshoot_alpha = s% overshoot_alpha diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index 855cc6a5f..c48e79ac5 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -384,7 +384,7 @@ integer :: RTI_smooth_iterations ! brunt - logical :: calculate_Brunt_B, calculate_Brunt_N2 + logical :: calculate_Brunt_B, calculate_Brunt_N2, use_eos_partials_for_Brunt real(dp) :: min_magnitude_brunt_B, brunt_N2_coefficient integer :: num_cells_for_smooth_brunt_B real(dp) :: threshold_for_smooth_brunt_B From 6849d59d9bced61036ae737d29981260a9cb1e7d Mon Sep 17 00:00:00 2001 From: Debraheem Date: Sat, 11 Jan 2025 19:12:22 -0500 Subject: [PATCH 02/10] [ci skip] fix typo --- star/defaults/controls.defaults | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 92013b656..e082b6aa5 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -7716,7 +7716,7 @@ ! calculate_Brunt_N2 ! ~~~~~~~~~~~~~~~~~~ - ! If true calculate B using equation 6 in MESA VI, otherwise use equation 8 + ! If true calculate B using equation 6 in MESA II, otherwise use equation 8 ! :: From 11ccb7a135379f5ad11a850e5ef687705a62b52d Mon Sep 17 00:00:00 2001 From: Debraheem Date: Sat, 11 Jan 2025 22:23:33 -0500 Subject: [PATCH 03/10] fix equations --- star/private/brunt.f90 | 60 +++++++++++++----------------------------- 1 file changed, 19 insertions(+), 41 deletions(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index c0154e017..b4576f8ad 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -352,7 +352,6 @@ end subroutine get_brunt_B subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) - ! Brassard from Mike Montgomery (MHM) use star_utils, only: get_face_values use interp_1d_def @@ -363,6 +362,7 @@ subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) real(dp), allocatable, dimension(:) :: T_face, rho_face, chiT_face, chiRho_face ! real(dp) :: mass_corr_factor, delta_lnP, delta_lnMbar, B_cell_centered + real(dp), allocatable, dimension(:,:) :: xa_face! (species, k) integer :: nz, species, k, op_err logical, parameter :: dbg = .false. @@ -373,7 +373,14 @@ subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) nz = s% nz species = s% species - allocate(T_face(nz), rho_face(nz), chiT_face(nz), chiRho_face(nz)) ! ,B_cell_centered(nz) + allocate(T_face(nz), rho_face(nz), chiT_face(nz), chiRho_face(nz), xa_face(species, nz)) ! ,B_cell_centered(nz) + + ! can we paralleize this? + do k = 1, species + call get_face_values(s, s% xa(k, :), xa_face(k, :), ierr) + if (ierr /= 0) return + end do + !!! call get_face_values(s, s% chiT, chiT_face, ierr) if (ierr /= 0) return @@ -390,39 +397,16 @@ subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) !$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) do k=nzlo,nzhi op_err = 0 - call get_brunt_B(& - s, species, nz, k, T_face(k), rho_face(k), chiT_face(k), chiRho_face(k), op_err) + call get_brunt_B_from_eos_partials(& + s, species, nz, k, T_face(k), rho_face(k), chiT_face(k), chiRho_face(k), xa_face(:,:), op_err) if (op_err /= 0) ierr = op_err end do - - ! Probably not needed. - ! Because I initially thought B was returned on the cell not face, the items below - ! were intended to apply corrections and move back to face centered. - !!!$OMP END PARALLEL DO - ! - !! store center cell info for B into new variable for cell centered_B - !B_cell_centered = s% brunt_B - !! Average cell_centered B onto faces and return to s% Brunt_B - !call get_face_values(s, B_cell_centered, s% brunt_B, ierr) - !if (ierr /= 0) return - ! - !! Add mass correction term if enabled - !!$OMP PARALLEL DO PRIVATE(k,op_err) SCHEDULE(dynamic,2) - !do k=nzlo,nzhi - !if (s% use_mass_corrections) then - ! delta_lnMbar = log(s% mass_correction(k-1)) - log(s% mass_correction(k)) - ! delta_lnP = s% lnPeos(k-1) - s% lnPeos(k) - ! mass_corr_factor = chiRho_face(k) * delta_lnMbar / (delta_lnP *chiT_face(k)) - ! s% brunt_B(k) = s% brunt_B(k) - mass_corr_factor - !end if - !end do - - !!!$OMP END PARALLEL DO + !$OMP END PARALLEL DO end subroutine do_brunt_B_eos_partials_form - subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, chiT_face, chiRho_face, ierr) + subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, chiT_face, chiRho_face, xa_face, ierr) use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnPgas use eos_support, only: get_eos @@ -432,10 +416,11 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch integer, intent(out) :: ierr real(dp) :: logRho_face, logT_face, Prad_face + real(dp), dimension(species,nz), intent(in):: xa_face! (species, nz) real(dp), dimension(num_eos_basic_results) :: res, d_eos_dlnd, d_eos_dlnT real(dp), dimension(num_eos_d_dxa_results, species) :: d_eos_dxa real(dp) :: mass_corr_factor, delta_lnP, delta_lnMbar, Ppoint, dlnP_dm, alfa - real(dp) :: chiT, B_term, total_derivative_P_dX, total_derivative_X_dP + real(dp) :: chiT, B_term, spatial_derivative_dX_dlnP integer :: i logical, parameter :: dbg = .false. @@ -457,7 +442,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! Call the EOS to get the required partial derivatives call get_eos( & - s, 0, s% xa(:, k), & + s, 0, xa_face(:,k), & rho_face, logRho_face, T_face, logT_face, & res, d_eos_dlnd, d_eos_dlnT, & d_eos_dxa, ierr) @@ -471,16 +456,9 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! Compute the Brunt B composition term do i = 1, species - ! Total derivative of P with respect to X_i - total_derivative_P_dX = d_eos_dxa(i_lnPgas, i) & - + d_eos_dlnd(i_lnPgas) * d_eos_dxa(i, i) & - + d_eos_dlnT(i_lnPgas) * d_eos_dxa(i_lnPgas, i) - - ! Total derivative of X_i with respect to P - total_derivative_X_dP = 1.0 / total_derivative_P_dX - - ! Accumulate the contribution to B_term - B_term = B_term - (d_eos_dxa(i_lnPgas, i) * total_derivative_X_dP) + ! spatial derivative of X_i with respect to P + spatial_derivative_dX_dlnP = (s% xa(i,k) - s% xa(i,k-1))/(s% lnPeos(k) - s% lnPeos(k-1)) + B_term = B_term - (d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP) end do ! Final calculation of B using chiT From 53c1d442a336ccc3c5c03c87ffd6c2634870f085 Mon Sep 17 00:00:00 2001 From: EbF Date: Wed, 21 May 2025 13:02:42 -0400 Subject: [PATCH 04/10] fix brunt bug for flat eos --- star/private/brunt.f90 | 32 ++++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index b2b5a0b4e..e2edfd965 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -414,6 +414,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch real(dp), dimension(num_eos_d_dxa_results, species) :: d_eos_dxa real(dp) :: mass_corr_factor, delta_lnP, delta_lnMbar, Ppoint, dlnP_dm, alfa real(dp) :: chiT, B_term, spatial_derivative_dX_dlnP + !real(dp) :: num, deriv, den integer :: i logical, parameter :: dbg = .false. @@ -426,7 +427,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch logT_face = log10(T_face) logRho_face = log10(rho_face) - Prad_face = crad * T_face**4 / 3 + Prad_face = crad * pow4(T_face) / 3d0 if (is_bad_num(logT_face) .or. is_bad_num(logRho_face)) then ierr = -1 @@ -445,7 +446,23 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch chiT = d_eos_dlnT(i_lnPgas) ! Initialize B_term to accumulate the contribution from each species - B_term = 0.0 + B_term = 0d0 + + ! Compute pressure difference across adjacent cells + delta_lnP = s% lnPeos(k-1) - s% lnPeos(k) + + if (.false.) then ! for debugging, if we want to know when Peos is unchanging. + if (abs(delta_lnP) < 1d-30) then + write (*,'(A,I0,A,F12.5)') '*** ZERO ΔlnP at k=', k, ', ΔlnP=', delta_lnP + end if + end if + + if (abs(delta_lnP) < 1d-15) then + ! no resolvable composition gradient ⇒ B stays zero + s%brunt_B(k) = 0d0 + return + end if + ! Compute the Brunt B composition term do i = 1, species @@ -454,6 +471,17 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch B_term = B_term - (d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP) end do +! do i = 1, species +! num = s% xa(i,k) - s% xa(i,k-1) +! den = delta_lnP +! deriv = num/den +! if (is_bad_num(deriv)) then +! write (*,'(A,I0,A,I0,A,F12.5,A,F12.5)') & +! 'Bad dX/dlnP at k=', k, ', i=', i, ': num=', num, ', den=', den +! end if +! B_term = B_term - d_eos_dxa(i_lnPgas,i) * deriv +! end do + ! Final calculation of B using chiT s% brunt_B(k) = B_term / chiT From 57749fb5e85060833b2efb44ba92c983199a24ed Mon Sep 17 00:00:00 2001 From: EbF Date: Wed, 21 May 2025 23:59:49 -0400 Subject: [PATCH 05/10] further refactoring --- star/private/brunt.f90 | 63 +++++++++++++++++++----------------------- 1 file changed, 28 insertions(+), 35 deletions(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index e2edfd965..53fccbbac 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -412,9 +412,9 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch real(dp), dimension(species,nz), intent(in):: xa_face! (species, nz) real(dp), dimension(num_eos_basic_results) :: res, d_eos_dlnd, d_eos_dlnT real(dp), dimension(num_eos_d_dxa_results, species) :: d_eos_dxa - real(dp) :: mass_corr_factor, delta_lnP, delta_lnMbar, Ppoint, dlnP_dm, alfa - real(dp) :: chiT, B_term, spatial_derivative_dX_dlnP - !real(dp) :: num, deriv, den + real(dp) :: delta_lnMbar, Ppoint, dlnP_dm, alfa + real(dp) :: B_term, spatial_derivative_dX_dlnP, chiT + real(dp) :: delta_lnP_cell, comp, y, t integer :: i logical, parameter :: dbg = .false. @@ -423,7 +423,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ierr = 0 s% brunt_B(k) = 0 - if (k <= 1) return + if (k <= 1) return ! should we add a check for nz? logT_face = log10(T_face) logRho_face = log10(rho_face) @@ -447,57 +447,50 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! Initialize B_term to accumulate the contribution from each species B_term = 0d0 + comp = 0d0 ! compensator ! Compute pressure difference across adjacent cells - delta_lnP = s% lnPeos(k-1) - s% lnPeos(k) - - if (.false.) then ! for debugging, if we want to know when Peos is unchanging. - if (abs(delta_lnP) < 1d-30) then - write (*,'(A,I0,A,F12.5)') '*** ZERO ΔlnP at k=', k, ', ΔlnP=', delta_lnP - end if + delta_lnP_cell = s%lnPeos(k-1) - s%lnPeos(k) ! center difference + if (abs(delta_lnP_cell) < tiny(1d0)) then ! is this an okay check? + ! if the face pressure is flat, then the derivative is numerically zero + s% brunt_B(k) = 0d0 + return end if - if (abs(delta_lnP) < 1d-15) then - ! no resolvable composition gradient ⇒ B stays zero - s%brunt_B(k) = 0d0 - return - end if + ! Compute the Brunt B composition term +! do i = 1, species +! !X_face_p1 = xa_face(i, k) ! face between zones k and k+1 +! !X_face_m1 = xa_face(i, k-1) ! face between zones k-1 and k +! spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP_cell +! if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 +! B_term = B_term - d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP +! end do - ! Compute the Brunt B composition term do i = 1, species - ! spatial derivative of X_i with respect to P - spatial_derivative_dX_dlnP = (s% xa(i,k) - s% xa(i,k-1))/(s% lnPeos(k) - s% lnPeos(k-1)) - B_term = B_term - (d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP) + spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP_cell + if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 + y = -d_eos_dxa(i_lnPgas,i) * spatial_derivative_dX_dlnP - comp + t = B_term + y + comp = (t - B_term) - y + B_term = t end do -! do i = 1, species -! num = s% xa(i,k) - s% xa(i,k-1) -! den = delta_lnP -! deriv = num/den -! if (is_bad_num(deriv)) then -! write (*,'(A,I0,A,I0,A,F12.5,A,F12.5)') & -! 'Bad dX/dlnP at k=', k, ', i=', i, ': num=', num, ', den=', den -! end if -! B_term = B_term - d_eos_dxa(i_lnPgas,i) * deriv -! end do - ! Final calculation of B using chiT s% brunt_B(k) = B_term / chiT - - delta_lnP = s% lnPeos(k-1) - s% lnPeos(k) - if (delta_lnP > -1d-50) then + ! this block is always executed and assume HSE, is this an error? + if (delta_lnP_cell > -1d-50) then alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k)) Ppoint = alfa*s% Peos(k) + (1-alfa)*s% Peos(k-1) dlnP_dm = -s% cgrav(k)*s% m(k)/(pi4*pow4(s% r(k))*Ppoint) - delta_lnP = dlnP_dm*s% dm_bar(k) + delta_lnP_cell = dlnP_dm*s% dm_bar(k) end if ! add term accounting for the composition-related gradient in gravitational mass if (s% use_mass_corrections) then delta_lnMbar = log(s% mass_correction(k-1)) - log(s% mass_correction(k)) - s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP/chiT_face + s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP_cell/chiT_face end if ! Check for bad numbers From 0992ecb9ad4fa3e3dba9a1a9e55be3473a2b1127 Mon Sep 17 00:00:00 2001 From: EbF Date: Thu, 22 May 2025 00:22:12 -0400 Subject: [PATCH 06/10] fix mistakes, more cleanup --- star/private/brunt.f90 | 59 ++++++++++++++++++------------------------ 1 file changed, 25 insertions(+), 34 deletions(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index 53fccbbac..0a74c0926 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -414,7 +414,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch real(dp), dimension(num_eos_d_dxa_results, species) :: d_eos_dxa real(dp) :: delta_lnMbar, Ppoint, dlnP_dm, alfa real(dp) :: B_term, spatial_derivative_dX_dlnP, chiT - real(dp) :: delta_lnP_cell, comp, y, t + real(dp) :: delta_lnP !, comp, y, t integer :: i logical, parameter :: dbg = .false. @@ -443,54 +443,45 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch if (ierr /= 0) return ! Extract chiT - chiT = d_eos_dlnT(i_lnPgas) + chiT = d_eos_dlnT(i_lnPgas) ! should we be using chiT_face instead? ! Initialize B_term to accumulate the contribution from each species B_term = 0d0 comp = 0d0 ! compensator ! Compute pressure difference across adjacent cells - delta_lnP_cell = s%lnPeos(k-1) - s%lnPeos(k) ! center difference - if (abs(delta_lnP_cell) < tiny(1d0)) then ! is this an okay check? - ! if the face pressure is flat, then the derivative is numerically zero - s% brunt_B(k) = 0d0 - return + delta_lnP = s%lnPeos(k-1) - s%lnPeos(k) ! center difference + if (delta_lnP > -1d-50) then ! always applied. + alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k)) + Ppoint = alfa*s% Peos(k) + (1-alfa)*s% Peos(k-1) + dlnP_dm = -s% cgrav(k)*s% m(k)/(pi4*pow4(s% r(k))*Ppoint) + delta_lnP = dlnP_dm*s% dm_bar(k) end if - ! Compute the Brunt B composition term -! do i = 1, species -! !X_face_p1 = xa_face(i, k) ! face between zones k and k+1 -! !X_face_m1 = xa_face(i, k-1) ! face between zones k-1 and k -! spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP_cell -! if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 -! B_term = B_term - d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP -! end do - - - do i = 1, species - spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP_cell - if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 - y = -d_eos_dxa(i_lnPgas,i) * spatial_derivative_dX_dlnP - comp - t = B_term + y - comp = (t - B_term) - y - B_term = t - end do + ! Compute the Brunt B composition term + do i = 1, species + spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP + if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 + B_term = B_term - d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP + end do + +! might be necessary to use a compensator for large sums... +! do i = 1, species +! spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP +! if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 +! y = -d_eos_dxa(i_lnPgas,i) * spatial_derivative_dX_dlnP - comp +! t = B_term + y +! comp = (t - B_term) - y +! B_term = t +! end do ! Final calculation of B using chiT s% brunt_B(k) = B_term / chiT - ! this block is always executed and assume HSE, is this an error? - if (delta_lnP_cell > -1d-50) then - alfa = s% dq(k-1)/(s% dq(k-1) + s% dq(k)) - Ppoint = alfa*s% Peos(k) + (1-alfa)*s% Peos(k-1) - dlnP_dm = -s% cgrav(k)*s% m(k)/(pi4*pow4(s% r(k))*Ppoint) - delta_lnP_cell = dlnP_dm*s% dm_bar(k) - end if - ! add term accounting for the composition-related gradient in gravitational mass if (s% use_mass_corrections) then delta_lnMbar = log(s% mass_correction(k-1)) - log(s% mass_correction(k)) - s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP_cell/chiT_face + s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP/ChiT end if ! Check for bad numbers From b9f553faff2a639e1201a8fe6e70693255e453c0 Mon Sep 17 00:00:00 2001 From: EbF Date: Thu, 22 May 2025 00:23:10 -0400 Subject: [PATCH 07/10] remove uninitialized variable --- star/private/brunt.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index 0a74c0926..077b49665 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -447,7 +447,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! Initialize B_term to accumulate the contribution from each species B_term = 0d0 - comp = 0d0 ! compensator + !comp = 0d0 ! compensator ! Compute pressure difference across adjacent cells delta_lnP = s%lnPeos(k-1) - s%lnPeos(k) ! center difference From b9fc80e41f925ca56c2d73e71e083e2f5c6f79f4 Mon Sep 17 00:00:00 2001 From: EbF Date: Thu, 22 May 2025 00:30:44 -0400 Subject: [PATCH 08/10] open the throttle --- star/private/brunt.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index 077b49665..5dac8d8f8 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -461,7 +461,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! Compute the Brunt B composition term do i = 1, species spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP - if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 + !if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 B_term = B_term - d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP end do From 6a210a535200ee8fcd82e9f63802453759304ce2 Mon Sep 17 00:00:00 2001 From: EbF Date: Thu, 22 May 2025 01:16:16 -0400 Subject: [PATCH 09/10] [ci skip] clean up --- star/private/brunt.f90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index 5dac8d8f8..7231b6c59 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -354,7 +354,6 @@ subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) real(dp), allocatable, dimension(:) :: T_face, rho_face, chiT_face, chiRho_face - ! real(dp) :: mass_corr_factor, delta_lnP, delta_lnMbar, B_cell_centered real(dp), allocatable, dimension(:,:) :: xa_face! (species, k) integer :: nz, species, k, op_err logical, parameter :: dbg = .false. @@ -366,9 +365,9 @@ subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) nz = s% nz species = s% species - allocate(T_face(nz), rho_face(nz), chiT_face(nz), chiRho_face(nz), xa_face(species, nz)) ! ,B_cell_centered(nz) + allocate(T_face(nz), rho_face(nz), chiT_face(nz), chiRho_face(nz), xa_face(species, nz)) - ! can we paralleize this? + ! should we parallelize this? do k = 1, species call get_face_values(s, s% xa(k, :), xa_face(k, :), ierr) if (ierr /= 0) return @@ -461,14 +460,12 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! Compute the Brunt B composition term do i = 1, species spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP - !if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 B_term = B_term - d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP end do ! might be necessary to use a compensator for large sums... ! do i = 1, species ! spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP -! if (abs(spatial_derivative_dX_dlnP) < 1d-12) spatial_derivative_dX_dlnP = 0d0 ! y = -d_eos_dxa(i_lnPgas,i) * spatial_derivative_dX_dlnP - comp ! t = B_term + y ! comp = (t - B_term) - y From 7ab6d68926ad0fd17e19fad83fd99a36a3b9d9e5 Mon Sep 17 00:00:00 2001 From: EbF Date: Thu, 12 Jun 2025 16:04:06 -0400 Subject: [PATCH 10/10] add broken finite difference for eos partials --- star/private/brunt.f90 | 252 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 246 insertions(+), 6 deletions(-) diff --git a/star/private/brunt.f90 b/star/private/brunt.f90 index 7231b6c59..f05250f4b 100644 --- a/star/private/brunt.f90 +++ b/star/private/brunt.f90 @@ -20,7 +20,7 @@ module brunt use star_private_def - use const_def, only: dp, pi4, crad + use const_def, only: dp, pi4, crad, ln10 use utils_lib implicit none @@ -399,8 +399,9 @@ subroutine do_brunt_B_eos_partials_form(s, nzlo, nzhi, ierr) end subroutine do_brunt_B_eos_partials_form subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, chiT_face, chiRho_face, xa_face, ierr) - use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnPgas + use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnPgas, i_chiT use eos_support, only: get_eos + use chem_def, only: ihe4, ic12 type (star_info), pointer :: s integer, intent(in) :: species, nz, k @@ -414,14 +415,18 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch real(dp) :: delta_lnMbar, Ppoint, dlnP_dm, alfa real(dp) :: B_term, spatial_derivative_dX_dlnP, chiT real(dp) :: delta_lnP !, comp, y, t + logical :: patch_eos_partials + integer :: i - logical, parameter :: dbg = .false. + logical, parameter :: dbg = .true. include 'formats' ierr = 0 s% brunt_B(k) = 0 + patch_eos_partials = (s%eos_frac_PC(k)+s%eos_frac_Skye(k)+s%eos_frac_ideal(k) > 0d0) + if (k <= 1) return ! should we add a check for nz? logT_face = log10(T_face) @@ -441,6 +446,14 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch d_eos_dxa, ierr) if (ierr /= 0) return + !Bill/Adam's hack for eos partials from finite difference when none are provided. + !Can be removed when EOS fully provides composition partials + ! temp block, set to false currently as this needs to be inside the newton solver. + if (s% fix_d_eos_dxa_partials .and. .false.) then + call fix_d_eos_dxa_partials(s, k, ierr) + if (ierr /= 0) return + end if + ! Extract chiT chiT = d_eos_dlnT(i_lnPgas) ! should we be using chiT_face instead? @@ -460,9 +473,38 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! Compute the Brunt B composition term do i = 1, species spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP - B_term = B_term - d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP + if (s% fix_d_eos_dxa_partials .and. patch_eos_partials) then + ! use the finite-differences to estimate eos partials + B_term = B_term - s%dlnPeos_dxa_for_partials(i,k) * spatial_derivative_dX_dlnP + else + ! use eos derivatives (if provided), otherwise returns 0 + B_term = B_term - d_eos_dxa(i_lnPgas, i) * spatial_derivative_dX_dlnP + end if + +! ! Debugging block +! if (dbg .and. k >= 500 .and. k <= nz) then +! write(*,'(A,I5)') 'Diagnostic at zone k =', k +! write(*,'(A,1P,E12.5)') ' T_face = ', T_face +! write(*,'(A,1P,E12.5)') ' rho_face = ', rho_face +! write(*,'(A,1P,E12.5)') ' chiT = ', chiT +! +! ! He-4 +! spatial_derivative_dX_dlnP = (s%xa(ihe4,k-1) - s%xa(ihe4,k)) / delta_lnP +! write(*,'(A,1P,E12.5)') ' dlnP/dX_He4 = ', d_eos_dxa(i_lnPgas, ihe4) +! write(*,'(A,1P,E12.5)') ' dX_He4/dlnP = ', spatial_derivative_dX_dlnP +! write(*,'(A,1P,E12.5)') ' contrib_He4 = ', -d_eos_dxa(i_lnPgas, ihe4)*spatial_derivative_dX_dlnP +! +! ! C-12 +! spatial_derivative_dX_dlnP = (s%xa(ic12,k-1) - s%xa(ic12,k)) / delta_lnP +! write(*,'(A,1P,E12.5)') ' dlnP/dX_C12 = ', d_eos_dxa(i_lnPgas, ic12) +! write(*,'(A,1P,E12.5)') ' dX_C12/dlnP = ', spatial_derivative_dX_dlnP +! write(*,'(A,1P,E12.5)') ' contrib_C12 = ', -d_eos_dxa(i_lnPgas, ic12)*spatial_derivative_dX_dlnP +! +! write(*,*) +! end if end do + ! might be necessary to use a compensator for large sums... ! do i = 1, species ! spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP @@ -478,7 +520,7 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch ! add term accounting for the composition-related gradient in gravitational mass if (s% use_mass_corrections) then delta_lnMbar = log(s% mass_correction(k-1)) - log(s% mass_correction(k)) - s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP/ChiT + s% brunt_B(k) = s% brunt_B(k) - chiRho_face*delta_lnMbar/delta_lnP/chiT end if ! Check for bad numbers @@ -492,8 +534,206 @@ subroutine get_brunt_B_from_eos_partials(s, species, nz, k, T_face, rho_face, ch call mesa_error(__FILE__, __LINE__, 'get_brunt_B_from_eos_partials') end if end if - end subroutine get_brunt_B_from_eos_partials + contains + + subroutine fix_d_eos_dxa_partials(s, k, ierr) + ! revise composition partials + ! subroutine can be removed when EOS fully provides composition partials + + use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results, i_lnE, i_lnPgas + use eos_support, only: get_eos + use star_utils, only: lookup_nameofvar + type (star_info), pointer :: s + integer, intent(in) :: k + integer, intent(out) :: ierr + + integer :: j + + logical, parameter :: debug = .false. + + ! these vars are for faking composition derivatives + real(dp), dimension(num_eos_basic_results) :: & + res, dres_dlnd, dres_dlnT + real(dp) :: dres_dxa(num_eos_d_dxa_results, s% species) + real(dp) :: dxa + real(dp) :: xa_start_1(s% species) + real(dp) :: frac_without_dxa + real(dp) :: lnE_with_xa_start, lnPgas_with_xa_start + + integer :: i_var, i_var_sink + + real(dp), parameter :: dxa_threshold = 1d-4 + + logical, parameter :: checking = .true. + + include 'formats' + + ierr = 0 + + ! some EOSes have composition partials and some do not + ! those currently without dx partials are PC & Skye & ideal + frac_without_dxa = s% eos_frac_PC(k) + s% eos_frac_Skye(k) + s% eos_frac_ideal(k) + + if (debug .and. k == s% solver_test_partials_k) then + write(*,2) 's% eos_frac_PC(k)', k, s% eos_frac_PC(k) + write(*,2) 's% eos_frac_Skye(k)', k, s% eos_frac_Skye(k) + write(*,2) 's% eos_frac_ideal(k)', k, s% eos_frac_ideal(k) + write(*,2) 'frac_without_dxa', k, frac_without_dxa + end if + + if (k == s% solver_test_partials_k .and. s% solver_iter == s% solver_test_partials_iter_number) then + i_var = lookup_nameofvar(s, s% solver_test_partials_var_name) + if (i_var > s% nvar_hydro) then + i_var_sink = lookup_nameofvar(s, s% solver_test_partials_sink_name) + end if + end if + + ! if we're on an EOS where there aren't composition partials, + ! approximate derivatives with finite differences + if (frac_without_dxa > 0) then + s% dlnE_dxa_for_partials (:,k) = 0d0 + s% dlnPeos_dxa_for_partials(:,k) = 0d0 + + do j=1, s% species + !write (*,*) 'k', k + !write (*,*) 'j', j + !write (*,*) 's% xa_start', s% xa_start(j,k) + dxa = s% xa(j,k) - s% xa_start(j,k) + + if (debug .and. k == s% solver_test_partials_k .and. & + s% solver_iter == s% solver_test_partials_iter_number) & + write(*,2) 'dxa', j, dxa + + if (abs(dxa) >= dxa_threshold) then + + ! first, get eos with xa_start + + call get_eos( & + s, k, s% xa_start(:,k), & + rho_face, logRho_face, T_face, logT_face, & + res, dres_dlnd, dres_dlnT, dres_dxa, ierr) + if (ierr /= 0) then + if (s% report_ierr) write(*,2) 'failed in get_eos with xa_start', k + return + end if + + + if (is_bad(res(i_lnPgas)) .or. is_bad(res(i_lnE)) .or. is_bad(dres_dxa(i_lnPgas,j)) .or. is_bad(dres_dxa(i_lnE,j))) then + write(*,*) '--- BAD IN GET_EOS (perturbed) ---' + write(*,*) 'k =',k,' j =',j + write(*,*) 'rho_face =',rho_face,' T_face =',T_face + write(*,*) 'xa_start_1(j) =',xa_start_1(j),' dxa =',dxa + write(*,*) 'res(i_lnPgas) =',res(i_lnPgas) + write(*,*) 'res(i_lnE) =',res(i_lnE) + write(*,*) 'dres_dxa(lnP) =',dres_dxa(i_lnPgas,j) + write(*,*) 'dres_dxa(lnE) =',dres_dxa(i_lnE,j) + end if + + lnE_with_xa_start = res(i_lnE) + lnPgas_with_xa_start = res(i_lnPgas) + + ! now, get eos with 1 iso perturbed + + xa_start_1 = s% xa_start(:,k) + xa_start_1(j) = s% xa_start(j,k) + dxa + + call get_eos( & + s, k, xa_start_1, & + rho_face, logRho_face, T_face, logT_face, & + res, dres_dlnd, dres_dlnT, dres_dxa, ierr) + + + + if (is_bad(res(i_lnPgas)).or. is_bad(res(i_lnE))) ierr = -1 ! ierr /= 0 + if (ierr /= 0) then + + ! punt silently for now + s% dlnE_dxa_for_partials(:,k) = 0d0 + s% dlnPeos_dxa_for_partials(:,k) = 0d0 + ierr = 0 + return + + if (s% report_ierr) write(*,2) 'failed in get_eos with xa_start_1', k + return + end if + + + if (is_bad(dres_dxa(i_lnPgas,j))) dres_dxa(i_lnPgas,j) = 0d0 + if (is_bad(dres_dxa(i_lnE ,j))) dres_dxa(i_lnE ,j) = 0d0 + + ! fix up derivatives + + if (debug .and. k == s% solver_test_partials_k) & ! .and. s% solver_iter == s% solver_test_partials_iter_number) & + write(*,2) 'res(i_lnE) - lnE_with_xa_start', j, res(i_lnE) - lnE_with_xa_start + + s% dlnE_dxa_for_partials(j,k) = dres_dxa(i_lnE, j) + & + frac_without_dxa * (res(i_lnE) - lnE_with_xa_start) / dxa + if (checking) call check_dxa(j,k,s% dlnE_dxa_for_partials(j,k),'dlnE_dxa_for_partials') + + s% dlnPeos_dxa_for_partials(j,k) = s% Pgas(k)*dres_dxa(i_lnPgas,j)/s% Peos(k) + & + frac_without_dxa * (s% Pgas(k)/s% Peos(k)) * (res(i_lnPgas) - lnPgas_with_xa_start) / dxa + if (.false. .and. s% model_number == 1100 .and. k == 151 .and. & + j == 1 .and. is_bad(s% dlnPeos_dxa_for_partials(j,k))) then + write(*,2) 's% Pgas(k)', k, s% Pgas(k) + write(*,2) 'dres_dxa(i_lnPgas,j)', k, dres_dxa(i_lnPgas,j) + write(*,2) 's% Peos(k)', k, s% Peos(k) + write(*,2) 'frac_without_dxa', k, frac_without_dxa + write(*,2) 'res(i_lnPgas)', k, res(i_lnPgas) + write(*,2) 'lnPgas_with_xa_start', k, lnPgas_with_xa_start + write(*,2) 'dxa', k, dxa + write(*,2) 'dxa_threshold', k, dxa_threshold + write(*,2) 's% xa_start(j,k)', j, s% xa_start(j,k) + write(*,2) 'xa_start_1(j)', j, xa_start_1(j) + write(*,2) 's% eos_frac_PC(k)', k, s% eos_frac_PC(k) + write(*,2) 's% eos_frac_Skye(k)', k, s% eos_frac_Skye(k) + end if + if (checking) call check_dxa(j,k,s% dlnPeos_dxa_for_partials(j,k),'dlnPeos_dxa_for_partials') + + else + + if (k == s% solver_test_partials_k .and. s% solver_iter == s% solver_test_partials_iter_number) then + if (i_var_sink > 0 .and. i_var > s% nvar_hydro) then + if (dxa < dxa_threshold) then + if (j == i_var - s% nvar_hydro) then + write(*,*) 'fix_d_eos_dxa_partials: skipping dxa derivative fix for ', & + trim (s% solver_test_partials_var_name), & + ' (dxa < dxa_threshold): ', abs(dxa), ' < ', dxa_threshold + end if + if (j == i_var_sink - s% nvar_hydro) then + write(*,*) 'fix_d_eos_dxa_partials: skipping dxa derivative fix for ', & + trim (s% solver_test_partials_sink_name), & + ' (dxa < dxa_threshold): ', abs(dxa), ' < ', dxa_threshold + end if + end if + end if + end if + + end if + end do + + end if + + end subroutine fix_d_eos_dxa_partials + + subroutine check_dxa(j, k, dxa, str) + integer, intent(in) :: j, k + real(dp), intent(in) :: dxa + character (len=*), intent(in) :: str + include 'formats' + if (is_bad(dxa)) then +!$omp critical (get_brunt_B_from_eos_partials1) + ierr = -1 + if (s% report_ierr) then + write(*,3) 'get_brunt_B_from_eos_partials: bad ' // trim(str), j, k, dxa + end if + if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'get_brunt_B_from_eos_partials') +!$omp end critical (get_brunt_B_from_eos_partials1) + return + end if + end subroutine check_dxa + + end subroutine get_brunt_B_from_eos_partials end module brunt