Merge pull request #500 from MESAHub/restore_convective_bdy_weight

Restore a simplified version of the convective_bdy_weight control

Author: evbauer

docs/source/changelog.rst

@@ -80,6 +80,18 @@ species to ``X``, ``Y`` or ``Z``, in which case the checks are applied
 *individually* to all isotopes of hydrogen, helium or metals,
 respectively.
 
+Resolution control convective_bdy_weight has been reintroduced
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The option to add extra resolution at convective boundaries
+with the control ``convective_bdy_weight`` was removed after
+version 12115, but has now been reintroduced in a simplified form.
+This control no longer applies to newly nonconvective zones, but
+does add resolution at the location of convective boundaries.
+This was found to be important for smooth convective boundary
+evolution with convective premixing.
+
+
 Changes in r22.11.1
 ===================
 

star/defaults/controls.defaults

@@ -6265,6 +6265,20 @@
     mesh_dlog_other_dlogP_extra = 0.25d0
     mesh_dlog_photo_dlogP_extra = 1
 
+      ! convective_bdy_weight
+      ! ~~~~~~~~~~~~~~~~~~~~~
+      ! convective_bdy_dq_limit
+      ! ~~~~~~~~~~~~~~~~~~~~~~~
+      ! convective_bdy_min_dt_yrs
+      ! ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      ! Mesh function to enhance resolution near convective boundaries
+
+      ! ::
+
+      convective_bdy_weight = 0
+      convective_bdy_dq_limit = 3d-5
+      convective_bdy_min_dt_yrs = 1d-3
 
       ! max_rel_delta_IE_for_mesh_total_energy_balance
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

star/private/ctrls_io.f90

@@ -233,6 +233,7 @@ module ctrls_io
     E_function_weight, E_function_param, P_function_weight, &
     mesh_logX_species, &
     mesh_logX_min_for_extra, mesh_dlogX_dlogP_extra, mesh_dlogX_dlogP_full_on, mesh_dlogX_dlogP_full_off, &
+    convective_bdy_min_dt_yrs, convective_bdy_dq_limit, convective_bdy_weight, &
 
     mesh_dlog_eps_min_for_extra, mesh_dlog_eps_dlogP_full_on, mesh_dlog_eps_dlogP_full_off, &
     mesh_dlog_pp_dlogP_extra, mesh_dlog_cno_dlogP_extra, mesh_dlog_3alf_dlogP_extra, &
@@ -1488,6 +1489,10 @@ subroutine store_controls(s, ierr)
  s% mesh_dlogX_dlogP_full_on = mesh_dlogX_dlogP_full_on
  s% mesh_dlogX_dlogP_full_off = mesh_dlogX_dlogP_full_off
 
+ s% convective_bdy_min_dt_yrs = convective_bdy_min_dt_yrs
+ s% convective_bdy_weight = convective_bdy_weight
+ s% convective_bdy_dq_limit = convective_bdy_dq_limit
+
  s% mesh_dlog_eps_min_for_extra = mesh_dlog_eps_min_for_extra
  s% mesh_dlog_eps_dlogP_full_on = mesh_dlog_eps_dlogP_full_on
  s% mesh_dlog_eps_dlogP_full_off = mesh_dlog_eps_dlogP_full_off
@@ -3151,6 +3156,10 @@ subroutine set_controls_for_writing(s, ierr)
  mesh_dlogX_dlogP_full_on = s% mesh_dlogX_dlogP_full_on
  mesh_dlogX_dlogP_full_off = s% mesh_dlogX_dlogP_full_off
 
+ convective_bdy_min_dt_yrs = s% convective_bdy_min_dt_yrs
+ convective_bdy_weight = s% convective_bdy_weight
+ convective_bdy_dq_limit = s% convective_bdy_dq_limit
+
  mesh_dlog_eps_min_for_extra = s% mesh_dlog_eps_min_for_extra
  mesh_dlog_eps_dlogP_full_on = s% mesh_dlog_eps_dlogP_full_on
  mesh_dlog_eps_dlogP_full_off = s% mesh_dlog_eps_dlogP_full_off

star/private/mesh_functions.f90

@@ -76,6 +76,7 @@ integer function num_mesh_functions(s)
          i = 0
          if (s% use_other_mesh_functions) &
             call s% how_many_other_mesh_fcns(s% id, i)
