Merge pull request #646 from MESAHub/kap_derivative_testing

Automatic differentiation for cubic interpolation across type-1 and type-2 kap tables

Author: evbauer

docs/source/changelog.rst

@@ -21,10 +21,9 @@ New Features
 ``max_allowed_nz`` is now ignored if the value is less than or equal to zero.
 
 Kap
-~~~~~
+~~~
 
-`High Temperature Opacity Tables`
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**High Temperature Opacity Tables**
 
 Type 1 Rosseland-mean opacity tables from The Los Alamos
 OPLIB database (`Colgan et al. 2016 <https://ui.adsabs.harvard.edu/abs/2016ApJ...817..116C/abstract>`_) are now available (Farag et al. 2024).
@@ -52,8 +51,8 @@ implementation of these tables. For further details on these new OPLIB opacity t
 the Type 1 OPAL/OP tables as well as their effect on solar models can be found in
 in Farag et al. 2024.
 
-`Low Temperature Opacity Tables`
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+**Low Temperature Opacity Tables**
 
 Low temperature Rosseland-mean opacity tables for both (`AAG21, Asplund et al. 2021 <https://ui.adsabs.harvard.edu/abs/2021A%26A...653A.141A/abstract>`_),
 and (`MB22, Magg et al. 2022 <https://doi.org/10.1051/0004-6361/202142971>`_) 
@@ -65,6 +64,65 @@ options for :ref:`kap/defaults:kap_lowT_prefix`:
 + ``'lowT_fa05_mb22'``
 + ``'lowT_fa05_aag21'``
 
+**Opacity interpolation**
+
+We have updated the opacity interpolation scheme to provide much higher quality derivatives when doing cubic interpolation
+in composition.
+
+MESA interpolates across opacity tables in the :math:`X–Z` plane through the use of two consequtive 1D splines.
+MESA offers users the ability to choose linear or cubic interpolation for these splines, 
+while leaving the default as linear interpolation::
+
+  cubic_interpolation_in_X = .false.
+  cubic_interpolation_in_Z = .false.
+
+This choice of default was primarily due to the fact that
+the previous cubic composition interpolation scheme in MESA suffered from poor quality interpolated opacity derivatives with respect to
+density and temperature, which often disagreed with the numerical derivatives produced via nearest neighbor
+Richardson extrapolation. The figure below shows this comparison on a logarithmic scale, where in general red indicates poor quality
+derivatives and blue indicates high quality derivatives.
+
+.. figure:: changelog_plots/cubic_dfridr_dkapdT.png
+   :alt: old cubic relative kap derivative error
+
+   This figure shows the logarithmic relative error in the derivative :math:`\partial \kappa / \partial T` (:math:`X` = 0.625, :math:`Z` = 0.015),
+   for an OPAL opacity table grid using Grevesse & Sauval (1998) abundances, generated from MESA’s kap module, using the previous cubic interpolation scheme.
+   The OPLIB log(:math:`R`) = −8, 1.5 table boundaries are marked with a solid black line and the OPAL/OP log(:math:`R`) = 1.0 boundary is shown with a dashed line.
+   The approximate location of the Z-dependent transition to an electron conduction dominated opacity is marked with dot-dash blue curve.
+   Regions for Atomic, molecular, and compton scattering opacity are labeled and presented with their associated blending regions.
+
+
+While the opacity derivatives do not directly appear in the canonical equations of stellar structure, they do appear in the Jacobian matrix for MESA's implicit solver.
+Numerically unstable opacity derivatives can halt the progress of the solver and ultimately crash a calculation. 
+
+To improve the numerical stability of MESA's cubic opacity interpolation routines, we have implemented
+automatic differentiation into the opacity interpolating functions. Now, when using cubic interpolation, the opacity derivatives for an arbitrary mixture
+in the :math:`X–Z` plane are computed by taking the derivative of the interpolating function as opposed to the interpolant of the derivatives. This improvement
+has led to a significant reduction in the relative derivative error and an increase in the numerical accuracy of opacity derivatives computed with cubic interpolation. 
+
+.. figure:: changelog_plots/cubic_dfridr_dkapdT_ad.png
+   :alt: new cubic relative kap derivative error
+
+   Same as previous figure, but for new cubic interpolation scheme taking advantage of automatic differentiation.
+
+
+This new implementation of cubic interpolation in composition for opacity tables comes close to achieving the derivative quality of the linear interpolation
+option (shown below), while also providing more accurate opacity physics between opacity table grid points.
+
+.. figure:: changelog_plots/linear_dfridr_dkapdT.png
+   :alt: linear relative kap derivative error
+
+   Same as previous figure, but for linear interpolation instead of cubic.
+
+
+For this MESA release, linear interpolation remains the default method for interpolating in composition between opacity tables
+while we continue to investigate the residual areas where cubic interpolation appears to occasionally produce lower quality derivatives.
+However, adopting cubic interpolation has been shown to consistently increase the overall 
+opacity of a model, and can directly effect the structure of solar models, see Appendix B & C in Farag et al. 2024.
+We anticipate making cubic interpolation the default in a future MESA release version. 
+We encourage users to experiment with these different opacity interpolation routines and be mindful of the effect they can have on their stellar models.
+
+
 Chem
 ~~~~~
 New initial metal mass fractions ``initial_zfracs`` taken from photospheric estimates of the solar heavy element abundances in (AAG21, Asplund et al. 2021) and (MB22, Magg et al. 2022)

