diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 23edc597e..da2d9421a 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -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_max2infinity, the resulting w_div_wc will match - ! w_div_wcrit_max, while nothing is done when w_div_wcrit_max 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 diff --git a/star/private/hydro_rotation.f90 b/star/private/hydro_rotation.f90 index 19e4e007e..a140355e2 100644 --- a/star/private/hydro_rotation.f90 +++ b/star/private/hydro_rotation.f90 @@ -25,65 +25,66 @@ ! Routine eval_fp_ft for computing rotation corrections to the stellar structure equations. -! Following Endal & Sofia, 1976, ApJ 210:184. - +! Following the method of Kippenhahn & Thomas, 1970; Endal & Sofia 1976, as implemented in +! Paxton et al., 2019, using the updated fits from Fabry, Marchant & Sana, 2022 module hydro_rotation - use const_def, only: pi, pi4, ln10, two_thirds, one_third + use const_def, only: pi, pi4, ln10, two_thirds, one_third, one_sixth use star_utils, only: get_r_from_xh use star_private_def implicit none + real(dp), parameter :: log_term_power = 5.626d0 + contains - ! compute w_div_w_roche for a known angular frequency omega, rphi, and Mphi - real(dp) function w_div_w_roche_omega(rphi,Mphi,omega,cgrav, max_w, max_w2, w_div_wc_flag) result(w_roche) - real(dp), intent(in) :: rphi,Mphi,omega,cgrav, max_w, max_w2 + ! compute w_div_w_roche for a known angular frequency omega, rpsi, and Mpsi + real(dp) function w_div_w_roche_omega(rpsi, Mphi, omega, cgrav, max_w, max_w2, w_div_wc_flag) result(w_roche) + real(dp), intent(in) :: rpsi, Mphi, omega, cgrav, max_w, max_w2 logical, intent(in) :: w_div_wc_flag - real(dp) :: wr, wr_high, wr_low, dimless_rphi, new_dimless_rphi, rphi_lim1, rphi_lim2 + real(dp) :: wr, wr_high, wr_low, dimless_rpsi, new_dimless_rpsi, rpsi_lim1, rpsi_lim2 if (omega == 0d0) then w_roche = 0d0 return end if - dimless_rphi = rphi*pow(abs(omega), two_thirds)/pow(cgrav*Mphi,one_third) + dimless_rpsi = rpsi * pow(abs(omega), two_thirds) / pow(cgrav * Mphi, one_third) if (.not. w_div_wc_flag) then ! verify if w_div_w_roche is not above max, otherwise limit it to that wr = max_w - new_dimless_rphi = pow(wr,two_thirds)*(1-pow2(wr)/6d0+0.01726d0*pow4(wr)-0.03569d0*pow6(wr)) - if (dimless_rphi > new_dimless_rphi) then + new_dimless_rpsi = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr)) + if (dimless_rpsi > new_dimless_rpsi) then w_roche = wr return end if else - ! smoothly cap to max_w to get a continuous function - ! nothing is done when we are below max_w2, but between max_w2 and max_w we smoothly - ! produce an asymptote that would result in w_div_wc=max_w for jrot->infinity + ! smoothly cap to max_w2 to get a continuous function + ! nothing is done when we are below max_w, but between max_w and max_w2 we smoothly + ! produce an asymptote that would result in w_div_wc=max_w2 for jrot->infinity wr = max_w - rphi_lim1 = pow(wr,two_thirds)*(1-pow2(wr)/6d0+0.01726d0*pow4(wr)-0.03569d0*pow6(wr)) + rpsi_lim1 = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr)) wr = max_w2 - rphi_lim2 = pow(wr,two_thirds)*(1-pow2(wr)/6d0+0.01726d0*pow4(wr)-0.03569d0*pow6(wr)) + rpsi_lim2 = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr)) - if (abs(dimless_rphi) > rphi_lim1) then - dimless_rphi = & - 2*(rphi_lim2-rphi_lim1)/(1+exp(-2*(abs(dimless_rphi)-rphi_lim1)/(rphi_lim2-rphi_lim1)))-rphi_lim2+2*rphi_lim1 + if (abs(dimless_rpsi) > rpsi_lim1) then + dimless_rpsi = sigmoid(abs(dimless_rpsi), rpsi_lim1, rpsi_lim2) end if end if ! otherwise, bisect result wr_high = wr wr_low = 0 - do while (wr_high-wr_low>1d-6) - wr = 0.5d0*(wr_high+wr_low) - new_dimless_rphi = pow(wr,two_thirds)*(1-pow2(wr)/6d0+0.01726d0*pow4(wr)-0.03569d0*pow6(wr)) - if (dimless_rphi > new_dimless_rphi) then + do while (wr_high - wr_low > 1d-6) + wr = 0.5d0 * (wr_high + wr_low) + new_dimless_rpsi = pow(wr, two_thirds) * rpsi_from_re_factor(pow2(wr), pow4(wr), pow6(wr)) + if (dimless_rpsi > new_dimless_rpsi) then wr_low = wr else wr_high = wr @@ -97,19 +98,19 @@ real(dp) function w_div_w_roche_omega(rphi,Mphi,omega,cgrav, max_w, max_w2, w_di end function w_div_w_roche_omega - ! compute w_div_w_roche for a known specific angular momentum jrot, rphi, and Mphi - real(dp) function w_div_w_roche_jrot(rphi,Mphi,jrot,cgrav, max_w, max_w2, w_div_wc_flag) result(w_roche) - real(dp), intent(in) :: rphi,Mphi,jrot,cgrav, max_w, max_w2 + ! compute w_div_w_roche for a known specific angular momentum jrot, rpsi, and Mphi + real(dp) function w_div_w_roche_jrot(rpsi, Mphi, jrot, cgrav, max_w, max_w2, w_div_wc_flag) result(w_roche) + real(dp), intent(in) :: rpsi,Mphi,jrot,cgrav, max_w, max_w2 logical, intent(in) :: w_div_wc_flag real(dp) :: wr, wr_high, wr_low, dimless_factor, new_dimless_factor - real(dp) :: w2, w4, w6, lg_one_sub_w4, jr_lim1, jr_lim2, A, C + real(dp) :: w2, w4, w6, w_log_term, jr_lim1, jr_lim2 if (jrot == 0d0) then w_roche = 0d0 return end if - dimless_factor = abs(jrot)/sqrt(cgrav*Mphi*rphi) + dimless_factor = abs(jrot) / sqrt(cgrav * Mphi * rpsi) if (.not. w_div_wc_flag) then ! verify if w_div_w_roche is not above max, otherwise limit it to that @@ -117,50 +118,50 @@ real(dp) function w_div_w_roche_jrot(rphi,Mphi,jrot,cgrav, max_w, max_w2, w_div_ w2 = pow2(wr) w4 = pow4(wr) w6 = pow6(wr) - lg_one_sub_w4 = log(1d0-w4) - new_dimless_factor = two_thirds*wr*(1+17d0/60d0*w2-0.3436d0*w4-0.4055d0*w6-0.9277d0*lg_one_sub_w4) & - /(1d0-0.1076d0*w4-0.2336d0*w6-0.5583d0*lg_one_sub_w4) + w_log_term = log(1d0 - pow(wr, log_term_power)) + new_dimless_factor = two_thirds * wr * C(w2, w4, w6, w_log_term) / A(w4, w6, w_log_term) if (dimless_factor > new_dimless_factor) then w_roche = wr return end if else - ! smoothly cap to max_w to get a continuous function - ! nothing is done when we are below max_w2, but between max_w2 and max_w we smoothly - ! produce an asymptote that would result in w_div_wc=max_w for jrot->infinity + ! smoothly cap to max_w2 to get a continuous function + ! nothing is done when we are below max_w, but between max_w and max_w2 we smoothly + ! produce an asymptote that would result in w_div_wc=max_w2 for jrot->infinity wr = max_w - A = 1d0-0.1076d0*pow4(wr)-0.2336d0*pow6(wr)-0.5583d0*log(1d0-pow4(wr)) - C = 1d0+17d0/60d0*pow2(wr)-0.3436d0*pow4(wr)-0.4055d0*pow6(wr)-0.9277d0*log(1d0-pow4(wr)) - jr_lim1 = two_thirds*wr*C/A + w4 = pow4(wr) + w6 = pow6(wr) + w_log_term = log(1 - pow(wr, log_term_power)) + jr_lim1 = two_thirds * wr * C(pow2(wr), w4, w6, w_log_term) / A(w4, w6, w_log_term) wr = max_w2 - A = 1d0-0.1076d0*pow4(wr)-0.2336d0*pow6(wr)-0.5583d0*log(1d0-pow4(wr)) - C = 1d0+17d0/60d0*pow2(wr)-0.3436d0*pow4(wr)-0.4055d0*pow6(wr)-0.9277d0*log(1d0-pow4(wr)) - jr_lim2 = two_thirds*wr*C/A + w4 = pow4(wr) + w6 = pow6(wr) + w_log_term = log(1 - pow(wr, log_term_power)) + jr_lim2 = two_thirds * wr * C(pow2(wr), w4, w6, w_log_term) / A(w4, w6, w_log_term) if (abs(dimless_factor) > jr_lim1) then - dimless_factor = 2*(jr_lim2-jr_lim1)/(1+exp(-2*(abs(dimless_factor)-jr_lim1)/(jr_lim2-jr_lim1)))-jr_lim2+2*jr_lim1 + dimless_factor = sigmoid(abs(dimless_factor), jr_lim1, jr_lim2) end if end if ! otherwise, bisect result wr_high = wr wr_low = 0 - do while (wr_high-wr_low>1d-6) - wr = 0.5d0*(wr_high+wr_low) + do while (wr_high - wr_low > 1d-6) + wr = 0.5d0 * (wr_high + wr_low) w2 = pow2(wr) w4 = pow4(wr) w6 = pow6(wr) - lg_one_sub_w4 = log(1d0-w4) - new_dimless_factor = two_thirds*wr*(1+17d0/60d0*w2-0.3436d0*w4-0.4055d0*w6-0.9277d0*lg_one_sub_w4) & - /(1d0-0.1076d0*w4-0.2336d0*w6-0.5583d0*lg_one_sub_w4) + w_log_term = log(1d0 - pow(wr, log_term_power)) + new_dimless_factor = two_thirds * wr * C(w2, w4, w6, w_log_term) / A(w4, w6, w_log_term) if (dimless_factor > new_dimless_factor) then wr_low = wr else wr_high = wr end if end do - w_roche = 0.5d0*(wr_high+wr_low) + w_roche = 0.5d0 * (wr_high + wr_low) if (jrot < 0d0) then w_roche = -w_roche @@ -168,51 +169,6 @@ real(dp) function w_div_w_roche_jrot(rphi,Mphi,jrot,cgrav, max_w, max_w2, w_div_ end function w_div_w_roche_jrot - subroutine eval_i_rot(s,k,r00,w_div_w_crit_roche, i_rot) - use auto_diff_support - type (star_info), pointer :: s - integer, intent(in) :: k ! just for debugging - real(dp), intent(in) :: r00,w_div_w_crit_roche - type(auto_diff_real_star_order1), intent(out) :: i_rot - - type(auto_diff_real_2var_order1) :: ir, r, re, w, w2, w4, w6, lg_one_sub_w4, B, A - - include 'formats' - - i_rot = 0d0 - if (s% use_other_eval_i_rot) then - call s% other_eval_i_rot(s% id,k,r00,w_div_w_crit_roche, i_rot) - else if (s% simple_i_rot_flag) then - i_rot = (2d0/3d0)*r00*r00 - i_rot% d1Array(i_lnR_00) = 2*i_rot% val - i_rot% d1Array(i_w_div_wc_00) = 0d0 - else - ! Compute i_rot following Paxton et al. 2019 (ApJs, 243, 10) - w = w_div_w_crit_roche - w%d1val1 = 1d0 - - r = r00 - r%d1val2 = r00 ! Makes the independent variable lnR - - w2 = pow2(w) - w4 = pow4(w) - w6 = pow6(w) - lg_one_sub_w4 = log(1d0-w4) - re = r*(1d0+w2/6d0-0.0002507d0*w4+0.06075d0*w6) - B = (1d0+w2/5d0-0.2735d0*w4-0.4327d0*w6-3d0/2d0*0.5583d0*lg_one_sub_w4) - A = (1d0-0.1076d0*w4-0.2336d0*w6-0.5583d0*lg_one_sub_w4) - - ir = two_thirds*pow2(re)*B/A - - i_rot = 0d0 - i_rot = ir%val - i_rot% d1Array(i_w_div_wc_00) = ir%d1val1 - i_rot% d1Array(i_lnR_00) = ir%d1val2 - end if - - end subroutine eval_i_rot - - subroutine set_i_rot(s, skip_w_div_w_crit_roche) type (star_info), pointer :: s logical, intent(in) :: skip_w_div_w_crit_roche @@ -257,7 +213,6 @@ subroutine set_j_rot(s) end do end subroutine set_j_rot - subroutine set_omega(s, str) type (star_info), pointer :: s character (len=*) :: str @@ -268,7 +223,6 @@ subroutine set_omega(s, str) end do end subroutine set_omega - subroutine check_omega(s, str) type (star_info), pointer :: s character (len=*) :: str @@ -288,7 +242,6 @@ subroutine check_omega(s, str) call mesa_error(__FILE__,__LINE__,'check_omega') end subroutine check_omega - subroutine update1_i_rot_from_xh(s, k) type (star_info), pointer :: s integer, intent(in) :: k @@ -331,7 +284,6 @@ subroutine use_xh_to_update_i_rot_and_j_rot(s) call set_j_rot(s) end subroutine use_xh_to_update_i_rot_and_j_rot - subroutine get_rotation_sigmas(s, nzlo, nzhi, dt, ierr) type (star_info), pointer :: s integer, intent(in) :: nzlo, nzhi @@ -370,7 +322,6 @@ subroutine get_rotation_sigmas(s, nzlo, nzhi, dt, ierr) end subroutine get_rotation_sigmas - subroutine get1_am_sig(s, nzlo, nzhi, am_nu, am_sig, dt, ierr) type (star_info), pointer :: s integer, intent(in) :: nzlo, nzhi @@ -426,7 +377,6 @@ subroutine get1_am_sig(s, nzlo, nzhi, am_nu, am_sig, dt, ierr) end subroutine get1_am_sig - subroutine set_uniform_omega(id, omega, ierr) use auto_diff_support integer, intent(in) :: id @@ -472,7 +422,6 @@ subroutine set_uniform_omega(id, omega, ierr) s% need_to_setvars = .true. end subroutine set_uniform_omega - subroutine set_rotation_info(s, skip_w_div_w_crit_roche, ierr) type (star_info), pointer :: s logical, intent(in) :: skip_w_div_w_crit_roche @@ -500,7 +449,6 @@ subroutine set_rotation_info(s, skip_w_div_w_crit_roche, ierr) end if end subroutine set_rotation_info - subroutine set_surf_avg_rotation_info(s) use star_utils, only: get_Lrad_div_Ledd type (star_info), pointer :: s @@ -654,7 +602,6 @@ subroutine set_surf_avg_rotation_info(s) end subroutine set_surf_avg_rotation_info - ! Input variables: ! N Number of meshpoints used by the model (arrays are this size) ! XM Mass coordinate [gram] @@ -665,8 +612,7 @@ end subroutine set_surf_avg_rotation_info ! Correction factor FT at each meshpoint ! Correction factor FP at each meshpoint ! r_polar, r_equatorial at each meshpoint - subroutine eval_fp_ft( & - id, nz, xm, r, rho, aw, ft, fp, r_polar, r_equatorial, report_ierr, ierr) + subroutine eval_fp_ft(id, nz, xm, r, rho, aw, ft, fp, r_polar, r_equatorial, report_ierr, ierr) use num_lib use auto_diff_support integer, intent(in) :: id @@ -682,8 +628,7 @@ subroutine eval_fp_ft( & logical :: dbg - type(auto_diff_real_1var_order1) :: A_omega,fp_numerator, ft_numerator, w, w2, w3, w4, w5, w6, lg_one_sub_w4, & - d_A_omega_dw, d_fp_numerator_dw, d_ft_numerator_dw, fp_temp, ft_temp + type (auto_diff_real_1var_order1) :: A_omega, w, w2, w4, w6, w_log_term, fp_temp, ft_temp include 'formats' @@ -693,48 +638,93 @@ subroutine eval_fp_ft( & dbg = .false. ! (s% model_number >= 5) -!$OMP PARALLEL DO PRIVATE(j, A_omega, fp_numerator, ft_numerator, fp_temp, ft_temp, w, w2, w3, w4, w5, w6, lg_one_sub_w4) SCHEDULE(dynamic,2) - do j=1, s% nz - !Compute fp, ft, re and rp using fits to the Roche geometry of a single star. - !by this point in the code, w_div_w_crit_roche is set - w = s% w_div_w_crit_roche(j) - w%d1val1 = 1d0 - - w2 = pow2(w) - w4 = pow4(w) - w6 = pow6(w) - lg_one_sub_w4 = log(1d0-w4) - A_omega = (1d0-0.1076d0*w4-0.2336d0*w6-0.5583d0*lg_one_sub_w4) - fp_numerator = (1d0-two_thirds*w2-0.06837d0*w4-0.2495d0*w6) - ft_numerator = (1d0+0.2185d0*w4-0.1109d0*w6) - !fits for fp, ft - fp_temp = fp_numerator/A_omega - ft_temp = ft_numerator/A_omega - !re and rp can be derived analytically from w_div_wcrit - r_equatorial(j) = r(j)*(1d0+w2% val/6d0-0.0002507d0*w4% val+0.06075d0*w6% val) - r_polar(j) = r_equatorial(j)/(1d0+0.5d0*w2% val) - ! Be sure they are consistent with r_Phi - r_equatorial(j) = max(r_equatorial(j),r(j)) - r_polar(j) = min(r_polar(j),r(j)) - - fp(j) = 0d0 - ft(j) = 0d0 - fp(j)% val = fp_temp% val - ft(j)% val = ft_temp% val - if (s% w_div_wc_flag) then - fp(j)% d1Array(i_w_div_wc_00) = fp_temp% d1val1 - ft(j)% d1Array(i_w_div_wc_00) = ft_temp% d1val1 - end if - !if (j == s% solver_test_partials_k) then - ! s% solver_test_partials_val = fp(j) - ! s% solver_test_partials_var = s% i_w_div_wc - ! s% solver_test_partials_dval_dx = s% dfp_rot_dw_div_wc(j) - !end if - end do +!$OMP PARALLEL DO PRIVATE(j, A_omega, fp_temp, ft_temp, w, w2, w4, w6, w_log_term) SCHEDULE(dynamic,2) + do j=1, s% nz + ! Compute fp, ft, re and rp using fits to the Roche geometry of a single star. + ! by this point in the code, w_div_w_crit_roche is set + w = abs(s% w_div_w_crit_roche(j)) + w% d1val1 = 1d0 + + w2 = pow2(w) + w4 = pow4(w) + w6 = pow6(w) + w_log_term = log(1d0 - pow(w, log_term_power)) + ! cannot use real function below because need derivatives + A_omega = 1d0 + 0.3293d0 * w4 - 0.4926d0 * w6 - 0.5560d0 * w_log_term + + ! fits for fp, ft; Fabry+2022, Eqs. A.10, A.11 + fp_temp = (1d0 - two_thirds * w2 - 0.2133d0 * w4 - 0.1068d0 * w6) / A_omega + ft_temp = (1d0 - 0.07955d0 * w4 - 0.2322d0 * w6) / A_omega + + ! re and rp can be derived analytically from w_div_wcrit + r_equatorial(j) = r(j) * re_from_rpsi_factor(w2% val, w4% val, w6% val) + r_polar(j) = r_equatorial(j) / (1d0 + 0.5d0 * w2% val) + + ! Be sure they are consistent with r_Psi + r_equatorial(j) = max(r_equatorial(j), r(j)) + r_polar(j) = min(r_polar(j), r(j)) + + fp(j) = 0d0 + ft(j) = 0d0 + fp(j)% val = fp_temp% val + ft(j)% val = ft_temp% val + if (s% w_div_wc_flag) then + fp(j)% d1Array(i_w_div_wc_00) = fp_temp% d1val1 + ft(j)% d1Array(i_w_div_wc_00) = ft_temp% d1val1 + end if + !if (j == s% solver_test_partials_k) then + ! s% solver_test_partials_val = fp(j) + ! s% solver_test_partials_var = s% i_w_div_wc + ! s% solver_test_partials_dval_dx = s% dfp_rot_dw_div_wc(j) + !end if + end do !$OMP END PARALLEL DO end subroutine eval_fp_ft + subroutine eval_i_rot(s, k, r00, w_div_w_crit_roche, i_rot) + use auto_diff_support + type (star_info), pointer :: s + integer, intent(in) :: k ! just for debugging + real(dp), intent(in) :: r00, w_div_w_crit_roche + type(auto_diff_real_star_order1), intent(out) :: i_rot + + type(auto_diff_real_2var_order1) :: ir, r, re, w, w2, w4, w6, w_log_term + + include 'formats' + + i_rot = 0d0 + if (s% use_other_eval_i_rot) then + call s% other_eval_i_rot(s% id, k, r00, w_div_w_crit_roche, i_rot) + else if (s% simple_i_rot_flag) then + i_rot = two_thirds * r00 * r00 + i_rot% d1Array(i_lnR_00) = 2d0 * i_rot% val + i_rot% d1Array(i_w_div_wc_00) = 0d0 + else + w = abs(w_div_w_crit_roche) + w% d1val1 = 1d0 + + r = r00 + r% d1val2 = r00 ! Makes the independent variable lnR + + w2 = pow2(w) + w4 = pow4(w) + w6 = pow6(w) + w_log_term = log(1d0 - pow(w, log_term_power)) + ! cannot use real functions below, since need derivatives + re = r * (1d0 + one_sixth * w2 - 0.005124d0 * w4 + 0.06562d0 * w6) + ! Compute i_rot following Fabry+2022, Eq. A.9 + ir = two_thirds * pow2(re) * (1d0 + w2 / 5d0 + 0.4140d0 * w4 - 0.8650d0 * w6 - 0.8370d0 * w_log_term) & + / (1d0 + 0.3293d0 * w4 - 0.4926d0 * w6 - 0.5560d0 * w_log_term) + + i_rot = 0d0 + i_rot = ir% val + i_rot% d1Array(i_w_div_wc_00) = ir% d1val1 + i_rot% d1Array(i_lnR_00) = ir% d1val2 + end if + + end subroutine eval_i_rot + subroutine compute_j_fluxes_and_extra_jdot(id, ierr) use auto_diff_support integer, intent(in) :: id @@ -801,6 +791,53 @@ subroutine compute_j_fluxes_and_extra_jdot(id, ierr) s% j_flux(s% nz) = 0d0 end subroutine compute_j_fluxes_and_extra_jdot + real(dp) function rpsi_from_re_factor(o2, o4, o6) + real(dp), intent(in) :: o2 + real(dp), intent(in) :: o4 + real(dp), intent(in) :: o6 + ! Fabry+2022, Eq. A.2 + rpsi_from_re_factor = 1d0 - one_sixth * o2 + 0.02025d0 * o4 - 0.03870d0 * o6 + end function rpsi_from_re_factor + + real(dp) function re_from_rpsi_factor(o2, o4, o6) + real(dp), intent(in) :: o2 + real(dp), intent(in) :: o4 + real(dp), intent(in) :: o6 + ! Fabry+2022, Eq. A.3 + re_from_rpsi_factor = 1d0 + one_sixth * o2 - 0.005124d0 * o4 + 0.06562d0 * o6 + end function re_from_rpsi_factor + + real(dp) function A(o4, o6, o_log_term) + real(dp), intent(in) :: o4 + real(dp), intent(in) :: o6 + real(dp), intent(in) :: o_log_term + ! Fabry+2022, Eq. A.6 + A = 1d0 + 0.3293d0 * o4 - 0.4926d0 * o6 - 0.5560d0 * o_log_term + end function A + + real(dp) function B(o2, o4, o6, o_log_term) + real(dp), intent(in) :: o2 + real(dp), intent(in) :: o4 + real(dp), intent(in) :: o6 + real(dp), intent(in) :: o_log_term + ! Fabry+2022, Eq. A.7 + B = 1d0 + o2 / 5d0 + 0.4140d0 * o4 - 0.8650d0 * o6 - 0.8370d0 * o_log_term + end function B + + real(dp) function C(o2, o4, o6, o_log_term) + real(dp), intent(in) :: o2 + real(dp), intent(in) :: o4 + real(dp), intent(in) :: o6 + real(dp), intent(in) :: o_log_term + ! Fabry+2022, Eq. A.12 + C = 1d0 + 17d0 / 60d0 * o2 + 0.4010d0 * o4 - 0.8606d0 * o6 - 0.9305d0 * o_log_term + end function C + + real(dp) function sigmoid(x, xmax1, xmax2) + real(dp), intent(in) :: x, xmax1, xmax2 + ! smoothly maps abs(x) = [xmax1, infty] to sigmoid = [xmax1, xmax2] + sigmoid = 2d0 * (xmax2 - xmax1) / (1d0 + exp(-2d0 * (abs(x) - xmax1) / (xmax2 - xmax1))) - xmax2 + 2d0 * xmax1 + end function sigmoid end module hydro_rotation diff --git a/star/public/star_lib.f90 b/star/public/star_lib.f90 index 101a77959..5a2331107 100644 --- a/star/public/star_lib.f90 +++ b/star/public/star_lib.f90 @@ -895,7 +895,7 @@ end subroutine star_set_u_flag ! simply adds or removes; doesn't reconverge the model. subroutine star_set_rotation_flag(id, rotation_flag, ierr) use set_flags, only: set_rotation_flag - use hydro_rotation, only: set_rotation_info, set_i_rot + use hydro_rotation, only: set_rotation_info integer, intent(in) :: id logical, intent(in) :: rotation_flag integer, intent(out) :: ierr diff --git a/star/test_suite/magnetic_braking/inlist_braking b/star/test_suite/magnetic_braking/inlist_braking index a4fa5f420..6698da965 100644 --- a/star/test_suite/magnetic_braking/inlist_braking +++ b/star/test_suite/magnetic_braking/inlist_braking @@ -38,7 +38,7 @@ &kap - use_type2_opacities = .true + use_type2_opacities = .true. Zbase = 1d-2 / ! end of kap namelist @@ -118,216 +118,3 @@ / ! end of controls namelist - - -&pgstar - !pause = .true. - !pgstar_show_log_age_in_years = .true. - pgstar_show_age_in_years = .true. - ! if true, the code waits for user to enter a RETURN on the command line - - file_white_on_black_flag = .true. - win_white_on_black_flag = .true. - - Grid2_win_flag = .true. - Grid2_num_cols = 10 ! divide plotting region into this many equal width cols - Grid2_num_plots =7 - Grid2_num_rows =8 - Grid2_win_width = 16 - !Grid2_win_aspect_ratio = 0.5 ! aspect_ratio = height/width - Grid2_xleft = 0.01 ! fraction of full window width for margin on left - Grid2_xright = 0.99 ! fraction of full window width for margin on right - Grid2_ybot = 0.12 ! fraction of full window width for margin on bottom - Grid2_ytop = 0.95 ! fraction of full window width for margin on top - show_TRho_Profile_eos_regions = .true. - - !Grid2_file_flag = .true. - Grid2_file_dir = 'png' - Grid2_file_prefix = 'Grid2' - Grid2_file_interval = 1 - - Grid2_file_width = 18 - - !Grid2_file_aspect_ratio = 0.7 - - - - Grid2_plot_name(1) = 'TRho_Profile' - Grid2_plot_row(1) = 1 - Grid2_plot_rowspan(1) =5 - Grid2_plot_col(1) = 1 - Grid2_plot_colspan(1) = 4 - Grid2_plot_pad_left(1) = 0.02 - Grid2_plot_pad_right(1) = 0.0 - Grid1_plot_pad_top(1) = 0.0 - Grid2_plot_pad_bot(1) = 0.15 - Grid2_txt_scale_factor(1) = 0.5 - - Grid2_plot_name(2) = 'HR' - Grid2_plot_row(2) = 5 - Grid2_plot_rowspan(2) = 3 - Grid2_plot_col(2) = 1 - Grid2_plot_colspan(2) = 2 - - Grid2_plot_pad_left(2) = 0.04 - Grid2_plot_pad_right(2) = 0.02 - Grid2_plot_pad_top(2) = 0.04 - Grid2_plot_pad_bot(2) = 0.00 - Grid2_txt_scale_factor(2) = 0.5 - - - Grid2_plot_name(3) = 'Abundance' - Grid2_plot_row(3) = 1 - Grid2_plot_rowspan(3) = 4 - Grid2_plot_col(3) = 5 - Grid2_plot_colspan(3) = 3 - - - Grid2_plot_pad_left(3) = 0.09 - Grid2_plot_pad_right(3) = 0.00 - Grid2_plot_pad_top(3) = 0.00 - Grid2_plot_pad_bot(3) = 0.05 - Grid2_txt_scale_factor(3) = 0.5 - Abundance_legend_max_cnt = 0 - Abundance_legend_txt_scale_factor =0.5 - - - show_TRho_annotation1 = .true. - show_TRho_annotation2 = .true. - show_TRho_annotation3 = .true. - show_TRho_degeneracy_line = .true. - - - Grid2_plot_name(4) = 'Text_Summary1' - Grid2_plot_row(4) = 7 - Grid2_plot_rowspan(4) = 2 - Grid2_plot_col(4) = 1 - Grid2_plot_colspan(4) = 7 - Grid2_plot_pad_left(4) = 0.0 - Grid2_plot_pad_right(4) = 0.0 - Grid2_plot_pad_top(4) = 0.2 - Grid2_plot_pad_bot(4) = -0.1 - Grid2_txt_scale_factor(4) = 0.18 - Text_Summary1_num_rows = 4 ! <= 20 - Text_Summary1_num_cols = 4 ! <= 20 - Text_Summary1_name(:,:) = '' - - Text_Summary1_name(1,1) = 'time_step' - Text_Summary1_name(1,2) = 'surf_avg_v_rot' - Text_Summary1_name(1,3) = 'surf_avg_v_div_v_crit' - Text_Summary1_name(1,4) = 'star_mass' - - Text_Summary1_name(2,1) = 'num_zones' - Text_Summary1_name(2,2) = 'star_mdot' - Text_Summary1_name(2,3) = 'star_age' - Text_Summary1_name(2,4) = 'model_number' - - Text_Summary1_name(3,1) = 'log_total_angular_momentum' - - - - - - - - - Grid2_plot_name(5) = 'History_Panels1' - Grid2_plot_row(5) = 1 - Grid2_plot_rowspan(5) = 9 - Grid2_plot_col(5) = 8 - Grid2_plot_colspan(5) = 3 - Grid2_plot_pad_left(5) = 0.075 - Grid2_plot_pad_right(5) = 0.04 - Grid2_plot_pad_top(5) = 0.0 - Grid2_plot_pad_bot(5) = 0.05 - Grid2_txt_scale_factor(5) = 0.5 - - History_Panels1_num_panels = 4 - - History_Panels1_xaxis_name = 'model_number' - History_Panels1_yaxis_name(1) = 'log_total_angular_momentum' - History_Panels1_yaxis_reversed(1) = .false. - History_Panels1_ymin(1) = 48d0 - History_Panels1_ymax(1) = -101d0 - History_Panels1_max_width = 500 - History_Panels1_dymin(1) = -1 - History_Panels1_other_yaxis_name(1) = 'log_R' - History_Panels1_other_yaxis_reversed(1) = .false. - History_Panels1_other_ymin(1) = -101d0 - History_Panels1_other_ymax(1) = -101d0 - History_Panels1_other_dymin(1) = -1 - - - - History_Panels1_yaxis_name(2) = 'surf_avg_v_rot' - History_Panels1_yaxis_reversed(2) = .false. - History_Panels1_ymin(2) = -101d0 - History_Panels1_ymax(2) = -101d0 - History_Panels1_dymin(2) = -1 - History_Panels1_other_yaxis_name(2) = 'surf_avg_v_div_v_crit' - History_Panels1_other_yaxis_reversed(2) = .false. - History_Panels1_other_ymin(2) = -101d0 - History_Panels1_other_ymax(2) = -101d0 - History_Panels1_other_dymin(2) = -1 - - - History_Panels1_yaxis_name(3) = 'star_mdot' - History_Panels1_yaxis_reversed(3) = .false. - History_Panels1_ymin(3) = -101d0 - History_Panels1_ymax(3) = -101d0 - History_Panels1_dymin(3) = -1 - History_Panels1_other_yaxis_name(3) = '' - History_Panels1_other_yaxis_reversed(3) = .false. - History_Panels1_other_ymin(3) = -101d0 - History_Panels1_other_ymax(3) = -101d0 - History_Panels1_other_dymin(3) = -1 - - History_Panels1_yaxis_name(4) = 'log_L' - History_Panels1_yaxis_reversed(4) = .false. - History_Panels1_ymin(4) = -101d0 - History_Panels1_ymax(4) = -101d0 - History_Panels1_dymin(4) = -1 - History_Panels1_other_yaxis_name(4) = 'log_L_div_Ledd' ! These are for profiles -> 'conv_L_div_L' ! 'log_Lrad_div_Ledd' ! 'log_L_div_Ledd' - History_Panels1_other_yaxis_reversed(4) = .false. - History_Panels1_other_ymin(4) = -101d0 - History_Panels1_other_ymax(4) = -101d0 - History_Panels1_other_dymin(4) = -1 - - - - Grid2_plot_name(6) = 'Kipp' - Grid2_plot_row(6) = 5 - Grid2_plot_rowspan(6) = 5 - Grid2_plot_col(6) = 5 - Grid2_plot_colspan(6) = 4 - Grid2_plot_pad_left(6) = 0.1 - Grid2_plot_pad_right(6) = 0.1 - Grid1_plot_pad_top(6) = 0.2 - Grid2_plot_pad_bot(6) = 0.2 - Grid2_txt_scale_factor(6) = 0.5 - - - - - - - Grid2_plot_name(7) = 'Profile_Panels1' - Grid2_plot_row(7) =5 - Grid2_plot_rowspan(7) =3 - Grid2_plot_col(7) =3 - Grid2_plot_colspan(7) = 2 - Profile_Panels1_num_panels = 1 - Profile_Panels1_xaxis_name = 'mass' - Profile_Panels1_yaxis_name(1) = 'omega' - !Profile_Panels1_other_yaxis_name(1) = 'engulfment_heating' - Profile_Panels1_other_ymin(:) = -101 ! only used if /= -101d0 - Profile_Panels1_other_ymax(:) = -101 ! only used if /= -101d0 - - - Grid2_plot_pad_left(7) = 0.04 - Grid2_plot_pad_right(7) = 0.02 - Grid2_plot_pad_top(7) = 0.04 - Grid2_plot_pad_bot(7) = 0.00 - Grid2_txt_scale_factor(7) = 0.4 - -/ ! end of pgstar namelist