fixed incorrect metal mass fraction array for Freedman lowT opacity t…

atm/test/test_output

@@ -85,17 +85,17 @@
 
  test_irradiated
 
-                              test_grey_irradiated: kap    4.3649748642304882D-03
+                              test_grey_irradiated: kap    4.7498368201485917D-03
                                             M/M_jupiter    1.4999999999999998D+00
                                             R/R_jupiter    1.5743735095633800D+00
                                                  R/Rsun    1.6178684913857289D-01
                                                   kap_v    4.0000000000000001D-03
                                                       P    1.0000000000000000D+06
 
-                                                    tau    2.9101126947637423D+00
+                                                    tau    3.1669763036977079D+00
                                                    Teff    5.7759999999999891D+03
                                                    T_eq    1.0000000000000000D+03
                                                   T_int    9.0000000000000000D+02
-                                                      T    1.2929349604416070D+03
-                                                   logT    3.1115766787614829D+00
+                                                      T    1.3124222808279656D+03
+                                                   logT    3.1180735947393834D+00
 

docs/source/changelog.rst

@@ -35,6 +35,8 @@ MESA no longer stops when reactions for which special rates are set are not in t
 Bug Fixes
 ---------
 
+Important bug fix for ``r26.4.1`` identified by Emily Sandford and Louis Siebenaler: the ``lowT_Freedman11`` opacity option used ``[M/H]`` labels as the metal mass fraction when interpolating in ``Z``, resulting in incorrect opacities. We recommend users who use these low-temperature opacities, such as in planet models, update to the latest MESA version or employ the fixes in :ref:`the known bugs entry <freedman_lowt_z_bug>` and `gh-993 <https://github.com/MESAHub/mesa/pull/993>`_.
+
 The parameter ``report_max_infall_inside_fe_core`` was ignored in versions r25.12.1 and r26.4.1 and always had it's default value. See `gh-981 https://github.com/MESAHub/mesa/pull/981`_.
 
 ``fe_core_infall_limit`` now obeys ``when_to_stop_rtol`` and ``when_to_stop_atol`` again (broken since r11532).

docs/source/kap/overview.rst

@@ -109,13 +109,19 @@ Tables based on the work of |Fergusson| include the effects of
 molecules and grains and cover the range
 :math:`2.7 \le \log T \le 4.5` and :math:`-8 \le \log R \le 1`.
 
-Tables based on the work of |Freedman| include the effects of
-molecules and cover the range :math:`1.88 \le \log T \le 4.5` and
-:math:`-8 \le \log R \le 9`.  The table set was privately communicated
-by R. S. Freedman in 2011.  Unlike other opacity sources, this is a 1D
-sequence of tables in :math:`Z` as opposed to a 2D grid of
-:math:`(X,Z)` values.  (The assumed H/He abundances scale with
-:math:`Z`.)
+Tables based on the work of |Freedman08| and |Freedman14| include the
+effects of molecules.  The MESA ``lowT_Freedman11`` tables cover the
+range :math:`1.88 \le \log T \le 4.5` and
+:math:`-8 \le \log R \le 9` after preprocessing.  The raw opacity data
+used to build these tables stop near :math:`T = 4000\,\rm K`,
+corresponding to :math:`\log T \simeq 3.60`; the higher temperature
+part of the MESA table is generated during preprocessing so the table
+can be used on the same grid as the other low temperature opacity
+tables.  The table set was privately communicated by R. S. Freedman in
+2011 and corresponds to the expanded opacity data later described by
+|Freedman14|.  Unlike other opacity sources, this is a 1D sequence of
+tables in :math:`Z` as opposed to a 2D grid of :math:`(X,Z)` values.
+(The assumed H/He abundances scale with :math:`Z`.)
 
 
 Tables from ÆSOPUS (|AESOPUS|) include variation factors for the CNO
@@ -182,7 +188,9 @@ degeneracy, the additional correction formulae of |B20| are applied.
 
 .. |Fergusson| replace:: `Ferguson et al. (2005) <https://ui.adsabs.harvard.edu/abs/2005ApJ...623..585F/abstract>`__
 
-.. |Freedman| replace:: `Freedman et al. (2008) <https://ui.adsabs.harvard.edu/abs/2008ApJS..174..504F/abstract>`__
+.. |Freedman08| replace:: `Freedman et al. (2008) <https://ui.adsabs.harvard.edu/abs/2008ApJS..174..504F/abstract>`__
+
+.. |Freedman14| replace:: `Freedman et al. (2014) <https://ui.adsabs.harvard.edu/abs/2014ApJS..214...25F/abstract>`__
 
 .. |AESOPUS| replace:: `Marigo & Aringer 2009 <https://ui.adsabs.harvard.edu/abs/2009A%26A...508.1539M/abstract>`__
 
@@ -193,4 +201,3 @@ degeneracy, the additional correction formulae of |B20| are applied.
 .. |OPLIB| replace:: `Colgan 2016 <https://ui.adsabs.harvard.edu/abs/2016ApJ...817..116C/abstract>`__
 
 .. |Cassisi| replace:: `Cassisi et al. (2007) <https://ui.adsabs.harvard.edu/abs/2007ApJ...661.1094C/abstract>`__
-

docs/source/known_bugs.rst

@@ -13,6 +13,27 @@ issue, but it may not be complete.
 r26.4.1
 =======
 
