Merge pull request #774 from MESAHub/rotation_update

Author: Debraheem

docs/source/changelog.rst

@@ -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```).
 

star/defaults/controls.defaults

@@ -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'
@@ -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
@@ -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
 
 
@@ -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|.
 
@@ -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)
 
       ! ::
@@ -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
@@ -8299,7 +8263,7 @@
     drag_coefficient = 0d0
     min_q_for_drag = 0d0
 
- 
+
       ! v_drag_factor
       ! ~~~~~~~~~~~~~
       ! v_drag
@@ -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.
 

star/private/ctrls_io.f90

@@ -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, &
@@ -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
@@ -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

star/private/hydro_eqns.f90

@@ -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'
@@ -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