Skip to content

Commit

Permalink
Add option to parse HDF5 data in ZhuTable opacity constructor.
Browse files Browse the repository at this point in the history 
+ Use filesystem header to check extention in ctor.
+ Refactor old ASCII parser to a private function (called if .txt).
+ Add ZhuTable::Save method and use to add HDF5 file.
+ Update unit test to parse HDF5 data instead.

Thanks to Jonah Miller for the suggestion to use the Save method to
create the HDF5 file.
  • Loading branch information
@RyanWollaeger
RyanWollaeger committed Oct 30, 2024
1 parent 4836c90 commit 61cf90f
Show file tree
Hide file tree
Showing 4 changed files with 112 additions and 88 deletions.
Binary file added singularity-opac/photons/zhu_data/opacitysolar09dustq3p5amax0p1new.hdf5
View file
Binary file not shown.
186 changes: 100 additions & 86 deletions singularity-opac/photons/zhu_grey_table_opacity_photons.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#define SINGULARITY_OPAC_PHOTONS_ZHU_GREY_TABLE_OPACITY_PHOTONS_

#include <fstream>
#include <filesystem>

#include <singularity-opac/base/opac_error.hpp>
#include <spiner/databox.hpp>
Expand All @@ -43,107 +44,49 @@ 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<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)
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<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::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 {
OPAC_ERROR("photons::ZhuTableOpacity: 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 += 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 {}

Expand Down Expand Up @@ -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<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)
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<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::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<pc> dist_;
Expand Down
3 changes: 3 additions & 0 deletions test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
11 changes: 9 additions & 2 deletions test/test_zhu_grey_table_opacities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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();
}
}

0 comments on commit 61cf90f

Please sign in to comment.