Skip to content
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

Jump to bottom

w_div_w_crit_roche update #637

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

w_div_w_crit_roche update #637

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
02a9cf2
remove obsolete controls, reinterpret sigmoid
Mar 22, 2024
6f603b7
some clean up and use recent fp/ft fits
Mar 22, 2024
fe9ff90
use abs cause non-even power used down the line
Mar 27, 2024
92a991d
[ci skip] some clean-up in mag braking test case
Mar 27, 2024
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: 4 additions & 41 deletions star/defaults/controls.defaults
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4048,57 +4048,20 @@

! ::

w_div_wcrit_max = 0.90d0
w_div_wcrit_max = 0.89d0


! w_div_wcrit_max2
! ~~~~~~~~~~~~~~~~

! When w_div_wc_flag is true, rather than a hard limit on w_div_wcrit
! we use w_div_wcrit_max2<w_div_wcrit_max to provide a smooth transition.
! we use w_div_wcrit_max<w_div_wcrit_max2 to provide a smooth transition.
! In the limit of j_rot->infinity, the resulting w_div_wc will match
! w_div_wcrit_max, while nothing is done when w_div_wcrit_max<w_div_wcrit_max2
! w_div_wcrit_max2, while nothing is done when w_div_wcrit_max<w_div_wcrit_max

! ::

w_div_wcrit_max2 = 0.89d0


! FP_min
! ~~~~~~
! FT_min
! ~~~~~~

! Lower limits for rotational distortion corrections factors FP and FT.
! Used for the calculation when fitted_fp_ft_i_rot = .false., otherwise the
! limits are set using w_div_wcrit_max

! ::

FP_min = 0.75d0
FT_min = 0.95d0


! FP_error_limit
! ~~~~~~~~~~~~~~

! If calculate an fp < this, treat it as an error.
! Used for the calculation when fitted_fp_ft_i_rot = .false.

! ::

FP_error_limit = 0d0


! FT_error_limit
! ~~~~~~~~~~~~~~

! If calculate an ft < this, treat it as an error.
! Used for the calculation when fitted_fp_ft_i_rot = .false.

! ::

FT_error_limit = 0d0
w_div_wcrit_max2 = 0.90d0


! D_mix_rotation_max_logT_full_on
Expand Down
9 changes: 0 additions & 9 deletions star/private/ctrls_io.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,6 @@ module ctrls_io
min_q_for_adjust_J_lost, min_J_div_delta_J, max_mdot_redo_cnt, mdot_revise_factor, &
implicit_mdot_boost, min_years_dt_for_redo_mdot, surf_omega_div_omega_crit_limit, surf_omega_div_omega_crit_tol, &
w_div_wcrit_max, w_div_wcrit_max2, &
fp_min, ft_min, fp_error_limit, ft_error_limit, &
D_mix_rotation_max_logT_full_on, D_mix_rotation_min_logT_full_off, &
set_uniform_am_nu_non_rot, uniform_am_nu_non_rot, &
set_min_am_nu_non_rot, min_am_nu_non_rot, min_center_Ye_for_min_am_nu_non_rot, &
Expand Down Expand Up @@ -1670,10 +1669,6 @@ subroutine store_controls(s, ierr)
s% surf_omega_div_omega_crit_tol = surf_omega_div_omega_crit_tol
s% w_div_wcrit_max = w_div_wcrit_max
s% w_div_wcrit_max2 = w_div_wcrit_max2
s% fp_min = fp_min
s% ft_min = ft_min
s% fp_error_limit = fp_error_limit
s% ft_error_limit = ft_error_limit

s% D_mix_rotation_max_logT_full_on = D_mix_rotation_max_logT_full_on
s% D_mix_rotation_min_logT_full_off = D_mix_rotation_min_logT_full_off
Expand Down Expand Up @@ -3344,10 +3339,6 @@ subroutine set_controls_for_writing(s, ierr)
surf_omega_div_omega_crit_tol = s% surf_omega_div_omega_crit_tol
w_div_wcrit_max = s% w_div_wcrit_max
w_div_wcrit_max2 = s% w_div_wcrit_max2
fp_min = s% fp_min
ft_min = s% ft_min
fp_error_limit = s% fp_error_limit
ft_error_limit = s% ft_error_limit

D_mix_rotation_max_logT_full_on = s% D_mix_rotation_max_logT_full_on
D_mix_rotation_min_logT_full_off = s% D_mix_rotation_min_logT_full_off
Expand Down
39 changes: 23 additions & 16 deletions star/private/hydro_eqns.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -602,10 +602,9 @@ subroutine do1_w_div_wc_eqn(s, k, nvar, ierr)
integer, intent(in) :: k, nvar
integer, intent(out) :: ierr
integer :: i_equ_w_div_wc, i_w_div_wc
real(dp) :: wwc, dimless_rphi, dimless_rphi_given_wwc, w1, w2
real(dp) :: jr_lim1, jr_lim2, A, C
real(dp) :: wwc, wwc4, wwc6, wwc_log_term, dimless_rphi, dimless_rphi_given_wwc, w1, w2, jr_lim1, jr_lim2
type(auto_diff_real_star_order1) :: &
w_d_wc00, r00, jrot00, resid_ad, A_ad, C_ad, &
w_d_wc00, w4_d_wc00, w6_d_wc00, r00, w_log_term_d_wc00, jrot00, resid_ad, A_ad, C_ad, &
jrot_ratio, sigmoid_jrot_ratio
logical :: test_partials
include 'formats'
Expand All @@ -619,34 +618,42 @@ subroutine do1_w_div_wc_eqn(s, k, nvar, ierr)
i_w_div_wc = s% i_w_div_wc

