diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 1bbe08907..097e22213 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -7737,6 +7737,15 @@ calculate_Brunt_B = .true. calculate_Brunt_N2 = .true. + ! calculate_Brunt_N2 + ! ~~~~~~~~~~~~~~~~~~ + + ! If true calculate B using equation 6 in MESA II, 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 c55214a9b..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 @@ -62,6 +62,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 @@ -332,4 +343,397 @@ 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) + 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), allocatable, dimension(:,:) :: xa_face! (species, k) + 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), xa_face(species, nz)) + + ! should we parallelize 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 + + 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_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 + !$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, xa_face, ierr) + 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 + 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(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) :: 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 = .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) + logRho_face = log10(rho_face) + Prad_face = crad * pow4(T_face) / 3d0 + + 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, xa_face(:,k), & + rho_face, logRho_face, T_face, logT_face, & + res, d_eos_dlnd, d_eos_dlnT, & + 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? + + ! 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) ! 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 + spatial_derivative_dX_dlnP = (s%xa(i,k-1) - s%xa(i,k)) / delta_lnP + 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 +! 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 + + ! 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 + 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 + + 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 diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 8487d4884..375acf6bc 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -146,8 +146,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, & @@ -1134,6 +1134,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 @@ -2825,6 +2826,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 e4bc1c8c4..95ac4a56c 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