diff --git a/ipcress2spiner/CMakeLists.txt b/ipcress2spiner/CMakeLists.txt index fc4ca8a..00edbdd 100644 --- a/ipcress2spiner/CMakeLists.txt +++ b/ipcress2spiner/CMakeLists.txt @@ -24,7 +24,6 @@ add_executable(ipcress2spiner main.cpp ) - target_include_directories(ipcress2spiner PUBLIC $ @@ -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() diff --git a/ipcress2spiner/generate_files.hpp b/ipcress2spiner/generate_files.hpp index db0f826..27d384c 100644 --- a/ipcress2spiner/generate_files.hpp +++ b/ipcress2spiner/generate_files.hpp @@ -22,12 +22,10 @@ #include #include #include -#include 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 &group_bounds, Spiner::DataBox &opacity) { diff --git a/ipcress2spiner/interpolate_opacity.hpp b/ipcress2spiner/interpolate_opacity.hpp index 4f0e9b2..549c04f 100644 --- a/ipcress2spiner/interpolate_opacity.hpp +++ b/ipcress2spiner/interpolate_opacity.hpp @@ -18,43 +18,32 @@ #define _IPCRESS2SPINER_INTERPOLATE_HPP_ inline double interpolate_mg_opacity_data(const std::vector &T_grid, const std::vector &rho_grid, const std::vector &group_bounds, const std::vector &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; @@ -62,17 +51,17 @@ inline double interpolate_mg_opacity_data(const std::vector &T_grid, con // 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: "< &T_grid, con return interp_opac; } +inline double interpolate_gray_opacity_data(const std::vector &T_grid, const std::vector &rho_grid, const std::vector &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: "< temperature_points(nt, 0.0); std::vector group_bounds(nhnu, 0.0); std::vector 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 mg_opac_keywords{ramg_keyword, rsmg_keyword, pmg_keyword}; - std::array mg_fields{ SP5::Fields::ramg, SP5::Fields::rsmg, SP5::Fields::pmg}; + std::array mg_opac_keywords{ramg_keyword, rsmg_keyword, rtmg_keyword, pmg_keyword, + ragray_keyword, rgray_keyword, pgray_keyword}; + std::array 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 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(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 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(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 "< 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 "<, std::vector> build_opacity_spiner_databox(const std::vector &temperature_points, const std::vector &density_points, const std::vector &group_bounds, const std::vector &ramg, const bool log_T_rho_hnu, const bool log_interpolation, const bool verbose_mode = false) { +inline std::pair, std::vector> build_multigroup_opacity_spiner_databox(const std::vector &temperature_points, const std::vector &density_points, const std::vector &group_bounds, const std::vector &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(); @@ -29,11 +29,15 @@ inline std::pair, std::vector> 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 new_temperature(temperature_start, temperature_end, temperature_points.size()); - Spiner::RegularGrid1D 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 new_temperature(temperature_start, temperature_end, n_points_multiplier*temperature_points.size()); + Spiner::RegularGrid1D new_density(density_start, density_end, n_points_multiplier*density_points.size()); Spiner::RegularGrid1D new_group_bounds(hnu_start, hnu_end, group_bounds.size()); Spiner::RegularGrid1D new_group_centers(hnu_midpoints_start, hnu_midpoints_end, group_bounds.size()-1); - Spiner::DataBox opacity_databox(temperature_points.size(), density_points.size(), n_groups); + Spiner::DataBox 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); @@ -42,7 +46,6 @@ inline std::pair, std::vector> build_opacity_sp std::cout<<"T rho hnu opac(sq_cm/g)"<, std::vector> 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<, std::vector> build_opacity_sp return {opacity_databox, vector_new_group_bounds}; } +inline Spiner::DataBox build_gray_opacity_spiner_databox(const std::vector &temperature_points, const std::vector &density_points, const std::vector &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 new_temperature(temperature_start, temperature_end, n_points_multiplier*temperature_points.size()); + Spiner::RegularGrid1D new_density(density_start, density_end, n_points_multiplier*density_points.size()); + Spiner::DataBox 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)"< +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#ifndef SPINER_USE_HDF +#error "HDF5 must be enabled" +#endif + +#include +#include +#include +#include +#include + +int main() { + + const char * filename = "test_file.h5"; + const char * group_name = "/10001"; + //const char * dataset_name = "/10001/plank total gray opacity"; + herr_t status = H5_SUCCESS; + + hid_t file; + hid_t dataset; + hid_t dataspace; + + // Open HDF5 file read-only + file = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT); + + // Open dataset + dataset = H5Gopen2(file, group_name, H5P_DEFAULT); + + // test gray databox + { + std::cout<<"--- Testing gray databox interpolation--- "< gray_databox; + gray_databox.loadHDF(dataset, SP5::Fields::rgray); + + // make sure gray databox has the expected values + if (gray_databox.rank() != 2) { + std::cout<<"Rank check fails! Rank is: "< T_points{0.1, 0.6, .99, 1.0, 1.01, 1.4, 3.0, 4.2, 6.2, 9.0, 10.0}; + // gold values from Draco's ipcress interpreter + std::vector gold_values{4.000100e-01, 4.758533e-01, 4.949213e-01, 0.5, 6.848793e-01, + 1.518317e+00, 1.879426e+01, 5.707131e+01, 2.064356e+02, + 7.064486e+02, 1.000307e+03}; + + std::cout<<"T gold interpolated_value relative_diff"< tolerance) { + std::cout<<"Value doesn't match to "< mg_databox; + mg_databox.loadHDF(dataset, SP5::Fields::rsmg); + + // make sure mg databox has the expected values + if (mg_databox.rank() != 3) { + std::cout<<"Rank check fails! Rank is: "< hnu_points{0.02, 0.05, 0.085, 0.2, 0.5, 0.85, 2.0 , 5.0 , 8.5, 20.0, 50.0, 85.0}; + // gold values from Draco's ipcress interpreter + std::vector gold_values{0.3999584351001695, 0.39989004841500875, 0.399799302343056, + 0.39958473935770206, 0.39890319754138404, 0.39800204558786706, + 0.39588583373187775, 0.3892961262035159, 0.38087997690498276, + 0.3623440968644833, 0.3137366760288401, 0.2663122909516335}; + std::cout<<"hnu gold interpolated_value relative_diff"< tolerance) { + std::cout<<"Value doesn't match to "<