starspots Part 2 (#713)

* [ci skip] adding some stuff, not yet done * [ci skip] updates to starspots * [ci skip] run_stars_extras * [ci skip] fix controls formatting * [ci skip] remove redundant library * [ci skip] modify turb_support * [ci skip] switch run_stars_extra to plain again * [ci skip] move need_atm_Psurf to right place * [ci skip] mv starspot stuff in hydro_vars.f90 * oh I think I found it --------- Co-authored-by: Meridith Joyce <[email protected]>
MESAHub · Aug 9, 2024 · ae4228d · ae4228d
1 parent b54f0e5
commit ae4228d
Show file tree

Hide file tree

Showing 10 changed files with 199 additions and 288 deletions.
diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults
@@ -4640,20 +4640,34 @@
 ! =========
 
 
-      ! fspot
-      ! ~~~~~~~~~~~
+      ! do_starspots
+      ! ~~~~~~~~~~~~
+
+      ! If ``.true.``, switch on impedence 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|.
+
+      ! ::
 
-      ! Filling Factor of starspots (in [0,1])
+    do_starspots = .false.
 
+
+      ! fspot
+      ! ~~~~~
+
+      ! Filling Factor of starspots. Valid values between 0.0 and 1.0 (no spots to 100% coverage)
       ! ::
 
     fspot = 0d0
 
 
       ! xspot
-      ! ~~~~~~~~~~~
+      ! ~~~~~
 
-      ! Temperature contrast between the spotted and unspotted regions (in [0,1])
+      ! 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)
 
       ! ::
 
@@ -8291,7 +8305,7 @@
 
       ! Note that these parameters are not exactly the same
       ! as used by Paul Duffell.
-      ! His calibrated D is 2, where mesa has default D = 3 (see mesaIV paper).
+      ! His calibrated D is 2, where mesa has default D = 3 (see |MESA IV|).
       ! Users should try various values since the choice is not clear cut.
 
       ! ::

diff --git a/star/make/makefile_base b/star/make/makefile_base
@@ -205,6 +205,7 @@ SRCS = \
    star_solver.f90 \
    struct_burn_mix.f90 \
    winds.f90 \
+   starspots.f90 \
    gravity_darkening.f90 \
    mass_utils.f90 \
    eps_mdot.f90 \

diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90
@@ -492,7 +492,7 @@ module ctrls_io
     use_T_tau_gradr_factor, &
 
     ! starspots
-    fspot, xspot, &
+    do_starspots, fspot, xspot, &
 
     ! extra heat near surface to model irradiation
     irradiation_flux, column_depth_for_irradiation, &
@@ -1289,6 +1289,7 @@ subroutine store_controls(s, ierr)
  s% use_T_tau_gradr_factor = use_T_tau_gradr_factor
 
  ! starspots
+ s% do_starspots = do_starspots
  s% fspot = fspot
  s% xspot = xspot
 
@@ -2974,6 +2975,7 @@ subroutine set_controls_for_writing(s, ierr)
  use_T_tau_gradr_factor = s% use_T_tau_gradr_factor
 
  ! starspots
+ do_starspots = s% do_starspots
  fspot = s% fspot
  xspot = s% xspot
 

diff --git a/star/private/hydro_vars.f90 b/star/private/hydro_vars.f90
@@ -387,22 +387,27 @@ end subroutine set_Teff
 
 
       subroutine set_Teff_info_for_eqns(s, skip_partials, &
-            need_atm_Psurf, need_atm_Tsurf, r_surf, L_surf, Teff, &
+            need_atm_Psurf_in, need_atm_Tsurf_in, r_surf, L_surf, Teff, &
             lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
             lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
             ierr)
          use star_utils, only: set_phot_info
          use atm_lib, only: atm_Teff
+         use starspots, only: starspot_tweak_PT, starspot_restore_PT
          type (star_info), pointer :: s
          logical, intent(in) :: skip_partials, &
-            need_atm_Psurf, need_atm_Tsurf
+            need_atm_Psurf_in, need_atm_Tsurf_in
          real(dp), intent(out) :: r_surf, L_surf, Teff, &
             lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
             lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap
          integer, intent(out) :: ierr
+         logical :: need_atm_Psurf, need_atm_Tsurf
 
          include 'formats'
 
+         need_atm_Psurf = need_atm_Psurf_in
+         need_atm_Tsurf = need_atm_Tsurf_in
+
          ierr = 0
 
          r_surf = s% r(1)
