Smoother Delta_Pg

[BucheleL]

When looking at the evolution of delta_Pg (with the previous implementation optimized for solar-like oscillators not the new traditional option included in the release candidate) there can be discontinuities in the evolution which seem to be related to the mesh. Here's an example of a 0.8Msol RGB where the blue line shows the MESA default and the orange line shows my attempt at solving this issue. You can clearly see a jump in the MESA value around the dashed grey line.

To see what is happening, I've plotted (very zoomed in) propagation diagrams for the models on either side of this discontinuity, where the points on the $N^2$ line show the MESA grid points.

These show that the discontinuity is due to the value of $\omega_{max}$ crossing to the next grid point and expanding the integration area used to calculate delta_Pg. I saw similar behavior at the outer turning point as well. To try and smooth out these jumps, I've written some code that interpolates the exact point where $N^2 = \omega_{max}^2$ and begins and ends the integral there:

         delta_Pi1 = 0._dp
         if (.not. s% calculate_Brunt_N2) return
         entered_g_mode_cavity = .false.
         k = 0
         r = 0._dp
         N2 = 0._dp
         S2 = 0._dp
         dr = 0._dp
         rho = 0._dp
         omega_max2 = pow2(2*pi*s% nu_max/1d6)

         ! Modifed slightly to include partial cell with face at actual inner turning point to avoid jumps due to mesh points 
         ! Find k of turning point 

         k_tpu = 0
         k_tpl = 0

         do k = s% nz, 2, -1
            N2 = s% brunt_N2(k) ! brunt_N2 at cell face 
            S2 = 2*pow2(s% csound_face(k)/s% r(k)) ! l=1 Lamb Freq ^2 at cell face
            if ((omega_max2 < N2) .and. (omega_max2 <S2)) then ! Check if in g-mode cavity at nu_max
               k_tpu = k
               k_tpl = k + 1
               exit
            end if
         end do

         if (k_tpu == 0) return  ! No g-mode cavity at nu_max 

         if (k_tpu /= s% nz) then ! If cavity boundary is not at center of star  

         ! calculate exact r of inner turning point assuming N2 vs r is linear in log-log space 
            r_tp = (log10(omega_max2) - log10(s% brunt_N2(k_tpl)))*(log10(s% r(k_tpu))-log10(s% r(k_tpl)))&
                                                          /(log10(s% brunt_N2(k_tpu))-log10(s% brunt_N2(k_tpl))) + log10(s% r(k_tpl))
            r_tp = 10**(r_tp)

            dr = s% r(k_tpu) - r_tp
            N2 = omega_max2 ! Since exactly at boundary 
            rho = (s% rho(k_tpu) - s% rho(k_tpl))/(s% r(k_tpu) - s% r(k_tpl))*(r_tp - s% r(k_tpl)) + s% rho(k_tpl)
            delta_Pi1 = delta_Pi1 + sqrt(N2)/r_tp*dr

         end if

         ! we integrate at cell edges 
         do k = k_tpu, 2, -1 ! Integrate backwards so that integration is over the innermost g-cavity
               N2 = s% brunt_N2(k) ! brunt_N2 at cell face 
               S2 = 2*pow2(s% csound_face(k)/s% r(k)) ! l=1 Lamb Freq ^2 at cell face
            if ((omega_max2 < N2) .and. (omega_max2 <S2)) then ! Check if in g-mode cavity at nu_max
               entered_g_mode_cavity = .true.
               r = s% r(k) ! r evaulated at cell face 
               rho = s% rho_face(k) ! r evaulated at cell face
               dr = s% rmid(k-1) - s%rmid(k) ! dr evaulated at cell face 
               delta_Pi1 = delta_Pi1 + sqrt(N2)/r*dr
            else
               if (entered_g_mode_cavity) then
                   ! calculate exact r of outer turning point
                   k_tpu = k
                   k_tpl = k+1
                   r_tp = (log10(omega_max2) - log10(s% brunt_N2(k_tpl)))*(log10(s% r(k_tpu))-log10(s% r(k_tpl)))&
                                             /(log10(s% brunt_N2(k_tpu))-log10(s% brunt_N2(k_tpl))) + log10(s% r(k_tpl))
                   r_tp = 10**(r_tp)
                   dr = r_tp - s% r(k_tpl)
                   N2 = omega_max2 ! Since exactly at boundary 
                   rho = (s% rho(k_tpu) - s% rho(k_tpl))/(s% r(k_tpu) - s% r(k_tpl))*(r_tp - s% r(k_tpl)) + s% rho(k_tpl)
                   delta_Pi1 = delta_Pi1 + sqrt(N2)/r_tp*dr
                   exit ! avoid integrating over other g-mode cavities such as near-surface layers
               endif
            end if
         end do

         delta_Pi1 = sqrt(2._dp)*pi2/delta_Pi1

This code results in the orange line plotted in the first figure above, smoothing out the discontinuity. The only thing I'm not sure about in this code is if I'm handling the calculation of dr correctly for the new interpolated grid points since for the other points dr is calculated from rmid rather than r.

[Debraheem]

Thanks @BucheleL, this is very interesting! Are you able to test the standard form (from the pre release) on this model as well, just for reference?

[BucheleL]

Hi @Debraheem, sorry for the delay. I've run this test (albeit with a slight different model). It looks like the traditional form does better than the previous default.

Looking at the partial derivatives, it's also pretty good. Although there is still some jitter, but it's really difficult to say how significant this is.

I guess one could implement some interpolation at the points where $N^2 = 0$ instead of $N^2 = \omega_{max}^2$ but I would ask some folks who do g-mode pulsators to confirm if that is necessary at their data precision.