Merge pull request #612 from MESAHub/EbF/reimplement_fixed_atm_surface_options

reimplement fixed atm options

Author: Debraheem

docs/source/changelog.rst

@@ -99,6 +99,10 @@ starting models has been recalculated with a more stringent stopping condition,
 and now all pre-computed ZAMS models have a central hydrogen mass fraction very 
 near 0.697.
 
+The ``fixed_Teff``, ``fixed_Tsurf``, ``fixed_Psurf``,  and ``fixed_Psurf_and_Tsurf``
+atmosphere options have been reimplemented, although we caution users that their
+implementation might conflict with mlt_option = ``TDC``.
+
 Changes in r23.05.1
 ===================
 

star/defaults/controls.defaults

@@ -4173,7 +4173,8 @@
       ! ~~~~~~~~~~
 
       ! Controls how the surface temperature Tsurf and pressure Psurf are evaluated when
-      ! setting up outer boundary conditions
+      ! setting up outer boundary conditions. We caution that the use of ``'fixed_'`` atmosphere
+      ! options might conflict with mlt_option = ``TDC``.
 
       ! + ``'T_tau'``: 
       !   set Tsurf and Psurf by solving for the atmosphere structure given a T(tau) relation.
@@ -4198,6 +4199,27 @@
       !   See also the ``atm_irradiated_opacity``, ``atm_irradiated_errtol``, ``atm_irradiated_T_eq``,
       !   ``atm_irradiated_T_eq``, ``atm_irradiated_kap_v``, ``atm_irradiated_kap_v_div_kap_th``,
       !   ``atm_irradiated_P_surf`` and ``atm_irradiated_max_tries`` controls.
+      !
+      ! + ``'fixed_Teff'`` : 
+      ! set Tsurf from Eddington T(tau) relation
+      !     for current surface tau and Teff = ``atm_fixed_Teff``.
+      ! set Psurf = Radiation_Pressure(Tsurf)
+      !
+      ! + ``'fixed_Tsurf'`` :
+      ! get value of Tsurf from control parameter ``atm_fixed_Tsurf``.
+      ! set Teff from Eddington T(tau) relation for given Tsurf and tau=2/3
+      ! set Psurf = Radiation_Pressure(Tsurf)
+      !      
+      ! + ``'fixed_Psurf'`` :
+      ! get value of Psurf from control parameter ``atm_fixed_Psurf``.
+      ! set Tsurf from L and R using ``L = 4*pi*R^2*boltz_sigma*T^4``.
+      ! set Teff using Eddington T(tau) relation for tau=2/3 and T=Tsurf.
+      !
+      ! + ``'fixed_Psurf_and_Tsurf'`` :
+      ! get value of Psurf from control parameter ``atm_fixed_Psurf``.
+      ! get value of Tsurf from control parameter ``atm_fixed_Tsurf``.
+      ! see the ``conductive_flame`` test_suite for an example of 
+      ! this boundary condition implemented via the ``other_surface_PT`` hook.
 
       ! ::
 

star/private/hydro_vars.f90

