diff --git a/README.md b/README.md index 459588f..7f24b2a 100644 --- a/README.md +++ b/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 diff --git a/singularity-opac/photons/example_ascii/kap_plaw.txt b/singularity-opac/photons/example_ascii/kap_plaw.txt new file mode 100644 index 0000000..59568c1 --- /dev/null +++ b/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 diff --git a/singularity-opac/photons/example_ascii/preproc_ascii_opac.py b/singularity-opac/photons/example_ascii/preproc_ascii_opac.py new file mode 100644 index 0000000..ad9a0ed --- /dev/null +++ b/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") + diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp index d4ac76b..a6a7511 100644 --- a/singularity-opac/photons/mean_opacity_photons.hpp +++ b/singularity-opac/photons/mean_opacity_photons.hpp @@ -17,6 +17,8 @@ #define SINGULARITY_OPAC_PHOTONS_MEAN_OPACITY_PHOTONS_ #include +#include +#include #include #include @@ -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,21 +113,19 @@ 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 @@ -114,7 +133,7 @@ class MeanOpacity { 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(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(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(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(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 lkappaPlanck_; - Spiner::DataBox lkappaRosseland_; + Spiner::DataBox lkappa_; const char *filename_; + enum OpacityAveraging { + Rosseland = 0, + Planck = 1 + }; }; #undef EPS diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 21be197..22d3860 100644 --- a/test/CMakeLists.txt +++ b/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 diff --git a/test/test_mean_opacities.cpp b/test/test_mean_opacities.cpp index 932d3d1..af2d291 100644 --- a/test/test_mean_opacities.cpp +++ b/test/test_mean_opacities.cpp @@ -695,3 +695,62 @@ TEST_CASE("Mean photon scattering opacities", "[MeanPhotonS]") { opac.Finalize(); } } + +TEST_CASE("ASCII-parsed Mean photon opacities", "[MeanPhotons]") { + const std::string grayname = "kap_plaw.txt"; + + WHEN("We initialize a mean photon opacity from an ASCII table") { + + constexpr Real rho_min = 1e-14; // g/cc. + constexpr Real temp_min = 1.0; // Kelvin. + constexpr Real ross_at_min = 0.001; // cm^2/g + constexpr Real plnk_at_min = 0.1; // cm^2/g + constexpr Real rho_max = 0.7943282347241912; // g/cc. + constexpr Real temp_max = 7943282.347242886; // Kelvin. + constexpr Real ross_at_max = 0.001; // cm^2/g + constexpr Real plnk_at_max = 0.1; // cm^2/g + + photons::MeanOpacity mean_opac_host = photons::MeanOpacity(grayname); + auto mean_opac = mean_opac_host.GetOnDevice(); + + THEN("The emissivity per nu omega is consistent with the emissity per nu") { + int n_wrong_h = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View n_wrong_d("wrong"); +#else + PortableMDArray n_wrong_d(&n_wrong_h, 1); +#endif + + // get a test value at min rho-T point + Real mross = mean_opac.RosselandMeanAbsorptionCoefficient(rho_min, temp_min); + Real mplnk = mean_opac.PlanckMeanAbsorptionCoefficient(rho_min, temp_min); + + // check min rho-T point + if (FractionalDifference(rho_min * ross_at_min, mross) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(rho_min * plnk_at_min, mplnk) > EPS_TEST) { + n_wrong_d() += 1; + } + + // get a test value at max rho-T point + mross = mean_opac.RosselandMeanAbsorptionCoefficient(rho_max, temp_max); + mplnk = mean_opac.PlanckMeanAbsorptionCoefficient(rho_max, temp_max); + + // check min rho-T point + if (FractionalDifference(rho_max * ross_at_max, mross) > EPS_TEST) { + n_wrong_d() += 1; + } + if (FractionalDifference(rho_max * plnk_at_max, mplnk) > EPS_TEST) { + n_wrong_d() += 1; + } + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(n_wrong_h, n_wrong_d); +#endif + REQUIRE(n_wrong_h == 0); + } + + mean_opac.Finalize(); + } +}