Merge pull request #54 from RyanWollaeger/zhu_ph_opac

Add ASCII opacity parser to mean_opacity_photons class.

README.md

@@ -59,6 +59,20 @@ A number of options are avaialable for compiling:
 | SINGULARITY_USE_HDF5              | ON      | Enables HDF5. Required for Spiner opacities.                                         |
 | SINGULARITY_KOKKOS_IN_TREE        | OFF     | Force cmake to use Kokkos source included in tree.                                   |
 
+### Loading ASCII Data
+
+Currently, the MeanOpacity class, defined in singularity-opac/photons/mean_opacity_photons.hpp, supports
+loading grey Rosseland and Planck opacity data in an ASCII format.  An example of this format is
+provided by singularity-opac/photons/example_ascii/kap_plaw.txt.  The 1st row of the header has the
+number of density points, NRho, then the number of temperature points, NT.  The 2nd (3rd) row of the header
+has min and max density (temperature) bounds.  These bounds are inclusive, so the opacity data in the file
+should have evaluations at these min and max values.  The rest of the ASCII file is a two-column table, where
+the 1st (2nd) column is Rosseland (Planck) opacity.  The number of rows in each column is NRhoxNT, where
+density is the slow index and temperature is the fast index (thus the row index = temperature index
++ NT x (density index), indexing from 0).  Each opacity is assumed to be evaluated on log-spaced
+density and temperature grids, where these grids are defined by NRho, NT, and the (again inclusive) min and
+max bounds the header.
+
 ## Copyright
 
 © 2021. Triad National Security, LLC. All rights reserved.  This

singularity-opac/photons/example_ascii/kap_plaw.txt

@@ -0,0 +1,19 @@
+ NRho 4 NT 4
+ rho_min 1.00000000e-14 rho_max 7.94328235e-01
+ T_min 1.00000000e+00 T_max 7.94328235e+06
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01
+1.00000000e-03 1.00000000e-01

singularity-opac/photons/example_ascii/preproc_ascii_opac.py

@@ -0,0 +1,83 @@
+# ======================================================================
+# © 2024. Triad National Security, LLC. All rights reserved.  This
+# program was produced under U.S. Government contract
+# 89233218CNA000001 for Los Alamos National Laboratory (LANL), which
+# is operated by Triad National Security, LLC for the U.S.
+# Department of Energy/National Nuclear Security Administration. All
+# rights in the program are reserved by Triad National Security, LLC,
+# and the U.S. Department of Energy/National Nuclear Security
+# Administration. The Government is granted for itself and others
+# acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+# license in this material to reproduce, prepare derivative works,
+# distribute copies to the public, perform publicly and display
+# publicly, and to permit others to do so.
+# ======================================================================
+
+#--------------------------------------------------------------------------------------------------#
+# Regrid (interpolate) Rosseland and Planck opacity data to a log-log rho-T grid.
+#--------------------------------------------------------------------------------------------------#
+import numpy as np
+from scipy.interpolate import interp2d
+import argparse
+
+#-- parse command line arguments
+parser = argparse.ArgumentParser(description="Re-grid opacity data to be log-log in density-temperature.")
+parser.add_argument("fname", type=str, help="Opacity file to re-grid.")
+parser.add_argument("--is_loglog", action="store_true", help="Avoid regridding if already log-log.")
+parser.add_argument("--fname_new", type=str, default="kappa_new.txt", help="File name for new file created.")
+args = parser.parse_args()
+
+NRho = -1
+NT = -1
+with open(args.fname, "r") as ff:
+    svec = ff.readline().split(" ")
+    NRho = int(svec[2])
+    NT = int(svec[4])
+
+print("NRho = ", NRho)
+print("NT = ", NT)
+
+opac = np.loadtxt(args.fname, skiprows=1)
+
+print("opac.shape = ", opac.shape)
+assert np.size(opac, 0) == NT * NRho, "np.size(opac, 0) != NT * NRho"
+
+#-- density, temperature grids
+Rho = np.unique(opac[:,0])
+T = np.unique(opac[:,1])
+
+assert np.size(Rho) == NRho, "np.size(Rho) != Rho"
+assert np.size(T) == NT, "np.size(T) != NT"
+
+r = np.logspace(np.log10(Rho[0]), np.log10(Rho[NRho - 1]), NRho)
+tt = np.logspace(np.log10(T[0]), np.log10(T[NT - 1]), NT)
+
+if (args.is_loglog):
+    assert np.max(abs(r - Rho) / Rho) < 1e-6, "np.max(abs(r - Rho) / Rho) >= 1e-6"
+    assert np.max(abs(tt - T) / T) < 1e-6, "np.max(abs(tt - T) / T) >= 1e-6"
+else:
+    print()
+    print("Interpolating data to log-log rho-T grid...")
+    #-- reshape to rho-T grid
+    ross_old_opac = np.reshape(opac[:,2], (NRho, NT))
+    plnk_old_opac = np.reshape(opac[:,3], (NRho, NT))
+    #-- linearly interpolate in log-log rho-T
+    ross_f2d = interp2d(np.log10(Rho), np.log10(T), ross_old_opac, kind='linear')
+    plnk_f2d = interp2d(np.log10(Rho), np.log10(T), plnk_old_opac, kind='linear')
+    ross_new_opac = ross_f2d(np.log10(r), np.log10(tt))
+    plnk_new_opac = plnk_f2d(np.log10(r), np.log10(tt))
+    #-- reset opacity array
+    for i in range(NRho):
+        for j in range(NT):
+            k = j + NT * i
+            opac[k,0] = r[i]
+            opac[k,1] = tt[j]
+            opac[k,2] = ross_new_opac[i,j]
+            opac[k,3] = plnk_new_opac[i,j]
+
+#-- save with new header containing min/max rho, T
+hdr = "NRho "+str(NRho)+" NT "+str(NT) + "\n" \
+    "rho_min {:.8e}".format(r[0])+" rho_max {:.8e}".format(r[NRho-1]) + "\n" \
+    "T_min {:.8e}".format(tt[0]) + " T_max {:.8e}".format(tt[NT-1])
+np.savetxt(args.fname_new, opac[:,2:], delimiter=" ", header=hdr, comments=" ", fmt="%.8e")
+

