Put Zhu ASCII parsing into mean_opacity_photons (BRR suggestion).

+ Remove Zhu grey table class and unit test.
+ Replace Zhu data table with toy power-law data table.
+ Generalize file-based ctor to take ASCII or HDF5 file.
+ Add sub-test for parsing toy power-law data into MeanOpacity.

singularity-opac/photons/example_ascii/create_mean_opac_ascii.py

@@ -0,0 +1,56 @@
+#--------------------------------------------------------------------------------------------------#
+# This script defaults to using the temperature and density bounds for opacity data found
+# in Zhaohuan Zhu et al (2021), "Global 3D Radiation Hydrodynamic Simulations of Proto-Jupiter’s
+# Convective EnvelopeConvective Envelope". Otherwise, it merely generates a power law for Rosseland
+# and Planck-averaged opacities, and produces a table in an ASCII format equivalent to that of Zhu.
+#--------------------------------------------------------------------------------------------------#
+import numpy as np
+import argparse
+
+#-- parse command line arguments
+parser = argparse.ArgumentParser(description="Create power-law grey opacity ASCII table (for testing).")
+parser.add_argument("fname", type=str, help="base opacity file name to create (.txt is appended).")
+parser.add_argument("--Nrho", type=int, default=4, help="Number of density points (assumes log space)")
+parser.add_argument("--rho", type=float, nargs=2, default=[1e-14, 0.7943282347241912],
+                    help="min, max densities [g/cc]")
+parser.add_argument("--NT", type=int, default=4, help="Number of temperature points (assumes log space)")
+parser.add_argument("--T", type=float, nargs=2, default=[1.0, 7943282.347242886],
+                    help="min, max temperatures [K]")
+parser.add_argument("--rho_ref", type=float, default=1e-6, help="Reference density [g/cc]")
+parser.add_argument("--rho_exp", type=float, default=0.0, help="Exponent of density / ref. density")
+parser.add_argument("--T_ref", type=float, default=11604, help="Reference density [K]")
+parser.add_argument("--T_exp", type=float, default=0.0, help="Exponent of temperature / ref. temperature")
+parser.add_argument("--Pkap_coef", type=float, default=0.01, help="Planck opacity coefficient [cm^2/g]")
+parser.add_argument("--Rkap_coef", type=float, default=0.01, help="Rosseland opacity coefficient [cm^2/g]")
+args = parser.parse_args()
+
+#-- implement simple power law
+def power_law_kappa(r, t):
+    '''power_law_kappa (cgs): r=density (1st arg), t=temperature (2nd arg)'''
+    kap_coef = [args.Rkap_coef, args.Pkap_coef]
+    kap_pwrs = (r / args.rho_ref)**(args.rho_exp) * (t / args.T_ref)**(args.T_exp)
+    return [kapc * kap_pwrs for kapc in kap_coef]
+
+#-- set rho and T arrays
+rho = np.logspace(np.log10(args.rho[0]), np.log10(args.rho[1]), args.Nrho)
+T = np.logspace(np.log10(args.T[0]), np.log10(args.T[1]), args.NT)
+
+#-- initialize the opacity array - columns: density, temperature, Rosseland, Planck
+Nrow = args.Nrho * args.NT
+opac = np.zeros((Nrow, 4))
+
+#-- populate opacity array
+for i in range(args.Nrho):
+    for j in range(args.NT):
+        k = j + args.NT * i
+        opac[k,0] = rho[i]
+        opac[k,1] = T[j]
+        opac[k,2:] = power_law_kappa(rho[i], T[j])
+
+#-- save opacity array
+np.savetxt(args.fname+'.txt', opac, delimiter=" ",
+           header="rho "+str(args.Nrho)+" T "+str(args.NT), comments=" ", fmt="%.8e")
+
+#--------------------------------------------------------------------------------------------------#
+# End of create_mean_opac_ascii.py
+#--------------------------------------------------------------------------------------------------#

singularity-opac/photons/example_ascii/kap_plaw.txt

