diff --git a/singularity-opac/photons/zhu_data/opacitysolar09dustq3p5amax0p1new.hdf5 b/singularity-opac/photons/zhu_data/opacitysolar09dustq3p5amax0p1new.hdf5 new file mode 100644 index 0000000..4d73f67 Binary files /dev/null and b/singularity-opac/photons/zhu_data/opacitysolar09dustq3p5amax0p1new.hdf5 differ diff --git a/singularity-opac/photons/zhu_grey_table_opacity_photons.hpp b/singularity-opac/photons/zhu_grey_table_opacity_photons.hpp index 1ee9a32..3afc81f 100644 --- a/singularity-opac/photons/zhu_grey_table_opacity_photons.hpp +++ b/singularity-opac/photons/zhu_grey_table_opacity_photons.hpp @@ -17,6 +17,7 @@ #define SINGULARITY_OPAC_PHOTONS_ZHU_GREY_TABLE_OPACITY_PHOTONS_ #include +#include #include #include @@ -43,100 +44,28 @@ class ZhuTableOpacity { ZhuTableOpacity() = default; // construct Planck/Rosseland DataBox from Zhu ascii file - ZhuTableOpacity(const int amax, const bool use_planck_absorb = false) + ZhuTableOpacity(const std::string filename, const bool use_planck_absorb = false) : opac_type_(use_planck_absorb) { - // TODO: replace this parsing with DataBox HDF loadHDF (requires pre-process of ASCII data) - - // select file based on dust grain size index - std::string filename = "no_file"; - switch(amax) { - case -1: - filename = "opacitysolar09dustq3p5amax0p1new.txt"; - break; - case -3: - filename = "opacitysolar09dustq3p5amax0p001new.txt"; - break; - case 0: - filename = "opacitysolar09dustq3p5amax1new.txt"; - break; - case 10: - filename = "opacitysolar09dustq3p5amax10new.txt"; - break; - default: - OPAC_ERROR("photons::ZhuTableOpacity: unrecognized amax ctor value."); - } - // get number of density and temperature points std::ifstream ff(filename.c_str()); const bool fexists = ff.good(); - int NRho = -1; - int NT = -1; if (fexists) { - // line read from file - std::string fline; - - // read 1-line 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 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) - kappa_.resize(NRho, NT, 2); - - // 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); - kappa_.setRange(1, lTMin, lTMax, NT); - kappa_.setRange(2, lRhoMin, lRhoMax, NRho); - - // fill tables - for (int iRho = 0; iRho < NRho; ++iRho) { - const Real lRho_i = kappa_.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, " "); - - // 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::ZhuTableOpacity: invalid rho"); - } - - // check for consistent temperature [K] on table row - const Real lT_i = kappa_.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::ZhuTableOpacity: invalid T"); - } - - // populate Rosseland opacity [cm^2/g] - fl_tok = std::strtok(nullptr, " "); - kappa_(iRho, iT, 0) = std::stod(fl_tok); - - // populate Planck opacity [cm^2/g] - fl_tok = std::strtok(nullptr, " "); - kappa_(iRho, iT, 1) = std::stod(fl_tok); - } + 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 + } else if (extension == ".hdf5" || extension == ".h5") { + hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + kappa_.loadHDF(file, "Zhu Table"); +#endif + } else { + OPAC_ERROR("photons::ZhuTableOpacity: unrecognized file extension"); } } else { @@ -144,6 +73,20 @@ class ZhuTableOpacity { } } +#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 += kappa_.saveHDF(file, "Zhu Table"); + status += H5Fclose(file); + + if (status != H5_SUCCESS) { + OPAC_ERROR("photons::ZhuTableOpacity: HDF5 error\n"); + } + } +#endif + ZhuTableOpacity GetOnDevice() { return *this; } inline void Finalize() noexcept {} @@ -283,6 +226,77 @@ class ZhuTableOpacity { return std::pow(10., lx); } + // 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(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) + kappa_.resize(NRho, NT, 2); + + // 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); + kappa_.setRange(1, lTMin, lTMax, NT); + kappa_.setRange(2, lRhoMin, lRhoMax, NRho); + + // fill tables + for (int iRho = 0; iRho < NRho; ++iRho) { + const Real lRho_i = kappa_.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, " "); + + // 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::ZhuTableOpacity: invalid rho"); + } + + // check for consistent temperature [K] on table row + const Real lT_i = kappa_.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::ZhuTableOpacity: invalid T"); + } + + // populate Rosseland opacity [cm^2/g] + fl_tok = std::strtok(nullptr, " "); + kappa_(iRho, iT, 0) = std::stod(fl_tok); + + // populate Planck opacity [cm^2/g] + fl_tok = std::strtok(nullptr, " "); + kappa_(iRho, iT, 1) = std::stod(fl_tok); + } + } + } + // WARNING: which only one of the two possibilities can be used for now int opac_type_; PlanckDistribution dist_; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c8c0463..2860324 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -55,6 +55,9 @@ PRIVATE file(CREATE_LINK ${CMAKE_SOURCE_DIR}/singularity-opac/photons/zhu_data/opacitysolar09dustq3p5amax0p1new.txt ${CMAKE_CURRENT_BINARY_DIR}/opacitysolar09dustq3p5amax0p1new.txt SYMBOLIC) +file(CREATE_LINK + ${CMAKE_SOURCE_DIR}/singularity-opac/photons/zhu_data/opacitysolar09dustq3p5amax0p1new.hdf5 + ${CMAKE_CURRENT_BINARY_DIR}/opacitysolar09dustq3p5amax0p1new.hdf5 SYMBOLIC) target_link_libraries(${PROJECT_NAME}_unit_tests PRIVATE diff --git a/test/test_zhu_grey_table_opacities.cpp b/test/test_zhu_grey_table_opacities.cpp index 9ea55ad..93704ff 100644 --- a/test/test_zhu_grey_table_opacities.cpp +++ b/test/test_zhu_grey_table_opacities.cpp @@ -39,7 +39,8 @@ TEST_CASE("Zhu table photon opacities", "[GrayPhotons]") { WHEN("We initialize a gray Zhu tabular photon opacity") { // values are directly copied from parsed table for a_grain=-1 - constexpr int a_grain = -1; + const std::string fbase = "opacitysolar09dustq3p5amax0p1new"; + const std::string fname = fbase + ".hdf5"; constexpr Real rho_min = 1e-14; // g/cc. constexpr Real temp_min = 1.0; // Kelvin. constexpr Real ross_at_min = 0.0030420184427588414; // cm^2/g @@ -49,7 +50,8 @@ TEST_CASE("Zhu table photon opacities", "[GrayPhotons]") { constexpr Real nu = 3e9; // Hz. UHF microwave - photons::Opacity opac_host = photons::ZhuTable(a_grain); + // not all opacities have a save method, so not using variant for host here + photons::ZhuTable opac_host = photons::ZhuTable(fname); photons::Opacity opac = opac_host.GetOnDevice(); // Check constants from mean opacity @@ -73,6 +75,11 @@ TEST_CASE("Zhu table photon opacities", "[GrayPhotons]") { REQUIRE(n_wrong == 0); } + // uncomment this save method to creat hdf5 file from ZhuTable DataBox + // const std::string hdfname = fbase + ".hdf5"; + // opac_host.Save(hdfname); + + // finished opac.Finalize(); } }