docs/source/changelog_plots/cubic_dfridr_dkapdT.png

@@ -15,7 +15,9 @@ include $(MESA_DIR)/utils/makefile_header
 SRCS = \
    interp_1d_misc.f90 \
    interp_1d_pm.f90 \
+   interp_1d_pm_autodiff.f90 \
    interp_1d_mp.f90 \
+   interp_1d_mp_autodiff.f90\
    interp_1d_misc_sg.f90 \
    interp_1d_pm_sg.f90 \
    interp_1d_mp_sg.f90 \
@@ -28,9 +30,9 @@ SRCS = \
 # LIBRARIES
 
 ifeq ($(USE_SHARED), YES)
-     LIBS_OTHER = $(LIBS_MATRIX)
+     LIBS_OTHER = auto_diff $(LIBS_MATRIX)
      DEPS_OTHER = $(patsubst %,$(MESA_LIB_DIR)/lib%.$(LIB_SUFFIX),$(LIBS_OTHER))
-     LOAD_OTHER = -L$(MESA_LIB_DIR) $(LOAD_MATRIX)
+     LOAD_OTHER = -L$(MESA_LIB_DIR) -lauto_diff $(LOAD_MATRIX)
 endif
 
 #################################################################

docs/source/changelog_plots/cubic_dfridr_dkapdT_ad.png

@@ -26,7 +26,8 @@
       module interp_1d_misc
 
       use const_lib, only: dp
-
+      use auto_diff
+      
       implicit none
 
       contains
@@ -181,7 +182,75 @@ subroutine do_interp_values(init_x, nx, f1, nv, x, vals, ierr)
 
 
       end subroutine do_interp_values
+
+      subroutine do_interp_values_autodiff(init_x, nx, f1, nv, x, vals, ierr)
+         type(auto_diff_real_2var_order1), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
+         integer, intent(in) :: nx ! length of init_x vector
+         type(auto_diff_real_2var_order1), intent(in), pointer :: f1(:) ! =(4, nx)  ! data & interpolation coefficients
+         integer, intent(in) :: nv ! length of new x vector and vals vector
+         type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nv)  ! locations where want interpolated values
+            ! strictly monotonic in same way as init_x
+            ! values out of range of init_x's are clipped to boundaries of init_x's
+         type(auto_diff_real_2var_order1), intent(inout) :: vals(:) ! (nv)
+         integer, intent(out) :: ierr ! 0 means aok
+   
+         integer :: k_old, k_new
+         type(auto_diff_real_2var_order1) :: xk_old, xkp1_old, xk_new, delta
+         logical :: increasing
+         type(auto_diff_real_2var_order1), pointer :: f(:,:) ! (4, nx)  ! data & interpolation coefficients
+         
+         ierr = 0
+         
+         if (nx == 1) then            
+            vals(1:nv) = f1(1)
+            return
+         end if
+
+         f(1:4,1:nx) => f1(1:4*nx)
+
+         if(init_x(1) < init_x(nx)) then
+            increasing = .true.
+         end if
+         
+         k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
+                  
+         do k_new = 1, nv
 
