Skip to content

Add gray data to converted IPCRESS file #80

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
alexrlongne wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Add gray data to converted IPCRESS file #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2c3f409
Add gray opacity to the ipcress2spiner conversion tool
alexrlongne Aug 7, 2026
b26ac96
+ Add converted file read test
alexrlongne Aug 12, 2026
bced6d9
Remove unused fields from sp5
alexrlongne Aug 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion ipcress2spiner/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ add_executable(ipcress2spiner
main.cpp
)


target_include_directories(ipcress2spiner
PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
Expand All @@ -46,6 +45,13 @@ target_compile_definitions(ipcress2spiner PRIVATE
IPCRESS2SPINER_VERSION=\"${PROJECT_VERSION}\"
)
install(TARGETS ipcress2spiner DESTINATION ${CMAKE_INSTALL_BINDIR})

# Tests
enable_testing()
add_subdirectory(test)

else()
message("GANDOLF was not found, not building ipcress2spiner")


endif()
2 changes: 0 additions & 2 deletions ipcress2spiner/generate_files.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,10 @@
#include <utils/spiner/spiner/sp5.hpp>
#include <utils/spiner/spiner/databox.hpp>
#include <utils/spiner/spiner/interpolation.hpp>
#include <utils/spiner/spiner/singularity_opac_sp5.hpp>

herr_t saveMaterial(hid_t loc, hid_t matGroup, const int matid,
const std::string &sMatid,
const std::string &sp5_field_name,
const std::vector<double> &group_bounds,
Spiner::DataBox<double> &opacity) {


Expand Down
103 changes: 73 additions & 30 deletions ipcress2spiner/interpolate_opacity.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,61 +18,50 @@
#define _IPCRESS2SPINER_INTERPOLATE_HPP_

inline double interpolate_mg_opacity_data(const std::vector<double> &T_grid, const std::vector<double> &rho_grid, const std::vector<double> &group_bounds, const std::vector<double> &op_data,
double target_T, double target_rho, double target_hnu, const double log_interpolation) {
double target_T, double target_rho, double target_hnu) {

size_t n_T = T_grid.size();
size_t n_rho = rho_grid.size();
size_t n_groups = group_bounds.size()-1;

// always use log interpolation between temperatures, densities, and opacity table values
auto T_L = std::distance(T_grid.begin(), std::lower_bound(T_grid.begin(), T_grid.end(), target_T));
T_L = (target_T < T_grid[T_L]) ? T_L-1 : T_L;
T_L = (T_L == (T_grid.size()-1)) ? T_L-1: T_L;
auto T_R = T_L + 1;
// fraction of T_2 to use
double T_frac = -1.0; // defaulted to invalid value
if (log_interpolation) {
auto log_T_L = log(T_grid[T_L]);
auto log_T_R = log(T_grid[T_R]);
auto log_T_target = log(target_T);
T_frac = (log_T_target - log_T_L) / (log_T_R - log_T_L);
}
else {
T_frac = (target_T - T_grid[T_L]) / (T_grid[T_R] - T_grid[T_L]);
}
// fraction of T_R opacity to use
double log_T_L = log(T_grid[T_L]);
double log_T_R = log(T_grid[T_R]);
double log_T_target = log(target_T);
double T_frac = (log_T_target - log_T_L) / (log_T_R - log_T_L);

auto rho_L = std::distance(rho_grid.begin(), std::lower_bound(rho_grid.begin(), rho_grid.end(), target_rho));
rho_L = (target_rho < rho_grid[rho_L]) ? rho_L-1 : rho_L;
rho_L = (rho_L == (rho_grid.size()-1)) ? rho_L-1: rho_L;
auto rho_R = rho_L + 1;
// fraction of rho_2 to use
double rho_frac = -1.0; // defaulted to invalid value
if (log_interpolation) {
auto log_rho_L = log(rho_grid[rho_L]);
auto log_rho_R = log(rho_grid[rho_R]);
auto log_rho_target = log(target_rho);
rho_frac = (log_rho_target - log_rho_L) / (log_rho_R - log_rho_L);
}
else {
rho_frac = (target_rho - rho_grid[rho_L]) / (rho_grid[rho_R] - rho_grid[rho_L]);
}
// fraction of rho_R opacity to use
double log_rho_L = log(rho_grid[rho_L]);
double log_rho_R = log(rho_grid[rho_R]);
double log_rho_target = log(target_rho);
double rho_frac = (log_rho_target - log_rho_L) / (log_rho_R - log_rho_L);

auto group = std::distance(group_bounds.begin(), std::lower_bound(group_bounds.begin(), group_bounds.end(), target_hnu));
group = (target_hnu < group_bounds[group]) ? group-1 : group;
group = (group == (group_bounds.size()-1)) ? group-1 : group;

// get the adjacent rows of the opacity index, looks like op_[T]_[rho]_[hnu] with L indicating
// left or index n and "R" indicating right index (n+1)
auto op_L_L = op_data[ T_L * (n_rho*n_groups) + rho_L*(n_groups) + group];
auto op_L_R = op_data[ T_L * (n_rho*n_groups) + rho_R*(n_groups) + group];
auto op_R_L = op_data[ T_R * (n_rho*n_groups) + rho_L*(n_groups) + group];
auto op_R_R = op_data[ T_R * (n_rho*n_groups) + rho_R*(n_groups) + group];
double op_L_L = log(op_data[ T_L * (n_rho*n_groups) + rho_L*(n_groups) + group]);
double op_L_R = log(op_data[ T_L * (n_rho*n_groups) + rho_R*(n_groups) + group]);
double op_R_L = log(op_data[ T_R * (n_rho*n_groups) + rho_L*(n_groups) + group]);
double op_R_R = log(op_data[ T_R * (n_rho*n_groups) + rho_R*(n_groups) + group]);

// reduce in temperature dimension
auto op_L =(1.0-T_frac) * op_L_L + T_frac*op_R_L;
auto op_R =(1.0-T_frac) * op_L_R + T_frac*op_R_R;
auto op_L =op_L_L + T_frac*(op_R_L - op_L_L);
auto op_R =op_L_R + T_frac*(op_R_R - op_L_R);

// reduce in density dimension and return
double interp_opac = (1.0-rho_frac) * op_L + rho_frac*op_R;
double interp_opac = exp(op_L + rho_frac*(op_R -op_L));
if (false) {
std::cout<<"Targets--T: "<<target_T<<" rho: "<<target_rho<<" hnu: "<<target_hnu<<" T frac:" <<T_frac<<" rho frac: "<<rho_frac<<std::endl;
std::cout<<"Interpolation box: "<<std::endl;
Expand All @@ -87,4 +76,58 @@ inline double interpolate_mg_opacity_data(const std::vector<double> &T_grid, con
return interp_opac;
}

inline double interpolate_gray_opacity_data(const std::vector<double> &T_grid, const std::vector<double> &rho_grid, const std::vector<double> &op_data,
double target_T, double target_rho) {

size_t n_T = T_grid.size();
size_t n_rho = rho_grid.size();

// always use log interpolation between temperatures, densities, and opacity table values
auto T_L = std::distance(T_grid.begin(), std::lower_bound(T_grid.begin(), T_grid.end(), target_T));
T_L = (target_T < T_grid[T_L]) ? T_L-1 : T_L;
T_L = (T_L == (T_grid.size()-1)) ? T_L-1: T_L;
auto T_R = T_L + 1;
// fraction of T_R opacity to use
double log_T_L = log(T_grid[T_L]);
double log_T_R = log(T_grid[T_R]);
double log_T_target = log(target_T);
double T_frac = (log_T_target - log_T_L) / (log_T_R - log_T_L);

auto rho_L = std::distance(rho_grid.begin(), std::lower_bound(rho_grid.begin(), rho_grid.end(), target_rho));
rho_L = (target_rho < rho_grid[rho_L]) ? rho_L-1 : rho_L;
rho_L = (rho_L == (rho_grid.size()-1)) ? rho_L-1: rho_L;
auto rho_R = rho_L + 1;
// fraction of rho_R opacity to use
double log_rho_L = log(rho_grid[rho_L]);
double log_rho_R = log(rho_grid[rho_R]);
double log_rho_target = log(target_rho);
double rho_frac = (log_rho_target - log_rho_L) / (log_rho_R - log_rho_L);

// get the adjacent rows of the opacity index, looks like op_[T]_[rho]_[hnu] with L indicating
// left or index n and "R" indicating right index (n+1)
double op_L_L = log(op_data[ T_L * n_rho + rho_L]);
double op_L_R = log(op_data[ T_L * n_rho + rho_R]);
double op_R_L = log(op_data[ T_R * n_rho + rho_L]);
double op_R_R = log(op_data[ T_R * n_rho + rho_R]);

// reduce in temperature dimension
double op_L =op_L_L + T_frac*(op_R_L - op_L_L);
double op_R =op_L_R + T_frac*(op_R_R - op_L_R);

// reduce in density dimension and return
double interp_opac = exp(op_L + rho_frac*(op_R - op_L));
if (false) {
std::cout<<"Targets--T: "<<target_T<<" rho: "<<target_rho<< "T frac:" <<T_frac<<" rho frac: "<<rho_frac<<std::endl;
std::cout<<"Interpolation box: "<<std::endl;
std::cout<<op_L_R<<"--------"<<op_R_R<<" |"<<std::endl;
std::cout<<"| |"<<std::endl;
std::cout<<"| |"<<std::endl;
std::cout<<"| "<<interp_opac<<" |"<<std::endl;
std::cout<<"| |"<<std::endl;
std::cout<<"| |"<<std::endl;
std::cout<<op_L_L<<"--------"<<op_R_L<<std::endl;
}
return interp_opac;
}

#endif // _IPCRESS2SPINER_INTERPOLATE_HPP_
80 changes: 50 additions & 30 deletions ipcress2spiner/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,13 +43,10 @@
#include "make_spiner_databox.hpp"

// ARL: My to-do list
// \TODO Thread log interpolation points and log interpolation mode options through command line
// \TODO Thread equal log spacing of rho and T points through command line
// \TODO Thread verbosity option through command-line options to print more info
// \TODO Add gray opacties
// \TODO Add ionization state?
// \TODO Add ionization state? Leave until singularity-eos and singularity-opac are merged
// \TODO clang-format
// \TODO Add tests
// \TODO Compare interpolation to draco's routines
// \TODO Delete char fields from new[]

// Function prototypes for Gandolf C function equivalents, prepend "C" to prevent name mangling
Expand Down Expand Up @@ -201,49 +198,72 @@ int main(int argc, char *argv[]) {
std::vector<double> temperature_points(nt, 0.0);
std::vector<double> group_bounds(nhnu, 0.0);
std::vector<double> density_points(nrho, 0.0);
// these post values are potentially reset by calls to gandolf?
int post_nt = nt;
int post_nrho = nrho;
int post_nhnu = nhnu;
int post_nmg = nmg;
int post_ngray = ngray;

std::string ramg_keyword = "ramg";
std::string rsmg_keyword = "rsmg";
std::string rtmg_keyword = "rtmg";
std::string pmg_keyword = "pmg";
std::string ragray_keyword = "ragray";
std::string rgray_keyword = "rgray";
std::string pgray_keyword = "pgray";

std::array<std::string, 3> mg_opac_keywords{ramg_keyword, rsmg_keyword, pmg_keyword};
std::array<std::string, 3> mg_fields{ SP5::Fields::ramg, SP5::Fields::rsmg, SP5::Fields::pmg};
std::array<std::string, 7> mg_opac_keywords{ramg_keyword, rsmg_keyword, rtmg_keyword, pmg_keyword,
ragray_keyword, rgray_keyword, pgray_keyword};
std::array<std::string, 7> mg_fields{ SP5::Fields::ramg, SP5::Fields::rsmg, SP5::Fields::rtmg, SP5::Fields::pmg,
SP5::Fields::ragray, SP5::Fields::rgray, SP5::Fields::pgray};

// TODO Thread these variables through the input
bool log_T_rho_hnu = true; // form a new grid from max and min evenly spaced in log10
bool log_interpolation = true; // interpolate logarithmically

for (size_t i=0;i<mg_opac_keywords.size();++i) {
auto key = mg_opac_keywords[i];
auto sp5_field_name = mg_fields[i];
std::cout<<"Interpolating and building databox for "<<mat_ID<<" and "<<key<<std::endl;
std::vector<double> opacity_data(nmg, 0.0);
c_ggetmg(filename_char, &mat_ID, key.data(),
temperature_points.data(), &nt, &post_nt,
density_points.data(), &nrho, &post_nrho,
group_bounds.data(), &nhnu, &post_nhnu,
opacity_data.data(), &nmg, &post_nmg, &ier);

auto [spiner_opacity_databox, new_group_bounds] = build_opacity_spiner_databox(
temperature_points, density_points, group_bounds, opacity_data, log_T_rho_hnu,
log_interpolation);

// only save the group bounds once for this material
if (i==0) {
const hsize_t dimensions[] = {static_cast<hsize_t>(new_group_bounds.size())};
std::string dataset_name("group bounds");
hid_t dataspace_id = H5Screate_simple(1, dimensions, nullptr);
hid_t dataset_id = H5Dcreate2(matGroup, dataset_name.c_str(), H5T_IEEE_F64LE, dataspace_id,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
const herr_t status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
new_group_bounds.data());
// multigroup opacities
if(key == "ramg" || key == "rsmg" || key=="rtmg" || key == "pmg") {
std::vector<double> opacity_data(nmg, 0.0);
c_ggetmg(filename_char, &mat_ID, key.data(),
temperature_points.data(), &nt, &post_nt,
density_points.data(), &nrho, &post_nrho,
group_bounds.data(), &nhnu, &post_nhnu,
opacity_data.data(), &nmg, &post_nmg, &ier);

auto [spiner_opacity_databox, new_group_bounds] = build_multigroup_opacity_spiner_databox(
temperature_points, density_points, group_bounds, opacity_data, log_T_rho_hnu);

// only save the group bounds once for this material
if (i==0) {
const hsize_t dimensions[] = {static_cast<hsize_t>(new_group_bounds.size())};
std::string dataset_name("group bounds");
hid_t dataspace_id = H5Screate_simple(1, dimensions, nullptr);
hid_t dataset_id = H5Dcreate2(matGroup, dataset_name.c_str(), H5T_IEEE_F64LE, dataspace_id,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
const herr_t status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
new_group_bounds.data());
}
std::cout<<"Saving multigroup databox for "<<mat_ID<<" and "<<key<<std::endl;
saveMaterial(file_loc, matGroup, mat_ID, sMatID, sp5_field_name, spiner_opacity_databox);
}
// gray opacites
else {
std::vector<double> opacity_data(ngray, 0.0);
c_ggetgray(filename_char, &mat_ID, key.data(),
temperature_points.data(), &nt, &post_nt,
density_points.data(), &nrho, &post_nrho,
opacity_data.data(), &ngray, &post_ngray, &ier);

auto spiner_opacity_databox = build_gray_opacity_spiner_databox(
temperature_points, density_points, opacity_data, log_T_rho_hnu);

std::cout<<"Saving gray databox for "<<mat_ID<<" and "<<key<<std::endl;
saveMaterial(file_loc, matGroup, mat_ID, sMatID, sp5_field_name, spiner_opacity_databox);
}
std::cout<<"Saving databox for "<<mat_ID<<" and "<<key<<std::endl;
saveMaterial(file_loc, matGroup, mat_ID, sMatID, sp5_field_name, new_group_bounds, spiner_opacity_databox);
}
status += H5Gclose(matGroup);
}
Expand Down
52 changes: 45 additions & 7 deletions ipcress2spiner/make_spiner_databox.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
#ifndef _IPCRESS2SPINER_MAKE_SPINER_DATABOX_HPP_
#define _IPCRESS2SPINER_MAKE_SPINER_DATABOX_HPP_

inline std::pair<Spiner::DataBox<double>, std::vector<double>> build_opacity_spiner_databox(const std::vector<double> &temperature_points, const std::vector<double> &density_points, const std::vector<double> &group_bounds, const std::vector<double> &ramg, const bool log_T_rho_hnu, const bool log_interpolation, const bool verbose_mode = false) {
inline std::pair<Spiner::DataBox<double>, std::vector<double>> build_multigroup_opacity_spiner_databox(const std::vector<double> &temperature_points, const std::vector<double> &density_points, const std::vector<double> &group_bounds, const std::vector<double> &full_multigroup_data, const bool log_T_rho_hnu, const bool verbose_mode = false) {

const auto n_groups = group_bounds.size() -1;
const double temperature_start = (log_T_rho_hnu) ? log10(temperature_points.front()) : temperature_points.front();
Expand All @@ -29,11 +29,15 @@ inline std::pair<Spiner::DataBox<double>, std::vector<double>> build_opacity_sp
const double hnu_midpoints_start = (log_T_rho_hnu) ? log10( 0.5*(group_bounds[0]+group_bounds[1])) : 0.5*(group_bounds[0] + group_bounds[1]);
const double hnu_midpoints_end = (log_T_rho_hnu) ? log10(0.5*(group_bounds[n_groups]+group_bounds[n_groups-1])) : 0.5*(group_bounds[n_groups]+group_bounds[n_groups-1]);

Spiner::RegularGrid1D<double> new_temperature(temperature_start, temperature_end, temperature_points.size());
Spiner::RegularGrid1D<double> new_density(density_start, density_end, density_points.size());
// ARL: Double number of interpolation points in density and temperature until non-uniform grids
// are supported

size_t n_points_multiplier = 2;
Spiner::RegularGrid1D<double> new_temperature(temperature_start, temperature_end, n_points_multiplier*temperature_points.size());
Spiner::RegularGrid1D<double> new_density(density_start, density_end, n_points_multiplier*density_points.size());
Spiner::RegularGrid1D<double> new_group_bounds(hnu_start, hnu_end, group_bounds.size());
Spiner::RegularGrid1D<double> new_group_centers(hnu_midpoints_start, hnu_midpoints_end, group_bounds.size()-1);
Spiner::DataBox<double> opacity_databox(temperature_points.size(), density_points.size(), n_groups);
Spiner::DataBox<double> opacity_databox(n_groups, n_points_multiplier*density_points.size(), n_points_multiplier*temperature_points.size());
opacity_databox.setRange(0, new_temperature);
opacity_databox.setRange(1, new_density);
opacity_databox.setRange(2, new_group_centers);
Expand All @@ -42,7 +46,6 @@ inline std::pair<Spiner::DataBox<double>, std::vector<double>> build_opacity_sp
std::cout<<"T rho hnu opac(sq_cm/g)"<<std::endl;
}

size_t flat_index=0;
for (size_t k=0; k<new_temperature.nPoints();++k) {
for (size_t j=0; j<new_density.nPoints();++j) {
for (size_t i=0; i<(new_group_bounds.nPoints()-1);++i) {
Expand All @@ -51,8 +54,8 @@ inline std::pair<Spiner::DataBox<double>, std::vector<double>> build_opacity_sp
new_group_center = 0.5*(std::pow(10, new_group_bounds.x(i)) + std::pow(10, new_group_bounds.x(i+1)));
double target_T = (log_T_rho_hnu) ? std::pow(10, new_temperature.x(k)) : new_temperature.x(k);
double target_rho = (log_T_rho_hnu) ? std::pow(10, new_density.x(j)) : new_density.x(j);
double interp_op = interpolate_mg_opacity_data(temperature_points, density_points, group_bounds, ramg, target_T, target_rho, new_group_center, log_interpolation);
opacity_databox(k,j,i) = interp_op;
double interp_op = interpolate_mg_opacity_data(temperature_points, density_points, group_bounds, full_multigroup_data, target_T, target_rho, new_group_center);
opacity_databox(i,j, k) = interp_op;
if(verbose_mode) {
std::cout<<target_T <<" "<<target_rho<<" "<<new_group_center<<" "<<opacity_databox(k,j,i)<<std::endl;
}
Expand All @@ -69,5 +72,40 @@ inline std::pair<Spiner::DataBox<double>, std::vector<double>> build_opacity_sp
return {opacity_databox, vector_new_group_bounds};
}

inline Spiner::DataBox<double> build_gray_opacity_spiner_databox(const std::vector<double> &temperature_points, const std::vector<double> &density_points, const std::vector<double> &full_gray_data, const bool log_T_rho_hnu, const bool verbose_mode = false) {

const double temperature_start = (log_T_rho_hnu) ? log10(temperature_points.front()) : temperature_points.front();
const double temperature_end = (log_T_rho_hnu) ? log10(temperature_points.back()) : temperature_points.back();
const double density_start = (log_T_rho_hnu) ? log10(density_points.front()) : density_points.front();
const double density_end = (log_T_rho_hnu) ? log10(density_points.back()) : density_points.back();

// ARL: For gray, use 10x the number of interpolation points in density and temperature as it's
// more difficult to match when the tables are coarse
size_t n_points_multiplier = 10;
Spiner::RegularGrid1D<double> new_temperature(temperature_start, temperature_end, n_points_multiplier*temperature_points.size());
Spiner::RegularGrid1D<double> new_density(density_start, density_end, n_points_multiplier*density_points.size());
Spiner::DataBox<double> opacity_databox(n_points_multiplier*density_points.size(), n_points_multiplier*temperature_points.size());
opacity_databox.setRange(0, new_temperature);
opacity_databox.setRange(1, new_density);

if(verbose_mode) {
std::cout<<"T rho opac(sq_cm/g)"<<std::endl;
}

for (size_t k=0; k<new_temperature.nPoints();++k) {
for (size_t j=0; j<new_density.nPoints();++j) {
double target_T = (log_T_rho_hnu) ? std::pow(10, new_temperature.x(k)) : new_temperature.x(k);
double target_rho = (log_T_rho_hnu) ? std::pow(10, new_density.x(j)) : new_density.x(j);
double interp_op = interpolate_gray_opacity_data(temperature_points, density_points, full_gray_data, target_T, target_rho);
opacity_databox(j,k) = interp_op;
if(verbose_mode) {
std::cout<<target_T <<" "<<target_rho<<" "<<opacity_databox(j,k)<<std::endl;
}
}
}

return opacity_databox;
}

#endif // _IPCRESS2SPINER_MAKE_SPINER_DATABOX_HPP_

Loading
Loading