singularity-opac/photons/mean_opacity_photons.hpp

@@ -17,6 +17,8 @@
 #define SINGULARITY_OPAC_PHOTONS_MEAN_OPACITY_PHOTONS_
 
 #include <cmath>
+#include <filesystem>
+#include <fstream>
 
 #include <ports-of-call/portability.hpp>
 #include <singularity-opac/base/radiation_types.hpp>
@@ -57,27 +59,46 @@ class MeanOpacity {
                                      NT, lNuMin, lNuMax, NNu, lambda);
   }
 
+  // construct Planck/Rosseland DataBox from ascii file
+  MeanOpacity(const std::string &filename)  : filename_(filename.c_str()) {
+
+    // get number of density and temperature points
+    std::ifstream ff(filename.c_str());
+    const bool fexists = ff.good();
+
+    if (fexists) {
+
+      std::filesystem::path filePath(filename);
+      std::string extension = filePath.extension().string();
+
+      if (extension == ".txt") {
+        loadASCII(ff);
 #ifdef SPINER_USE_HDF
-  MeanOpacity(const std::string &filename) : filename_(filename.c_str()) {
-    herr_t status = H5_SUCCESS;
-    hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
-    status += lkappaPlanck_.loadHDF(file, SP5::MeanOpac::PlanckMeanOpacity);
-    status +=
-        lkappaRosseland_.loadHDF(file, SP5::MeanOpac::RosselandMeanOpacity);
-    status += H5Fclose(file);
+      } else if (extension == ".hdf5" || extension == ".h5" || extension == ".sp5") {
+        herr_t status = H5_SUCCESS;
+        hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+        status += lkappa_.loadHDF(file, "Rosseland & Planck Mean Opacities");
+        status += H5Fclose(file);
 
-    if (status != H5_SUCCESS) {
-      OPAC_ERROR("photons::MeanOpacity: HDF5 error\n");
+        if (status != H5_SUCCESS) {
+          OPAC_ERROR("photons::MeanOpacity: HDF5 error\n");
+        }
+#endif
+      } else {
+        OPAC_ERROR("photons::MeanOpacity: unrecognized file extension");
+      }
+
+    } else {
+      OPAC_ERROR("photons::MeanOpacity: file does not exist");
     }
   }
 
+#ifdef SPINER_USE_HDF
   void Save(const std::string &filename) const {
     herr_t status = H5_SUCCESS;
     hid_t file =
         H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
-    status += lkappaPlanck_.saveHDF(file, SP5::MeanOpac::PlanckMeanOpacity);
-    status +=
-        lkappaRosseland_.saveHDF(file, SP5::MeanOpac::RosselandMeanOpacity);
+    status += lkappa_.saveHDF(file, "Rosseland & Planck Mean Opacities");
     status += H5Fclose(file);
 
     if (status != H5_SUCCESS) {
@@ -92,29 +113,27 @@ class MeanOpacity {
 
   MeanOpacity GetOnDevice() {
     MeanOpacity other;
-    other.lkappaPlanck_ = Spiner::getOnDeviceDataBox(lkappaPlanck_);
-    other.lkappaRosseland_ = Spiner::getOnDeviceDataBox(lkappaRosseland_);
+    other.lkappa_ = Spiner::getOnDeviceDataBox(lkappa_);
     return other;
   }
 
   void Finalize() {
-    lkappaPlanck_.finalize();
-    lkappaRosseland_.finalize();
+    lkappa_.finalize();
   }
 
   PORTABLE_INLINE_FUNCTION
   Real PlanckMeanAbsorptionCoefficient(const Real rho, const Real temp) const {
     Real lRho = toLog_(rho);
     Real lT = toLog_(temp);
-    return rho * fromLog_(lkappaPlanck_.interpToReal(lRho, lT));
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Planck));
   }
 
   PORTABLE_INLINE_FUNCTION
   Real RosselandMeanAbsorptionCoefficient(const Real rho,
                                           const Real temp) const {
     Real lRho = toLog_(rho);
     Real lT = toLog_(temp);
-    return rho * fromLog_(lkappaRosseland_.interpToReal(lRho, lT));
+    return rho * fromLog_(lkappa_.interpToReal(lRho, lT, Rosseland));
   }
 
  private:
@@ -125,17 +144,16 @@ class MeanOpacity {
                         Real lNuMax, const int NNu, Real *lambda = nullptr) {
     using PC = typename Opacity::PC;
 
-    lkappaPlanck_.resize(NRho, NT);
-    lkappaPlanck_.setRange(0, lTMin, lTMax, NT);
-    lkappaPlanck_.setRange(1, lRhoMin, lRhoMax, NRho);
-    lkappaRosseland_.copyMetadata(lkappaPlanck_);
+    lkappa_.resize(NRho, NT, 2);
+    lkappa_.setRange(1, lTMin, lTMax, NT);
+    lkappa_.setRange(2, lRhoMin, lRhoMax, NRho);
 
     // Fill tables
     for (int iRho = 0; iRho < NRho; ++iRho) {
-      Real lRho = lkappaPlanck_.range(1).x(iRho);
+      Real lRho = lkappa_.range(2).x(iRho);
       Real rho = fromLog_(lRho);
       for (int iT = 0; iT < NT; ++iT) {
-        Real lT = lkappaPlanck_.range(0).x(iT);
+        Real lT = lkappa_.range(1).x(iT);
         Real T = fromLog_(lT);
         Real kappaPlanckNum = 0.;
         Real kappaPlanckDenom = 0.;
@@ -172,24 +190,112 @@ class MeanOpacity {
                                   ? singularity_opac::robust::ratio(
                                         kappaRosselandDenom, kappaRosselandNum)
                                   : 0.;
-        lkappaPlanck_(iRho, iT) = toLog_(kappaPlanck);
-        lkappaRosseland_(iRho, iT) = toLog_(kappaRosseland);
-        if (std::isnan(lkappaPlanck_(iRho, iT)) ||
-            std::isnan(lkappaRosseland_(iRho, iT))) {
+
+        lkappa_(iRho, iT, Rosseland) = toLog_(kappaRosseland);
+        lkappa_(iRho, iT, Planck) = toLog_(kappaPlanck);
+        if (std::isnan(lkappa_(iRho, iT, Rosseland)) ||
+            std::isnan(lkappa_(iRho, iT, Planck))) {
           OPAC_ERROR("photons::MeanOpacity: NAN in opacity evaluations");
         }
       }
     }
   }
+
+  // ASCII opacity file reader
+  void loadASCII(std::ifstream &ff) {
+
+    int NRho = -1;
+    int NT = -1;
+
+    // line read from file
+    std::string fline;
+
+    // read 1st line of header to get sizes
+    std::getline(ff, fline);
+
+    // tokenize fline
+    char* cfline = const_cast<char*>(fline.c_str());
+    char* fl_tok = std::strtok(cfline, " ");
+
+    // move to next token to get number of density points
+    fl_tok = std::strtok(nullptr, " ");
+    NRho = std::stoi(fl_tok);
+
+    // move to next token to get number of temperature points
+    fl_tok = std::strtok(nullptr, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    NT = std::stoi(fl_tok);
+
+    // read 2nd line of header to get min/max density
+    std::getline(ff, fline);
+    // tokenize fline
+    cfline = const_cast<char*>(fline.c_str());
+    fl_tok = std::strtok(cfline, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real RhoMin = std::stod(fl_tok);
+    fl_tok = std::strtok(nullptr, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real RhoMax = std::stod(fl_tok);
+
+    // read 3nd line of header to get min/max temperature
+    std::getline(ff, fline);
+    // tokenize fline
+    cfline = const_cast<char*>(fline.c_str());
+    fl_tok = std::strtok(cfline, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real TMin = std::stod(fl_tok);
+    fl_tok = std::strtok(nullptr, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    const Real TMax = std::stod(fl_tok);
+
+    // reseize the Planck and Rosseland databox
+    lkappa_.resize(NRho, NT, 2);
+
+    // set rho-T ranges
+    const Real lTMin = toLog_(TMin);
+    const Real lTMax = toLog_(TMax);
+    const Real lRhoMin = toLog_(RhoMin);
+    const Real lRhoMax = toLog_(RhoMax);
+    lkappa_.setRange(1, lTMin, lTMax, NT);
+    lkappa_.setRange(2, lRhoMin, lRhoMax, NRho);
+
+    // fill tables
+    for (int iRho = 0; iRho < NRho; ++iRho) {
+      const Real lRho_i = lkappa_.range(2).x(iRho);
+      for (int iT = 0; iT < NT; ++iT) {
+
+        // get new file-line
+        std::getline(ff, fline);
+        cfline = const_cast<char*>(fline.c_str());
+        fl_tok = std::strtok(cfline, " ");
+
+        // populate Rosseland opacity [cm^2/g]
+        lkappa_(iRho, iT, Rosseland) = toLog_(std::stod(fl_tok));
+
+        // populate Planck opacity [cm^2/g]
+        fl_tok = std::strtok(nullptr, " ");
+        lkappa_(iRho, iT, Planck) = toLog_(std::stod(fl_tok));
+
+        if (std::isnan(lkappa_(iRho, iT, Rosseland)) ||
+            std::isnan(lkappa_(iRho, iT, Planck))) {
+          OPAC_ERROR("photons::MeanOpacity: NAN in parsed ASCII opacity");
+        }
+      }
+    }
+  }
+
   PORTABLE_INLINE_FUNCTION Real toLog_(const Real x) const {
     return std::log10(std::abs(x) + EPS);
   }
   PORTABLE_INLINE_FUNCTION Real fromLog_(const Real lx) const {
     return std::pow(10., lx);
   }
-  Spiner::DataBox<Real> lkappaPlanck_;
-  Spiner::DataBox<Real> lkappaRosseland_;
+  Spiner::DataBox<Real> lkappa_;
   const char *filename_;
+  enum OpacityAveraging {
+    Rosseland = 0,
+    Planck = 1
+  };
 };
 
 #undef EPS

test/CMakeLists.txt

@@ -50,6 +50,11 @@ PRIVATE
   test_variant.cpp
 )
 
+# Link data file directory for Zhu table test
+file(CREATE_LINK
+    ${CMAKE_SOURCE_DIR}/singularity-opac/photons/example_ascii/kap_plaw.txt
+    ${CMAKE_CURRENT_BINARY_DIR}/kap_plaw.txt SYMBOLIC)
+
 target_link_libraries(${PROJECT_NAME}_unit_tests
 PRIVATE
   catch2_define