+            xk_new = x(k_new)
+            if (increasing) then
+               if (xk_new > init_x(nx)) then
+                  xk_new = init_x(nx)
+               else if (xk_new < init_x(1)) then
+                  xk_new = init_x(1)
+               end if
+            else ! decreasing
+               if (xk_new < init_x(nx)) then
+                  xk_new = init_x(nx)
+               else if (xk_new > init_x(1)) then
+                  xk_new = init_x(1)
+               end if
+            end if
+            do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
+               k_old = k_old + 1
+               if (k_old >= nx) then
+                  k_old = k_old - 1
+                  xk_new = xkp1_old
+                  exit
+               end if
+               xk_old = xkp1_old
+               xkp1_old = init_x(k_old+1)
+            end do
+         
+            delta = xk_new - xk_old
+         
+            vals(k_new) =  &
+                  f(1, k_old) + delta*(f(2, k_old)  &
+                     + delta*(f(3, k_old) + delta*f(4, k_old)))
+            
+         end do
+
+   
+      end subroutine do_interp_values_autodiff
 
       subroutine do_interp_values_and_slopes(init_x, nx, f1, nv, x, vals, slopes, ierr)
          real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
@@ -576,10 +645,10 @@ subroutine minmod4(z, n, f1, f2, f3, f4)
          call minmod(z, n, z, f3)
          call minmod(z, n, z, f4)
       end subroutine minmod4
-      
-      
+
+
       subroutine median(z, n, f1, f2, f3)
-         real(dp), intent(inout) :: z(:)     
+         real(dp), intent(inout) :: z(:)
          integer, intent(in) :: n       ! length of vectors
          real(dp), intent(in) :: f1(:), f2(:), f3(:)
          real(dp), target :: tmp1_ary(n), tmp2_ary(n)
@@ -593,4 +662,61 @@ subroutine median(z, n, f1, f2, f3)
       end subroutine median
 
 
+      type(auto_diff_real_2var_order1) function minmod1_autodiff(f1, f2)
+         use auto_diff
+         type(auto_diff_real_2var_order1), intent(in) :: f1, f2
+         minmod1_autodiff = 0.5d0 * (sign(f1) + sign(f2)) * min(abs(f1), abs(f2))
+      end function minmod1_autodiff
+
+
+      type(auto_diff_real_2var_order1) function median1_autodiff(f1, f2, f3)
+         use auto_diff
+         type(auto_diff_real_2var_order1), intent(in) :: f1, f2, f3
+         median1_autodiff = f1 + minmod1_autodiff(f2 - f1, f3 - f1)
+      end function median1_autodiff
+
+
+      subroutine minmod_autodiff(z, n, f1, f2)
+         use auto_diff
+         type(auto_diff_real_2var_order1), intent(inout) :: z(:)
+         integer, intent(in) :: n       ! length of vectors
+         integer :: k
+         type(auto_diff_real_2var_order1), intent(in) :: f1(:), f2(:)
+         do k =1, n
+         z(k) = 0.5d0 * (sign(f1(k)) + sign(f2(k))) * min(abs(f1(k)), abs(f2(k)))
+         end do
+      end subroutine minmod_autodiff
+
+      subroutine minmod4_autodiff(z, n, f1, f2, f3, f4)
+         use auto_diff
+         type(auto_diff_real_2var_order1), intent(inout) :: z(:)
+         integer, intent(in) :: n       ! length of vectors
+         type(auto_diff_real_2var_order1), intent(in) :: f1(:), f2(:), f3(:), f4(:)
+         call minmod_autodiff(z, n, f1, f2)
+         call minmod_autodiff(z, n, z, f3)
+         call minmod_autodiff(z, n, z, f4)
+      end subroutine minmod4_autodiff
+
+
+      subroutine median_autodiff(z, n, f1, f2, f3)
+         use auto_diff
+         type(auto_diff_real_2var_order1), intent(inout) :: z(:)
+         integer, intent(in) :: n       ! length of vectors
+         integer :: k
+         type(auto_diff_real_2var_order1), intent(in) :: f1(:), f2(:), f3(:)
+         type(auto_diff_real_2var_order1), target :: tmp1_ary(n), tmp2_ary(n)
+         type(auto_diff_real_2var_order1), pointer :: tmp1(:), tmp2(:)
+         tmp1 => tmp1_ary
+         tmp2 => tmp2_ary
+         do k =1,n
+         tmp1(k) = f2(k) - f1(k)
+         tmp2(k) = f3(k) - f1(k)
+         end do
+         call minmod_autodiff(z, n, tmp1, tmp2)
+         do k =1,n
+         z(k) = z(k) + f1(k)
+         end do
+      end subroutine median_autodiff
+
+
       end module interp_1d_misc