@@ -0,0 +1,17 @@
+ rho 4 T 4
+1.00000000e-14 1.00000000e+00 1.00000000e-03 1.00000000e-01
+1.00000000e-14 1.99526231e+02 1.00000000e-03 1.00000000e-01
+1.00000000e-14 3.98107171e+04 1.00000000e-03 1.00000000e-01
+1.00000000e-14 7.94328235e+06 1.00000000e-03 1.00000000e-01
+4.29866235e-10 1.00000000e+00 1.00000000e-03 1.00000000e-01
+4.29866235e-10 1.99526231e+02 1.00000000e-03 1.00000000e-01
+4.29866235e-10 3.98107171e+04 1.00000000e-03 1.00000000e-01
+4.29866235e-10 7.94328235e+06 1.00000000e-03 1.00000000e-01
+1.84784980e-05 1.00000000e+00 1.00000000e-03 1.00000000e-01
+1.84784980e-05 1.99526231e+02 1.00000000e-03 1.00000000e-01
+1.84784980e-05 3.98107171e+04 1.00000000e-03 1.00000000e-01
+1.84784980e-05 7.94328235e+06 1.00000000e-03 1.00000000e-01
+7.94328235e-01 1.00000000e+00 1.00000000e-03 1.00000000e-01
+7.94328235e-01 1.99526231e+02 1.00000000e-03 1.00000000e-01
+7.94328235e-01 3.98107171e+04 1.00000000e-03 1.00000000e-01
+7.94328235e-01 7.94328235e+06 1.00000000e-03 1.00000000e-01

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,20 +59,44 @@ class MeanOpacity {
                                      NT, lNuMin, lNuMax, NNu, lambda);
   }
 
+  // construct Planck/Rosseland DataBox from Zhu-formatted 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") {
+        // assume this is one of the original Zhu et al (2021) ASCII files
+        loadZhuASCII(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 += lkappaPlanck_.loadHDF(file, SP5::MeanOpac::PlanckMeanOpacity);
+        status +=
+          lkappaRosseland_.loadHDF(file, SP5::MeanOpac::RosselandMeanOpacity);
+        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 =
@@ -181,6 +207,86 @@ class MeanOpacity {
       }
     }
   }
+
+  // if .txt file, assume it is the original Zhu dust opacity file
+  void loadZhuASCII(std::ifstream &ff) {
+
+    int NRho = -1;
+    int NT = -1;
+
+    // line read from file
+    std::string fline;
+
+    // read 1-line 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 density points
+    fl_tok = std::strtok(nullptr, " ");
+    fl_tok = std::strtok(nullptr, " ");
+    NT = std::stoi(fl_tok);
+
+    // reseize the Planck and Rosseland databoxes (number of types of opac=2)
+    lkappaRosseland_.resize(NRho, NT);
+    lkappaPlanck_.resize(NRho, NT);
+
+    // set rho-T rankges (Zhu tables are uniform in log-log rho-T space)
+    const Real lTMin = toLog_(1.0);
+    const Real lTMax = toLog_(7943282.347242886);
+    const Real lRhoMin = toLog_(1.0e-14);
+    const Real lRhoMax = toLog_(0.7943282347241912);
+    lkappaRosseland_.setRange(1, lTMin, lTMax, NT);
+    lkappaRosseland_.setRange(2, lRhoMin, lRhoMax, NRho);
+    lkappaPlanck_.setRange(1, lTMin, lTMax, NT);
+    lkappaPlanck_.setRange(2, lRhoMin, lRhoMax, NRho);
+
+    // fill tables
+    for (int iRho = 0; iRho < NRho; ++iRho) {
+      const Real lRho_i = lkappaRosseland_.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, " ");
+
+        // check for consistent density [g/cm^3] on table row
+        const Real Rho = std::stod(fl_tok);
+        if (std::abs(Rho - fromLog_(lRho_i)) > 1e-6 * std::abs(Rho)) {
+          OPAC_ERROR("photons::MeanOpacity: invalid rho");
+        }
+
+        // check for consistent temperature [K] on table row
+        const Real lT_i = lkappaRosseland_.range(1).x(iT);
+        fl_tok = std::strtok(nullptr, " ");
+        const Real T = std::stod(fl_tok);
+        if (std::abs(T - fromLog_(lT_i)) > 1e-6 * std::abs(T)) {
+          OPAC_ERROR("photons::MeanOpacity: invalid T");
+        }
+
+        // populate Rosseland opacity [cm^2/g]
+        fl_tok = std::strtok(nullptr, " ");
+        lkappaRosseland_(iRho, iT) = toLog_(std::stod(fl_tok));
+
+        // populate Planck opacity [cm^2/g]
+        fl_tok = std::strtok(nullptr, " ");
+        lkappaPlanck_(iRho, iT) = toLog_(std::stod(fl_tok));
+
+        if (std::isnan(lkappaPlanck_(iRho, iT)) ||
+            std::isnan(lkappaRosseland_(iRho, iT))) {
+          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);
   }