@@ -435,11 +440,21 @@ subroutine set_Teff_info_for_eqns(s, skip_partials, &
                lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
                ierr)
          else
+            ! starspot YREC routine
+            if (s% do_starspots) then
+               need_atm_Psurf = .true.
+               need_atm_Tsurf = .true.
+               call starspot_tweak_PT(s)
+            end if
             call get_surf_PT( &
                s, skip_partials, need_atm_Psurf, need_atm_Tsurf, &
                lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
                lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
                ierr)
+            ! starspot YREC routine
+            if (s% do_starspots) then
+               call starspot_restore_PT(s)
+            end if
          end if
          if (ierr /= 0) then
             if (s% report_ierr) then
@@ -785,7 +800,6 @@ subroutine get_surf_PT( &
          real(dp) :: kap_surf
          real(dp) :: M_surf
 
-
          include 'formats'
 
          ! Set up stellar surface parameters

diff --git a/star/private/starspots.f90 b/star/private/starspots.f90
@@ -0,0 +1,111 @@
+! ***********************************************************************
+!
+!   Copyright (C) 2010-2019  Meridith Joyce & The MESA Team
+!
+!   MESA is free software; you can use it and/or modify
+!   it under the combined terms and restrictions of the MESA MANIFESTO
+!   and the GNU General Library Public License as published
+!   by the Free Software Foundation; either version 2 of the License,
+!   or (at your option) any later version.
+!
+!   You should have received a copy of the MESA MANIFESTO along with
+!   this software; if not, it is available at the mesa website:
+!   http://mesa.sourceforge.net/
+!
+!   MESA is distributed in the hope that it will be useful,
+!   but WITHOUT ANY WARRANTY; without even the implied warranty of
+!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+!   See the GNU Library General Public License for more details.
+!
+!   You should have received a copy of the GNU Library General Public License
+!   along with this software; if not, write to the Free Software
+!   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+!
+! ***********************************************************************
+
+module starspots
+
+   use star_private_def
+   use const_def
+   use utils_lib
+
+   implicit none
+
+   private
+
+   public :: starspot_tweak_gradr
+   public :: starspot_tweak_PT
+   public :: starspot_restore_PT
+
+   real(dp) :: L_init
+
+contains
+
+! -------------
+! parameterized YREC routines
+! MESA models a pressure contrast rather than temperature contrast
+! -------------
+
+   subroutine starspot_tweak_gradr(s, P, gradr, gradr_spot)
+      ! adjusts the gradient of the radius to account for starspots.
+      ! This subroutine is called at the beginning of Get_results()
+      ! in turb_support.f90
+      ! ------------------------------------------------------------
+      use auto_diff_support
+      type(star_info), pointer :: s
+      type(auto_diff_real_star_order1), intent(in) :: P
+      type(auto_diff_real_star_order1), intent(in) :: gradr
+      type(auto_diff_real_star_order1), intent(out) :: gradr_spot
+
+      real(dp) :: mu_ideal_gas, R2, Teff_local, PB_i
+      type(auto_diff_real_star_order1) :: xspot_of_r ! xspot4
+
+      if (.not. s%doing_relax .and. .not. s%doing_first_model_of_run) then
+         mu_ideal_gas = s%mu(1)  !1.00794d0 ! for hydrogen, 1 gram per mole
+         R2 = pow2(s%R(1))
+         Teff_local = pow(s%L(1)/(pi4*boltz_sigma*R2), 0.25d0)
+         PB_i = (cgas*s%rho(1)/mu_ideal_gas)*(1.0 - s%xspot)*Teff_local
+         xspot_of_r = (P - PB_i)/P
+         gradr_spot = gradr/(s%fspot*pow4(xspot_of_r) + 1d0 - s%fspot)
+      else
+         gradr_spot = gradr
+      end if
+   end subroutine starspot_tweak_gradr
+
+   subroutine starspot_tweak_PT(s)
+      ! saves the surface luminosity and adjusts it and the effective
+      ! temperature to account for starspots.
+      ! This subroutine is called at the beginning of get_surf_PT()
+      ! in hydro_vars.f90
+      ! ------------------------------------------------------------
+
+      type(star_info), pointer :: s
+
+      real(dp) :: alp
+
+      alp = 1d0 - s%fspot + s%fspot*pow4(s%xspot)
+
+      ! This is the surface-average value for luminosity
+      L_init = s%L(1)
+
+      ! Set the surface L to the unspotted, ambient L
+      s%L(1) = s%L(1)/alp
+
+      ! Now, set the Teff. Used in atm table lookup to set boundary conditions
+      s%Teff = pow(s%L(1)/(pi4*pow2(s%r(1))*boltz_sigma), 0.25_dp)
+
+   end subroutine starspot_tweak_PT
+
+   subroutine starspot_restore_PT(s)
+      ! restores the surface luminosity effective temeperature.
+      ! Called at the end of get_surf_PT()
+      ! ------------------------------------------------------------
+
+      type(star_info), pointer :: s
+
+      s%Teff = pow(L_init/(pi4*pow2(s%r(1))*boltz_sigma), 0.25_dp)
+      s%L(1) = L_init
+
+   end subroutine starspot_restore_PT
+
+end module starspots
diff --git a/star/private/turb_support.f90 b/star/private/turb_support.f90
@@ -105,13 +105,14 @@ end subroutine get_gradT
 
    subroutine do1_mlt_eval( &
          s, k, MLT_option, gradL_composition_term, &
-         gradr, grada, scale_height, mixing_length_alpha, &
+         gradr_in, grada, scale_height, mixing_length_alpha, &
          mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr)
       use chem_def, only: ih1
+      use starspots, only: starspot_tweak_gradr
       type (star_info), pointer :: s
       integer, intent(in) :: k
       character (len=*), intent(in) :: MLT_option
-      type(auto_diff_real_star_order1), intent(in) :: gradr, grada, scale_height
+      type(auto_diff_real_star_order1), intent(in) :: gradr_in, grada, scale_height
       real(dp), intent(in) :: gradL_composition_term, mixing_length_alpha
       integer, intent(out) :: mixing_type
       type(auto_diff_real_star_order1), intent(out) :: &
@@ -120,9 +121,11 @@ subroutine do1_mlt_eval( &
 
       real(dp) :: cgrav, m, XH1, gradL_old, grada_face_old
       integer :: iso, old_mix_type
-      type(auto_diff_real_star_order1) :: r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp
+      type(auto_diff_real_star_order1) :: gradr, r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp
       include 'formats'
       ierr = 0
+
+      gradr = gradr_in
 
       cgrav = s% cgrav(k)
       m = s% m_grav(k)
@@ -145,7 +148,12 @@ subroutine do1_mlt_eval( &
             iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
             s% alpha_semiconvection, s% thermohaline_coeff, &
             mixing_type, gradT, Y_face, mlt_vc, D, Gamma, ierr)
-      else         
+      else
+         ! starspot YREC routine
+         if (s% do_starspots) then
+            dV = 0d0 ! dV = 1/rho - 1/rho_start and we assume rho = rho_start.
+            call starspot_tweak_gradr(s, P, gradr_in, gradr)
+         end if
          call Get_results(s, k, MLT_option, &
             r, L, T, P, opacity, rho, dV, chiRho, chiT, Cp, gradr, grada, scale_height, &
             iso, XH1, cgrav, m, gradL_composition_term, mixing_length_alpha, &
@@ -228,7 +236,6 @@ subroutine Get_results(s, k, MLT_option, &  ! NOTE: k=0 is a valid arg
       end if
 
 
-
       ! check if this particular k can be done with TDC
       using_TDC = .false.
       if (s% MLT_option == 'TDC') using_TDC = .true.

diff --git a/star/test_suite/starspots/README.rst b/star/test_suite/starspots/README.rst
@@ -34,6 +34,6 @@ A plot of the HR diagram for the problem is included below:
    :width: 100%
 
 
-Last-Updated: 06Aug2024 by Meridith Joyce
+Last-Updated: 09Aug2024 by Meridith Joyce
 
 Last-Run: 06Aug2024 (MESA 7890d305) by pmocz on C916PXT6XW in 629 seconds using 8 threads.
diff --git a/star/test_suite/starspots/inlist_starspots b/star/test_suite/starspots/inlist_starspots
@@ -61,10 +61,8 @@
    star_history_name = 'history.data'
 
    max_age = 24.63d9 !26d9 !0.8d9
-
-   use_other_mlt_results = .true.
-   use_other_surface_PT = .true.
-
+
+   do_starspots = .true.
    fspot = 0.34d0
    xspot = 0.85d0