@@ -781,13 +781,19 @@ subroutine get_surf_PT( &
          real(dp) :: T_surf
          real(dp) :: P_surf_atm
          real(dp) :: P_surf
+         real(dp) :: Pextra
+         real(dp) :: kap_surf
+         real(dp) :: M_surf
+
 
          include 'formats'
 
          ! Set up stellar surface parameters
 
          L_surf = s% L(1)
          R_surf = s% r(1)
+         kap_surf = s% opacity(1)
+         M_surf = s% m(1)
          Teff = s% Teff
 
          ! Initialize partials
@@ -823,38 +829,164 @@ subroutine get_surf_PT( &
             endif
 
          else
-
             ! Evaluate temperature and pressure based on atm_option
+            ! (yes, we might want to translate atm_option into an integer,
+            ! but these string comparisons are cheap)
+
+            ! The first few are special, 'trivial-atmosphere' options
+
+            select case (s% atm_option)
+         
+            case ('fixed_Teff')
+
+               ! set Tsurf from Eddington T-tau relation
+               !     for current surface tau and Teff = `atm_fixed_Teff`.
+               ! set Psurf = Radiation_Pressure(Tsurf)
+
+               Teff = s% atm_fixed_Teff
+               Teff4 = Teff*Teff*Teff*Teff
+               T_surf = pow(3._dp/4._dp*Teff4*(tau_surf + 2._dp/3._dp), 0.25_dp)
+               lnT_surf = log(T_surf)
+               lnP_surf = Radiation_Pressure(T_surf)
+
+               if (.not. skip_partials) then
+                  dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
+                  dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
+               endif
 
-            if (.false. .and. 1 == s% solver_test_partials_k .and. &
-                  s% solver_iter == s% solver_test_partials_iter_number) then
-               star_debugging_atm_flag = .true.
-            end if
+            case ('fixed_Tsurf')
 
-            call get_atm_PT( &
-                 s, tau_surf, L_surf, R_surf, s% m(1), s% cgrav(1), skip_partials, &
-                 Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
-                 lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
-                 ierr)
-            if (ierr /= 0) then
-               if (s% report_ierr) then
-                  write(*,1) 'tau_surf', tau_surf
-                  write(*,1) 'L_surf', L_surf
-                  write(*,1) 'R_surf', R_surf
-                  write(*,1) 'Teff', Teff
-                  write(*,1) 's% m(1)', s% m(1)
-                  write(*,1) 's% cgrav(1)', s% cgrav(1)
-                  write(*,*) 'failed in get_atm_PT'
+               ! set Teff from Eddington T-tau relation for given
+               ! Tsurf and tau=2/3 set Psurf =
+               ! Radiation_Pressure(Tsurf)
+
+               T_surf = s% atm_fixed_Tsurf
+               lnT_surf = log(T_surf)
+               T_surf4 = T_surf*T_surf*T_surf*T_surf
+               Teff = pow(4._dp/3._dp*T_surf4/(tau_surf + 2._dp/3._dp), 0.25_dp)
+               lnP_surf = Radiation_Pressure(T_surf)
+
+               if (.not. skip_partials) then
+                  dlnT_dL = 0._dp; dlnT_dlnR = 0._dp; dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
+                  dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
+               endif
+
+            case ('fixed_Psurf')
+
+               ! set Teff from L and R using L = 4*pi*R^2*boltz_sigma*T^4.
+               ! set Tsurf using Eddington T-tau relation
+
+               if (L_surf < 0._dp) then
+                  if (s% report_ierr) then
+                     write(*,2) 'get_surf_PT: L_surf <= 0', s% model_number, L_surf
+                     call mesa_error(__FILE__,__LINE__)
+                  end if
+                  s% retry_message = 'L_surf < 0'
+                  ierr = -1
+                  return
                end if
-               return
-            end if
 
-            if (.false. .and. 1 == s% solver_test_partials_k .and. &
-                  s% solver_iter == s% solver_test_partials_iter_number) then
-               s% solver_test_partials_val = atm_test_partials_val
-               s% solver_test_partials_dval_dx = atm_test_partials_dval_dx
-            end if
+               lnP_surf = safe_log(s% atm_fixed_Psurf)
+               if (R_surf <= 0._dp) then
+                  T_surf4 = 1._dp
+                  T_surf = 1._dp
+                  lnT_surf = 0._dp
+                  if (.not. skip_partials) then
+                     dlnT_dlnR = 0._dp
+                     dlnT_dL = 0._dp
+                  endif
+                  Teff = s% T(1)
+               else
+                  Teff = pow(L_surf/(4._dp*pi*R_surf*R_surf*boltz_sigma), 0.25_dp)
+                  T_surf4 = 3d0/4d0*pow(Teff, 4.d0)*(tau_surf + 2._dp/3._dp)
+                  T_surf = pow(T_surf4, 0.25_dp)
+                  lnT_surf = log(T_surf)
+                  if (.not. skip_partials) then
+                     dlnT_dlnR = -0.5_dp
+                     dlnT_dL = 1._dp/(4._dp*L_surf)
+                  endif
+               end if
+
+               if (.not. skip_partials) then
+                  dlnT_dlnM = 0._dp; dlnT_dlnkap = 0._dp
+                  dlnP_dL = 0._dp; dlnP_dlnR = 0._dp; dlnP_dlnM = 0._dp; dlnP_dlnkap = 0._dp
+               endif
+
+            case ('fixed_Psurf_and_Tsurf')
+
+               lnP_surf = safe_log(s% atm_fixed_Psurf)
+               T_surf = s% atm_fixed_Tsurf
+               lnT_surf = log(T_surf)
+               T_surf4 = T_surf*T_surf*T_surf*T_surf
+               Teff = pow(4d0/3d0*T_surf4/(tau_surf + 2d0/3d0), 0.25d0)
+
+               if (.not. skip_partials) then
+                  dlnT_dL = 0; dlnT_dlnR = 0; dlnT_dlnM = 0; dlnT_dlnkap = 0
+                  dlnP_dL = 0; dlnP_dlnR = 0; dlnP_dlnM = 0; dlnP_dlnkap = 0
+               endif
+
+            case default
+
+                  ! Everything else -- the 'non-trivial atmospheres' ---
+                  ! gets passed to atm_support
+
+               if (.false. .and. 1 == s% solver_test_partials_k .and. &
+                     s% solver_iter == s% solver_test_partials_iter_number) then
+                  star_debugging_atm_flag = .true.
+               end if
 
+               call get_atm_PT( &
+                  s, tau_surf, L_surf, R_surf, s% m(1), s% cgrav(1), skip_partials, &
+                  Teff, lnT_surf, dlnT_dL, dlnT_dlnR, dlnT_dlnM, dlnT_dlnkap, &
+                  lnP_surf, dlnP_dL, dlnP_dlnR, dlnP_dlnM, dlnP_dlnkap, &
+                  ierr)
+               if (ierr /= 0) then
+                  if (s% report_ierr) then
+                     write(*,1) 'tau_surf', tau_surf
+                     write(*,1) 'L_surf', L_surf
+                     write(*,1) 'R_surf', R_surf
+                     write(*,1) 'Teff', Teff
+                     write(*,1) 's% m(1)', s% m(1)
+                     write(*,1) 's% cgrav(1)', s% cgrav(1)
+                     write(*,*) 'failed in get_atm_PT'
+                  end if
+                  return
+               end if
+
+               if (.false. .and. 1 == s% solver_test_partials_k .and. &
+                     s% solver_iter == s% solver_test_partials_iter_number) then
+                  s% solver_test_partials_val = atm_test_partials_val
+                  s% solver_test_partials_dval_dx = atm_test_partials_dval_dx
+               end if
+            end select
+         end if
+         
+         ! if using fixed surface, calculate Pextra.
+         if (s% atm_option == 'fixed_Tsurf' .or. s% atm_option == 'fixed_Psurf_and_Tsurf' &
+            .or. s% atm_option == 'fixed_Psurf' .or. s% atm_option == 'fixed_Teff') then
+            ! add extra pressure for fixed atmosphere options
+            if (s% Pextra_factor /= 0._dp) then
+               P_surf_atm = exp(lnP_surf)
+               Pextra = s% Pextra_factor*(kap_surf/tau_surf)*(L_surf/M_surf)/(6._dp*pi*clight*s% cgrav(1))
+               P_surf = P_surf_atm + Pextra
+               if (P_surf < 1E-50_dp) then
+                  lnP_surf = -50*ln10
+                  if (.not. skip_partials) then
+                     dlnP_dL = 0._dp
+                     dlnP_dlnR = 0._dp
+                     dlnP_dlnM = 0._dp
+                     dlnP_dlnkap = 0._dp
+                  endif
+               else
+                  lnP_surf = log(P_surf)
+                  if (.not. skip_partials) then
+                     dlnP_dL = dlnP_dL*P_surf_atm/P_surf
+                     dlnP_dlnR = dlnP_dlnR*P_surf_atm/P_surf
+                     dlnP_dlnM = dlnP_dlnM*P_surf_atm/P_surf
+                     dlnP_dlnkap = dlnP_dlnkap*P_surf_atm/P_surf
+                  endif
+               end if
+            end if
          end if
 
          ! Check outputs

star/test_suite/conductive_flame/README.rst

@@ -29,7 +29,8 @@ controlled from the inlist.
 
 The inner boundary is at r = 0.  The outer boundary has a fixed
 temperature and a fixed pressure equal to the initial pressure of the
-material.  This is achieved via the ``use_other_surface_PT`` hook.
+material.  This is achieved via the ``use_other_surface_PT`` hook, but
+can also be done using the ``fixed_Psurf_and_Tsurf`` atmosphere option.
 
 After an initial transient, the entire flame structure, approximately
 isobaric, propagates into the upstream fuel with a unique speed and
@@ -55,5 +56,5 @@ routine ``flame_properties`` in the ``run_star_extras.f90``.
 .. |Schwab2020| replace:: `Schwab et al. (2020) <https://ui.adsabs.harvard.edu/abs/2020ApJ...891....5S/abstract>`__
 
 
-Last-Updated: 2021-06-21 (mesa b2364463) by Josiah Schwab
+Last-Updated: 2021-06-21 (mesa b2364463) by Josiah Schwab, + documentation 2024-01-22 EbF