+         if (s% convective_bdy_weight > 0) i=i+1
          if (s% E_function_weight > 0) i=i+1
          if (s% P_function_weight > 0) i=i+1
          if (s% T_function1_weight > 0) i=i+1
@@ -175,6 +176,9 @@ subroutine set_mesh_function_data( &
          if (s% omega_function_weight > 0 .and. s% rotation_flag) then
             i = i+1; names(i) = 'omega_function'
          end if
+         if (s% convective_bdy_weight > 0) then
+            i = i+1; names(i) = 'convective_bdy'
+         end if
          do k=1,num_xa_function
             if (do_mass_function(s, s% xa_function_species(k), s% xa_function_weight(k), j)) then
                i = i+1; names(i) = trim(s% xa_function_species(k))
@@ -263,6 +267,9 @@ subroutine set_mesh_function_data( &
                   vals(k,i) = s% omega_function_weight*log10(max(1d-99,abs(s% omega(k))))
                end do
 
+            else if (names(i) == 'convective_bdy') then
+               call do_conv_bdy(i)
+
             else
                do k=1,num_xa_function
                   call do1_xa_function(k,i)
@@ -274,6 +281,37 @@ subroutine set_mesh_function_data( &
 
          contains
 
+
+         subroutine do_conv_bdy(i)
+            integer, intent(in) :: i
+            integer :: k, j
+
+            vals(1:nz,i) = 0
+
+            if (s% dt < s% convective_bdy_min_dt_yrs*secyer) return
+
+            !do k=1,nz
+            !   if (s% dq(k) < s% convective_bdy_dq_limit) cycle
+            !   if (s% cz_bdy_dq(k) /= 0) vals(k,i) = s% convective_bdy_weight
+            !end do
+
+            do j = 1, s% num_conv_boundaries
+               k = s% conv_bdy_loc(j)
+               if (s% dq(k) < s% convective_bdy_dq_limit) cycle
+               vals(k,i) = s% convective_bdy_weight
+               if (k < nz .and. s% top_conv_bdy(i)) then
+                  vals(k+1,i) = s% convective_bdy_weight
+               else if (k > 1 .and. .not. s% top_conv_bdy(i)) then
+                  vals(k-1,i) = s% convective_bdy_weight
+               end if
+            end do
+
+            do k=2,nz
+               vals(k,i) = vals(k,i) + vals(k-1,i)
+            end do
+            
+         end subroutine do_conv_bdy
+
          subroutine do1_xa_function(k,i)
             integer, intent(in) :: k,i
             real(dp) :: weight, param

star/test_suite/16M_conv_premix/inlist_16M_conv_premix

@@ -85,6 +85,8 @@
    ! RELAXED VALUE FOR TEST SUITE
    mesh_delta_coeff = 0.5 ! 0.2
 
+   convective_bdy_weight = 2
+
    ! Output
 
    photo_interval = 50 ! 500

star/test_suite/diffusion_smoothness/inlist_diffusion_smoothness

@@ -64,6 +64,8 @@
       limit_for_rel_error_in_energy_conservation = 1d-6 ! this makes d2N2 too large.  takes 63 steps. log_rel_run_E_err  -3.5859393528096946
       limit_for_rel_error_in_energy_conservation = 1d-5 ! this is smooth enough.  takes 58 steps. log_rel_run_E_err  -3.4820705862368895
       hard_limit_for_rel_error_in_energy_conservation = 1d-4
+
+      convective_bdy_weight = 1
 
       cool_wind_RGB_scheme = 'Reimers'
       cool_wind_AGB_scheme = 'Blocker'

star_data/private/star_controls.inc

@@ -597,6 +597,9 @@
          real(dp) :: mesh_dlogX_dlogP_full_on(num_mesh_logX)
          real(dp) :: mesh_dlogX_dlogP_full_off(num_mesh_logX)
 
+         real(dp) :: convective_bdy_weight, &
+            convective_bdy_dq_limit, convective_bdy_min_dt_yrs
+
          real(dp) :: mesh_dlog_eps_min_for_extra
          real(dp) :: mesh_dlog_eps_dlogP_full_on
          real(dp) :: mesh_dlog_eps_dlogP_full_off