+.. _freedman_lowt_z_bug:
+
+Kap: ``lowT_Freedman11`` used ``[M/H]`` labels as ``Z``
+-------------------------------------------------------
+
+The ``lowT_Freedman11`` opacity option used the Freedman table labels
+``0.01``, ``0.02``, ``0.04``, ``0.10``, ``0.20``, ``0.63``, and ``1.00`` as
+metal mass fractions.  These labels correspond to ``[M/H]`` not
+MESA's metal mass fraction ``Z``.  The opacity routines incorrectly interpolated in ``Z``,
+using the ``[M/H]`` labels as the interpolation grid.
+
+This has been fixed in the main branch after ``r26.4.1``.  The fix updates the
+``lowT_Freedman11`` ``Z`` grid to use the correct corresponding metal mass fractions.
+See `gh-993 <https://github.com/MESAHub/mesa/pull/993>`_.
+
+If applying these fixes to an existing checkout, rerun the opacity preprocessor
+to regenerate the tables; this is not done by a normal MESA install::
+
+   cd $MESA_DIR/kap/preprocessor
+   ./build_data_and_export
+
 .. _report_max_infall_inside_fe_core_bug:
 
 Controls: ``report_max_infall_inside_fe_core`` is ignored

kap/preprocessor/Makefile

@@ -59,6 +59,7 @@ $(call kap,$1,0.1000,-1)
 endef
 
 define kap-co-enhanced
+$(call kap,$1,0.000,-1)
 $(call kap,$1,0.001,-1)
 $(call kap,$1,0.004,-1)
 $(call kap,$1,0.010,-1)
@@ -69,13 +70,13 @@ $(call kap,$1,0.100,-1)
 endef
 
 define kap-freedman
-$(call kap,$1,0.01,0)
-$(call kap,$1,0.02,0)
-$(call kap,$1,0.04,0)
-$(call kap,$1,0.10,0)
-$(call kap,$1,0.20,0)
-$(call kap,$1,0.63,0)
-$(call kap,$1,1.00,0)
+$(call kap,$1,0.007452,0)
+$(call kap,$1,0.014760,0)
+$(call kap,$1,0.029024,0)
+$(call kap,$1,0.069840,0)
+$(call kap,$1,0.130292,0)
+$(call kap,$1,0.321456,0)
+$(call kap,$1,0.428844,0)
 endef
 
 define kap-oplib
@@ -131,10 +132,10 @@ build-data: $(BUILD_DIR_MODULE)/kap_input_data/ logT_points $(OBJ_OUT) | $(BUILD
 	$(call kap-fixed-metal,inlist_lowT_fa05_mb22)
 
 	cp -r $(BUILD_DIR_MODULE)/kap_input_data/wichita/AAG21_lowT_MESA_tables/*  $(BUILD_DIR_MODULE)/data/kap_data/
-	cp condtable.data $(BUILD_DIR_MODULE)/data/kap_data/
+	cp condtabl.data $(BUILD_DIR_MODULE)/data/kap_data/
 
 $(BUILD_DIR_MODULE)/AESOPUS_AGSS09.h5: $(BUILD_DIR_MODULE)/ags09/
-	cd $(BUILD_DIR_MODULE); python $(CURDIR)/AESOPUS/aesopus.py $(CURDIR)/AESOPUS/AGSS09.yaml
+	cd $(BUILD_DIR_MODULE); python $(CURDIR)/AESOPUS/aesopus.py $(CURDIR)/AESOPUS/AGSS09.yml
 
 build-aesopus: $(BUILD_DIR_MODULE)/AESOPUS_AGSS09.h5
 

kap/preprocessor/build_data_and_export

@@ -0,0 +1,11 @@
+#!/bin/bash
+
+set -e
+
+cd "$(dirname "$0")"
+
+kap_data_dir="../../build/kap-preprocessor/data/kap_data"
+
+rm -rf "$kap_data_dir"
+mkdir -p "$kap_data_dir"
+make build-data build-aesopus export-data

kap/preprocessor/src/freedman.f90

@@ -41,21 +41,21 @@ subroutine get_Freedman_fname(data_dir, Z, fname)
       character(len=*), intent(in) :: data_dir
       character(len=*), intent(out) :: fname
       integer :: iz
-      iz = floor(Z*1d5 + 0.1d0)
+      iz = nint(Z*1d6)
       select case (iz)
-      case (1000)
+      case (7452)
          fname = trim(data_dir)//'/m0.3.txt'
-      case (2000)
+      case (14760)
          fname = trim(data_dir)//'/p0.0.txt'
-      case (4000)
+      case (29024)
          fname = trim(data_dir)//'/p0.3.txt'
-      case (10000)
+      case (69840)
          fname = trim(data_dir)//'/p0.7.txt'
-      case (20000)
+      case (130292)
          fname = trim(data_dir)//'/p1.0.txt'
-      case (63000)
+      case (321456)
          fname = trim(data_dir)//'/p1.5.txt'
-      case (100000)
+      case (428844)
          fname = trim(data_dir)//'/p1.7.txt'
       case default
          write (*, *) 'get_Freedman_fname: unexpected Z value for Freedman data', Z

kap/public/kap_def.f90

@@ -618,7 +618,8 @@ subroutine kap_def_init(kap_cache_dir_in)
           kap_lowT_Xs(1:num_kap_lowT_Xs(i), i) = [ 0.00d0 ]
           num_kap_lowT_Zs(i) = 7
           kap_lowT_Zs(1:num_kap_lowT_Zs(i), i) = &
-             [ 0.01d0, 0.02d0, 0.04d0, 0.100d0, 0.200d0, 0.63d0, 1.00d0 ]
+             [ 0.007452d0, 0.014760d0, 0.029024d0, 0.069840d0, &
+               0.130292d0, 0.321456d0, 0.428844d0 ]
           num_kap_lowT_Xs_for_this_Z(1:num_kap_lowT_Zs(i), i) = num_kap_lowT_Xs(i)
        case DEFAULT
           num_kap_lowT_Xs(i) = 10