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

w_div_w_crit_roche update #637

Merged
merged 7 commits into from
Jan 31, 2025
Merged

w_div_w_crit_roche update #637

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 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
cec2f77
Merge branch 'main' into rotation_update
Jan 30, 2025
22b47b2
changlog entry
Jan 30, 2025
9f1366c
[ci skip] wording in defaults
Jan 30, 2025
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
4 changes: 3 additions & 1 deletion docs/source/changelog.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,15 @@ New Features

A pseudo drag term ``v_drag`` has been reintroduced for ``u_flag`` to damp spurious shocks.

``hydro_rotation`` now contains the more accurate deformation fits from Fabry+2022, A&A 661, A123

.. _Bug Fixes main:

Bug Fixes
---------

fixed small bug in star/private/create_initial_model.f90 that will have a small effect on creating initial models

fixed bug in ``star/private/hydro_rotation.f90`` where the sigmoid function to cap ``w_div_w_crit`` was incorrectly implemented. This only influences models with `w_div_wc_flag = .true.`

.. note:: Before releasing a new version of MESA, move `Changes in main` to a new section below with the version number as the title, and add a new `Changes in main` section at the top of the file (see ```changelog_template.rst```).

Expand Down
74 changes: 19 additions & 55 deletions star/defaults/controls.defaults
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2313,7 +2313,7 @@
! convection model ``TDC`` with an overshooting model, or combining ``TDC``
! with other models for chemical composition gradients
! (predictive mixing or convective premixing), rotation, etc. (MESA VI)

! ::

MLT_option = 'TDC'
Expand Down Expand Up @@ -4060,57 +4060,21 @@

! ::

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.
! 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 = 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.
! we use two limiting values w_div_wcrit_max < w_div_wcrit_max2 to provide a smooth transition.
! Nothing is done when w_div_wc < w_div_wcrit_max, if w_div_wcrit_max < w_div_wc < w_div_wcrit_max2,
! we apply a sigmoid, and in the limit of j_rot -> infinity, the resulting w_div_wc will match
! w_div_wcrit_max2 (being on the top of the sigmoid)

! ::

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 Expand Up @@ -4142,17 +4106,17 @@
! For numerical stability, turn off rotational part of ``D_mix`` at very low tau.

! ::

D_mix_rotation_min_tau_full_off = 0d0

! D_mix_rotation_max_tau_full_on
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! Use rotational components of ``D_mix`` for locations where tau >= this.
! For numerical stability, turn off rotational part of ``D_mix`` at very low tau.

! ::

D_mix_rotation_min_tau_full_on = 0d0


Expand Down Expand Up @@ -4673,8 +4637,8 @@
! do_starspots
! ~~~~~~~~~~~~

! If ``.true.``, switch on impedance of the surface flux due to magnetic pressure from starspots,
! parameterized in the style of an atmospheric boundary modification. First described by
! If ``.true.``, switch on impedance of the surface flux due to magnetic pressure from starspots,
! parameterized in the style of an atmospheric boundary modification. First described by
! `Somers et al. (2015; ApJ) <https://ui.adsabs.harvard.edu/abs/2015ApJ...807..174S>`__.
! Detailed discussion of this functionality can be found in |MESA VI|.

Expand All @@ -4695,8 +4659,8 @@
! xspot
! ~~~~~

! Temperature contrast between the spotted and unspotted regions.
! Valid values are between 1.0 (no contribution from magnetic pressure) and 0.5
! Temperature contrast between the spotted and unspotted regions.
! Valid values are between 1.0 (no contribution from magnetic pressure) and 0.5
! (half of the total pressure is due to magnetic pressure)

! ::
Expand Down Expand Up @@ -8263,11 +8227,11 @@

use_dPrad_dm_form_of_T_gradient_eqn = .false.
use_gradT_actual_vs_gradT_MLT_for_T_gradient_eqn = .false.


! Hydrodynamic drag
! =================

! drag_coefficient
! ~~~~~~~~~~~~~~~~
! min_q_for_drag
Expand Down Expand Up @@ -8299,7 +8263,7 @@
drag_coefficient = 0d0
min_q_for_drag = 0d0


! v_drag_factor
! ~~~~~~~~~~~~~
! v_drag
Expand All @@ -8312,7 +8276,7 @@

! Only when u_flag = .true.. Adds a pseudo drag term of the form
! -v_drag_factor*(v-v_drag)^2/r, can be used damp velocities in outer layers
! of a star, useful for smoothing out spurious shocks colliding in ejected layers.
! of a star, useful for smoothing out spurious shocks colliding in ejected layers.
! Effect is full on for q>q_for_v_drag_full_on and full off for
! q < q_for_v_drag_full_off.

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, &
D_mix_rotation_min_tau_full_off, D_mix_rotation_min_tau_full_on, &
set_uniform_am_nu_non_rot, uniform_am_nu_non_rot, &
Expand Down Expand Up @@ -1679,10 +1678,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 @@ -3365,10 +3360,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 @@ -601,10 +601,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
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 @@ -618,34 +617,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