wwc = s% w_div_wcrit_max
A = 1d0-0.1076d0*pow4(wwc)-0.2336d0*pow6(wwc)-0.5583d0*log(1d0-pow4(wwc))
C = 1d0+17d0/60d0*pow2(wwc)-0.3436d0*pow4(wwc)-0.4055d0*pow6(wwc)-0.9277d0*log(1d0-pow4(wwc))
jr_lim1 = two_thirds*wwc*C/A
wwc4 = pow4(wwc)
wwc6 = pow6(wwc)
wwc_log_term = log(1d0 - pow(wwc, log_term_power))
jr_lim1 = two_thirds * wwc * C(pow2(wwc), wwc4, wwc6, wwc_log_term) / A(wwc4, wwc6, wwc_log_term)

wwc = s% w_div_wcrit_max2
A = 1d0-0.1076d0*pow4(wwc)-0.2336d0*pow6(wwc)-0.5583d0*log(1d0-pow4(wwc))
C = 1d0+17d0/60d0*pow2(wwc)-0.3436d0*pow4(wwc)-0.4055d0*pow6(wwc)-0.9277d0*log(1d0-pow4(wwc))
jr_lim2 = two_thirds*wwc*C/A
wwc4 = pow4(wwc)
wwc6 = pow6(wwc)
wwc_log_term = log(1d0 - pow(wwc, log_term_power))
jr_lim2 = two_thirds * wwc * C(pow2(wwc), wwc4, wwc6, wwc_log_term) / A(wwc4, wwc6, wwc_log_term)

w_d_wc00 = wrap_w_div_wc_00(s, k)
A_ad = 1d0-0.1076d0*pow4(w_d_wc00)-0.2336d0*pow6(w_d_wc00)-0.5583d0*log(1d0-pow4(w_d_wc00))
C_ad = 1d0+17d0/60d0*pow2(w_d_wc00)-0.3436d0*pow4(w_d_wc00)-0.4055d0*pow6(w_d_wc00)-0.9277d0*log(1d0-pow4(w_d_wc00))
w4_d_wc00 = pow4(w_d_wc00)
w6_d_wc00 = pow6(w_d_wc00)
w_log_term_d_wc00 = log(1d0 - pow(w_d_wc00, log_term_power))
A_ad = 1d0 + 0.3293d0 * w4_d_wc00 - 0.4926d0 * w6_d_wc00 - 0.5560d0 * w_log_term_d_wc00
C_ad = 1d0 + 17d0 / 60d0 * pow2(w_d_wc00) + 0.4010d0 * w4_d_wc00 - 0.8606d0 * w6_d_wc00 &
- 0.9305d0 * w_log_term_d_wc00

r00 = wrap_r_00(s, k)
if (s% j_rot_flag) then
jrot00 = wrap_jrot_00(s, k)
jrot_ratio = jrot00/sqrt(s% cgrav(k)*s% m(k)*r00)
jrot_ratio = jrot00 / sqrt(s% cgrav(k) * s% m(k) * r00)
else
jrot_ratio = s% j_rot(k)/sqrt(s% cgrav(k)*s% m(k)*r00)
jrot_ratio = s% j_rot(k) / sqrt(s% cgrav(k) * s% m(k) * r00)
end if
if (abs(jrot_ratio% val) > jr_lim1) then
sigmoid_jrot_ratio = 2*(jr_lim2-jr_lim1)/(1+exp(-2*(abs(jrot_ratio)-jr_lim1)/(jr_lim2-jr_lim1)))-jr_lim2+2*jr_lim1
sigmoid_jrot_ratio = 2d0 * (jr_lim2-jr_lim1) &
/ (1d0 + exp(-2d0 * (abs(jrot_ratio) - jr_lim1) / (jr_lim2 - jr_lim1))) &
- jr_lim2 + 2 * jr_lim1
if (jrot_ratio% val < 0d0) then
sigmoid_jrot_ratio = -sigmoid_jrot_ratio
end if
resid_ad = (sigmoid_jrot_ratio - two_thirds*w_d_wc00*C_ad/A_ad)
resid_ad = (sigmoid_jrot_ratio - two_thirds * w_d_wc00 * C_ad / A_ad)
else
resid_ad = (jrot_ratio - two_thirds*w_d_wc00*C_ad/A_ad)
resid_ad = (jrot_ratio - two_thirds * w_d_wc00 * C_ad / A_ad)
end if

s% equ(i_equ_w_div_wc, k) = resid_ad% val
Expand Down
Loading