diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 5d07edb212..8815d243ab 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -71,6 +71,7 @@ extern "C" { int openmc_nuclide_name(int index, const char** name); int openmc_plot_geometry(); int openmc_id_map(const void* slice, int32_t* data_out); + int openmc_property_map(const void* slice, double* data_out); int openmc_reset(); int openmc_run(); void openmc_set_seed(int64_t new_seed); diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 8a9e4e74dc..2752add9cf 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -10,6 +10,8 @@ #include "hdf5.h" #include "openmc/position.h" #include "openmc/constants.h" +#include "openmc/cell.h" +#include "openmc/geometry.h" #include "openmc/particle.h" #include "openmc/xml_interface.h" @@ -52,7 +54,28 @@ struct RGBColor { }; typedef xt::xtensor ImageData; -typedef xt::xtensor IdData; + +struct IdData { + // Constructor + IdData(size_t h_res, size_t v_res); + + // Methods + void set_value(size_t y, size_t x, const Particle& p, int level); + + // Members + xt::xtensor data_; //!< 2D array of cell & material ids +}; + +struct PropertyData { + // Constructor + PropertyData(size_t h_res, size_t v_res); + + // Methods + void set_value(size_t y, size_t x, const Particle& p, int level); + + // Members + xt::xtensor data_; //!< 2D array of temperature & density data +}; enum class PlotType { slice = 1, @@ -66,24 +89,95 @@ enum class PlotBasis { }; enum class PlotColorBy { - cells = 1, - mats = 2 + cells = 0, + mats = 1 }; //=============================================================================== // Plot class //=============================================================================== -struct PlotBase { +class PlotBase { +public: + template T get_map() const; + // Members +public: Position origin_; //!< Plot origin in geometry Position width_; //!< Plot width in geometry PlotBasis basis_; //!< Plot basis (XY/XZ/YZ) - std::array pixels_; //!< Plot size in pixels + std::array pixels_; //!< Plot size in pixels int level_; //!< Plot universe level }; -class Plot : public PlotBase -{ +template +T PlotBase::get_map() const { + + size_t width = pixels_[0]; + size_t height = pixels_[1]; + + // get pixel size + double in_pixel = (width_[0])/static_cast(width); + double out_pixel = (width_[1])/static_cast(height); + + // size data array + T data(width, height); + + // setup basis indices and initial position centered on pixel + int in_i, out_i; + Position xyz = origin_; + switch(basis_) { + case PlotBasis::xy : + in_i = 0; + out_i = 1; + break; + case PlotBasis::xz : + in_i = 0; + out_i = 2; + break; + case PlotBasis::yz : + in_i = 1; + out_i = 2; + break; + } + + // set initial position + xyz[in_i] = origin_[in_i] - width_[0] / 2. + in_pixel / 2.; + xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.; + + // arbitrary direction + Direction dir = {0.7071, 0.7071, 0.0}; + + #pragma omp parallel + { + Particle p; + p.r() = xyz; + p.u() = dir; + p.coord_[0].universe = model::root_universe; + int level = level_; + int j{}; + + #pragma omp for + for (int y = 0; y < height; y++) { + p.r()[out_i] = xyz[out_i] - out_pixel * y; + for (int x = 0; x < width; x++) { + p.r()[in_i] = xyz[in_i] + in_pixel * x; + p.n_coord_ = 1; + // local variables + bool found_cell = find_cell(&p, 0); + j = p.n_coord_ - 1; + if (level >=0) {j = level + 1;} + if (found_cell) { + data.set_value(y, x, p, j); + Cell* c = model::cells[p.coord_[j].cell].get(); + } + } // inner for + } // outer for + } // omp parallel + + return data; +} + +class Plot : public PlotBase { public: // Constructor @@ -131,13 +225,6 @@ void draw_mesh_lines(Plot pl, ImageData& data); //! \param[out] image data associated with the plot object void output_ppm(Plot pl, const ImageData& data); -//! Get the rgb color for a given particle position in a plot -//! \param[in] particle with position for current pixel -//! \param[in] plot object -//! \param[out] rgb color -//! \param[out] cell or material id for particle position -void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id); - //! Initialize a voxel file //! \param[in] id of an open hdf5 file //! \param[in] dimensions of the voxel file (dx, dy, dz) diff --git a/openmc/capi/plot.py b/openmc/capi/plot.py index 78fcef912c..9ff398f5c1 100644 --- a/openmc/capi/plot.py +++ b/openmc/capi/plot.py @@ -1,7 +1,6 @@ -from ctypes import c_int, c_int32, c_double, Structure, POINTER +from ctypes import c_int, c_size_t, c_int32, c_double, Structure, POINTER from . import _dll -from .core import _DLLGlobal from .error import _error_handler import numpy as np @@ -31,7 +30,7 @@ class _Position(Structure): elif idx == 2: return self.z else: - raise IndexError("{} index is invalid for _Position".format(key)) + raise IndexError("{} index is invalid for _Position".format(idx)) def __setitem__(self, idx, val): if idx == 0: @@ -58,7 +57,7 @@ class _PlotBase(Structure): The width of the plot along the x, y, and z axes, respectively basis_ : c_int The axes basis of the plot view. - pixels_ : c_int[3] + pixels_ : c_size_t[3] The resolution of the plot in the horizontal and vertical dimensions level_ : c_int The universe level for the plot view @@ -74,9 +73,9 @@ class _PlotBase(Structure): basis : string One of {'xy', 'xz', 'yz'} indicating the horizontal and vertical axes of the plot. - h_res : float + h_res : int The horizontal resolution of the plot in pixels - v_res : float + v_res : int The vertical resolution of the plot in pixels level : int The universe level for the plot (default: -1 -> all universes shown) @@ -84,7 +83,7 @@ class _PlotBase(Structure): _fields_ = [('origin_', _Position), ('width_', _Position), ('basis_', c_int), - ('pixels_', 3*c_int), + ('pixels_', 3*c_size_t), ('level_', c_int)] def __init__(self): @@ -198,7 +197,7 @@ _dll.openmc_id_map.errcheck = _error_handler def id_map(plot): """ - Generate a 2-D map of (cell_id, material_id). Used for in-memory image + Generate a 2-D map of cell and material IDs. Used for in-memory image generation. Parameters @@ -215,6 +214,32 @@ def id_map(plot): """ img_data = np.zeros((plot.v_res, plot.h_res, 2), dtype=np.dtype('int32')) - _dll.openmc_id_map(POINTER(_PlotBase)(plot), - img_data.ctypes.data_as(POINTER(c_int32))) + _dll.openmc_id_map(plot, img_data.ctypes.data_as(POINTER(c_int32))) return img_data + + +_dll.openmc_property_map.argtypes = [POINTER(_PlotBase), POINTER(c_double)] +_dll.openmc_property_map.restype = c_int +_dll.openmc_property_map.errcheck = _error_handler + + +def property_map(plot): + """ + Generate a 2-D map of cell temperatures and material densities. Used for + in-memory image generation. + + Parameters + ---------- + plot : openmc.capi.plot._PlotBase + Object describing the slice of the model to be generated + + Returns + ------- + property_map : numpy.ndarray + A NumPy array with shape (vertical pixels, horizontal pixels, 2) of + OpenMC property ids with dtype float + + """ + prop_data = np.zeros((plot.v_res, plot.h_res, 2)) + _dll.openmc_property_map(plot, prop_data.ctypes.data_as(POINTER(c_double))) + return prop_data diff --git a/src/plot.cpp b/src/plot.cpp index 729b2db486..02fb42d76b 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -4,7 +4,8 @@ #include #include -#include "openmc/cell.h" +#include "xtensor/xview.hpp" + #include "openmc/constants.h" #include "openmc/file_utils.h" #include "openmc/geometry.h" @@ -28,7 +29,39 @@ namespace openmc { const RGBColor WHITE {255, 255, 255}; constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level -constexpr int NOT_FOUND {-1}; +constexpr int32_t NOT_FOUND {-2}; + +IdData::IdData(size_t h_res, size_t v_res) + : data_({v_res, h_res, 2}, NOT_FOUND) +{ } + +void +IdData::set_value(size_t y, size_t x, const Particle& p, int level) { + Cell* c = model::cells[p.coord_[level].cell].get(); + data_(y,x,0) = c->id_; + if (p.material_ == MATERIAL_VOID) { + data_(y,x,1) = MATERIAL_VOID; + return; + } else if (c->type_ != FILL_UNIVERSE) { + Material* m = model::materials[p.material_].get(); + data_(y,x,1) = m->id_; + } +} + +PropertyData::PropertyData(size_t h_res, size_t v_res) + : data_({v_res, h_res, 2}, NOT_FOUND) +{ } + +void +PropertyData::set_value(size_t y, size_t x, const Particle& p, int level) { + Cell* c = model::cells[p.coord_[level].cell].get(); + data_(y,x,0) = (p.sqrtkT_ * p.sqrtkT_) / K_BOLTZMANN; + if (c->type_ != FILL_UNIVERSE && p.material_ != MATERIAL_VOID) { + Material* m = model::materials[p.material_].get(); + data_(y,x,1) = m->density_gpcc_; + } +} + //============================================================================== // Global variables //============================================================================== @@ -100,60 +133,28 @@ void create_ppm(Plot pl) size_t width = pl.pixels_[0]; size_t height = pl.pixels_[1]; - double in_pixel = (pl.width_[0])/static_cast(width); - double out_pixel = (pl.width_[1])/static_cast(height); + ImageData data({width, height}, pl.not_found_); - ImageData data; - data.resize({width, height}); + // generate ids for the plot + auto ids = pl.get_map(); - int in_i, out_i; - Position r; - switch(pl.basis_) { - case PlotBasis::xy : - in_i = 0; - out_i = 1; - r.x = pl.origin_[0] - pl.width_[0] / 2.; - r.y = pl.origin_[1] + pl.width_[1] / 2.; - r.z = pl.origin_[2]; - break; - case PlotBasis::xz : - in_i = 0; - out_i = 2; - r.x = pl.origin_[0] - pl.width_[0] / 2.; - r.y = pl.origin_[1]; - r.z = pl.origin_[2] + pl.width_[1] / 2.; - break; - case PlotBasis::yz : - in_i = 1; - out_i = 2; - r.x = pl.origin_[0]; - r.y = pl.origin_[1] - pl.width_[0] / 2.; - r.z = pl.origin_[2] + pl.width_[1] / 2.; - break; - } - - Direction u {0.7071, 0.7071, 0.0}; - - #pragma omp parallel - { - Particle p; - p.r() = r; - p.u() = u; - p.coord_[0].universe = model::root_universe; - - #pragma omp for - for (int y = 0; y < height; y++) { - p.r()[out_i] = r[out_i] - out_pixel * y; - for (int x = 0; x < width; x++) { - // local variables - RGBColor rgb; - int id; - p.r()[in_i] = r[in_i] + in_pixel * x; - position_rgb(p, pl, rgb, id); - data(x,y) = rgb; - } - } - } + // assign colors + for (size_t y = 0; y < height; y++) { + for (size_t x = 0; x < width; x++) { + auto id = ids.data_(y, x, pl.color_by_); + // no setting needed if not found + if (id == NOT_FOUND) { continue; } + if (PlotColorBy::cells == pl.color_by_) { + data(x,y) = pl.colors_[model::cell_map[id]]; + } else if (PlotColorBy::mats == pl.color_by_) { + if (id == MATERIAL_VOID) { + data(x,y) = WHITE; + continue; + } + data(x,y) = pl.colors_[model::material_map[id]]; + } // color_by if-else + } // x for loop + } // y for loop // draw mesh lines if present if (pl.index_meshlines_mesh_ >= 0) {draw_mesh_lines(pl, data);} @@ -637,53 +638,6 @@ Plot::Plot(pugi::xml_node plot_node) set_mask(plot_node); } // End Plot constructor -//============================================================================== -// POSITION_RGB computes the red/green/blue values for a given plot with the -// current particle's position -//============================================================================== - - -void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) -{ - p.n_coord_ = 1; - - bool found_cell = find_cell(&p, 0); - - int j = p.n_coord_ - 1; - - if (settings::check_overlaps) {check_cell_overlap(&p);} - - // Set coordinate level if specified - if (pl.level_ >= 0) {j = pl.level_ + 1;} - - if (!found_cell) { - // If no cell, revert to default color - rgb = pl.not_found_; - id = NOT_FOUND; - } else { - if (PlotColorBy::mats == pl.color_by_) { - // Assign color based on material - const auto& c = model::cells[p.coord_[j].cell]; - if (c->type_ == FILL_UNIVERSE) { - // If we stopped on a middle universe level, treat as if not found - rgb = pl.not_found_; - id = NOT_FOUND; - } else if (p.material_ == MATERIAL_VOID) { - // By default, color void cells white - rgb = WHITE; - id = NOT_FOUND; - } else { - rgb = pl.colors_[p.material_]; - id = model::materials[p.material_]->id_; - } - } else if (PlotColorBy::cells == pl.color_by_) { - // Assign color based on cell - rgb = pl.colors_[p.coord_[j].cell]; - id = model::cells[p.coord_[j].cell]->id_; - } - } // endif found_cell -} - //============================================================================== // OUTPUT_PPM writes out a previously generated image to a PPM file //============================================================================== @@ -815,18 +769,17 @@ void draw_mesh_lines(Plot pl, ImageData& data) // CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D // geometry visualization. It works the same way as create_ppm by dragging a // particle across the geometry for the specified number of voxels. The first 3 -// int(4)'s in the binary are the number of x, y, and z voxels. The next 3 -// real(8)'s are the widths of the voxels in the x, y, and z directions. The -// next 3 real(8)'s are the x, y, and z coordinates of the lower left -// point. Finally the binary is filled with entries of four int(4)'s each. Each -// 'row' in the binary contains four int(4)'s: 3 for x,y,z position and 1 for +// int's in the binary are the number of x, y, and z voxels. The next 3 +// double's are the widths of the voxels in the x, y, and z directions. The +// next 3 double's are the x, y, and z coordinates of the lower left +// point. Finally the binary is filled with entries of four int's each. Each +// 'row' in the binary contains four int's: 3 for x,y,z position and 1 for // cell or material id. For 1 million voxels this produces a file of // approximately 15MB. // ============================================================================= void create_voxel(Plot pl) { - // compute voxel widths in each direction std::array vox; vox[0] = pl.width_[0]/(double)pl.pixels_[0]; @@ -836,13 +789,6 @@ void create_voxel(Plot pl) // initial particle position Position ll = pl.origin_ - pl.width_ / 2.; - // allocate and initialize particle - Direction u {0.7071, 0.7071, 0.0}; - Particle p; - p.r() = ll; - p.u() = u; - p.coord_[0].universe = model::root_universe; - // Open binary plot file for writing std::ofstream of; std::string fname = std::string(pl.path_plot_); @@ -861,7 +807,9 @@ void create_voxel(Plot pl) // Write current date and time write_attribute(file_id, "date_and_time", time_stamp().c_str()); hsize_t three = 3; - write_attribute(file_id, "num_voxels", pl.pixels_); + std::array pixels; + std::copy(pl.pixels_.begin(), pl.pixels_.end(), pixels.begin()); + write_attribute(file_id, "num_voxels", pixels); write_attribute(file_id, "voxel_width", vox); write_attribute(file_id, "lower_left", ll); @@ -874,43 +822,33 @@ void create_voxel(Plot pl) hid_t dspace, dset, memspace; voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace); - // move to center of voxels - ll.x += vox[0] / 2.; - ll.y += vox[1] / 2.; - ll.z += vox[2] / 2.; - - int data[pl.pixels_[1]][pl.pixels_[0]]; + PlotBase pltbase; + pltbase.width_ = pl.width_; + pltbase.origin_ = pl.origin_; + pltbase.basis_ = PlotBasis::xy; + pltbase.pixels_ = pl.pixels_; + pltbase.level_ = -1; // all universes for voxel files ProgressBar pb; - - RGBColor rgb; - int id; for (int z = 0; z < pl.pixels_[2]; z++) { + // update progress bar pb.set_value(100.*(double)z/(double)(pl.pixels_[2]-1)); - for (int y = 0; y < pl.pixels_[1]; y++) { - for (int x = 0; x < pl.pixels_[0]; x++) { - // get voxel color - position_rgb(p, pl, rgb, id); - // write to plot data - data[y][x] = id; - // advance particle in x direction - p.r().x += vox[0]; - } - // advance particle in y direction - p.r().y += vox[1]; - p.r().x = ll[0]; - } - // advance particle in z direction - p.r().z += vox[2]; - p.r().y = ll[1]; - p.r().x = ll[0]; + + // update z coordinate + pltbase.origin_.z = ll.z + z * vox[2]; + + // generate ids using plotbase + IdData ids = pltbase.get_map(); + + // select only cell ID data and flip the y-axis + xt::xtensor data1 = xt::flip(xt::view(ids.data_, xt::all(), xt::all(), 0), 0); + // Write to HDF5 dataset - voxel_write_slice(z, dspace, dset, memspace, &(data[0])); + voxel_write_slice(z, dspace, dset, memspace, &(data1(0,0))); } voxel_finalize(dspace, dset, memspace); file_close(file_id); - } void @@ -963,75 +901,29 @@ extern "C" int openmc_id_map(const void* plot, int32_t* data_out) return OPENMC_E_INVALID_ARGUMENT; } - size_t width = plt->pixels_[0]; - size_t height = plt->pixels_[1]; - - // get pixel size - double in_pixel = (plt->width_[0])/static_cast(width); - double out_pixel = (plt->width_[1])/static_cast(height); - - // size data array - IdData data({height, width, 2}, NOT_FOUND); - - // setup basis indices and initial position centered on pixel - int in_i, out_i; - Position xyz = plt->origin_; - switch(plt->basis_) { - case PlotBasis::xy : - in_i = 0; - out_i = 1; - break; - case PlotBasis::xz : - in_i = 0; - out_i = 2; - break; - case PlotBasis::yz : - in_i = 1; - out_i = 2; - break; - } - - // set initial position - xyz[in_i] = plt->origin_[in_i] - plt->width_[0] / 2. + in_pixel / 2.; - xyz[out_i] = plt->origin_[out_i] + plt->width_[1] / 2. - out_pixel / 2.; - - // arbitrary direction - Direction dir = {0.5, 0.5, 0.5}; - - #pragma omp parallel - { - Particle p; - p.r() = xyz; - p.u() = dir; - p.coord_[0].universe = model::root_universe; - int level = plt->level_; - int j{}; - - #pragma omp for - for (int y = 0; y < height; y++) { - p.r()[out_i] = xyz[out_i] - out_pixel * y; - for (int x = 0; x < width; x++) { - p.r()[in_i] = xyz[in_i] + in_pixel * x; - p.n_coord_ = 1; - // local variables - bool found_cell = find_cell(&p, 0); - j = p.n_coord_ - 1; - if (level >=0) {j = level + 1;} - if (found_cell) { - Cell* c = model::cells[p.coord_[j].cell].get(); - data(y,x,0) = c->id_; - if (c->type_ != FILL_UNIVERSE && p.material_ != MATERIAL_VOID) { - data(y,x,1) = model::materials[p.material_]->id_; - } - } - } // inner for - } // outer for - } // omp parallel + auto ids = plt->get_map(); // write id data to array - std::copy(data.begin(), data.end(), data_out); + std::copy(ids.data_.begin(), ids.data_.end(), data_out); return 0; } +extern "C" int openmc_property_map(const void* plot, double* data_out) { + + auto plt = reinterpret_cast(plot); + if (!plt) { + set_errmsg("Invalid slice pointer passed to openmc_id_map"); + return OPENMC_E_INVALID_ARGUMENT; + } + + auto props = plt->get_map(); + + // write id data to array + std::copy(props.data_.begin(), props.data_.end(), data_out); + + return 0; +} + + } // namespace openmc diff --git a/tests/regression_tests/plot/results_true.dat b/tests/regression_tests/plot/results_true.dat index 7f8fcc6a82..18251e8d0c 100644 --- a/tests/regression_tests/plot/results_true.dat +++ b/tests/regression_tests/plot/results_true.dat @@ -1 +1 @@ -30f435795c755b6fb28565a33c68e95415e5c09bca59f5549a6cbd2dd165bac9d90f5a599e12a5041cb927698abde562b7bac14175441c57b5c0dccd00544506 \ No newline at end of file +6b3eb36488b995b42233e6b7c0164cf6c92a1823b3c91985b97fd076f86c0aad93fbdb74e70026f5f12c61a6aee52a09588171a7fd839511ee7d9efc3952beb2 \ No newline at end of file diff --git a/tests/regression_tests/plot_voxel/results_true.dat b/tests/regression_tests/plot_voxel/results_true.dat index 8016217847..52222000f6 100644 --- a/tests/regression_tests/plot_voxel/results_true.dat +++ b/tests/regression_tests/plot_voxel/results_true.dat @@ -1 +1 @@ -76274390783bf0fd54a63e3ee141e072e797a564935a16f22a97788bc051e33fd7160dd311f0c50e45b0b498701f6735a37fabacbf9f36903381d9e18759a684 \ No newline at end of file +ff596b17d2cd33f964a952279b824292272ab73f57ded94812286e6bce12b5c54852d5ccd1dd137fa6312101568c06510cd89f822469d6c904a17f96259e4d03 \ No newline at end of file diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 97d18531bd..be88072d9f 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -419,3 +419,35 @@ def test_id_map(capi_init): ids = openmc.capi.plot.id_map(s) assert np.array_equal(expected_ids, ids) + +def test_property_map(capi_init): + expected_properties = np.array( + [[(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)], + [ (293.6, 6.55), (293.6, 10.29769), (293.6, 6.55)], + [(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)]], dtype='float') + + # create a plot object + s = openmc.capi.plot._PlotBase() + s.width = 1.26 + s.height = 1.26 + s.v_res = 3 + s.h_res = 3 + s.origin = (0.0, 0.0, 0.0) + s.basis = 'xy' + s.level = -1 + + properties = openmc.capi.plot.property_map(s) + assert np.allclose(expected_properties, properties, atol=1e-04) + + +def test_position(capi_init): + + pos = openmc.capi.plot._Position(1.0, 2.0, 3.0) + + assert tuple(pos) == (1.0, 2.0, 3.0) + + pos[0] = 1.3 + pos[1] = 2.3 + pos[2] = 3.3 + + assert tuple(pos) == (1.3, 2.3, 3.3)