crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1201 from pshriwise/capi_props

Browse source 
CAPI Property Exposure and Plotting Loop Isolation
This commit is contained in:
Paul Romano 2019-03-19 13:55:14 -05:00 • committed by GitHub
commit 522b6be8eb
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
7 changed files with 271 additions and 234 deletions
Download patch file Download diff file Expand all files Collapse all files

1
include/openmc/capi.h
View file

@ -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);

115
include/openmc/plot.h
View file

@ -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<RGBColor, 2> ImageData;
typedef xt::xtensor<int32_t, 3> 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<int32_t, 3> 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<double, 3> 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<class T> 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<int, 3> pixels_; //!< Plot size in pixels
std::array<size_t, 3> pixels_; //!< Plot size in pixels
int level_; //!< Plot universe level
};
class Plot : public PlotBase
{
template<class T>
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<double>(width);
double out_pixel = (width_[1])/static_cast<double>(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)

45
openmc/capi/plot.py
View file

@ -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

308
src/plot.cpp
View file

@ -4,7 +4,8 @@
#include <fstream>
#include <sstream>
#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<double>(width);
double out_pixel = (pl.width_[1])/static_cast<double>(height);
ImageData data({width, height}, pl.not_found_);
ImageData data;
data.resize({width, height});
// generate ids for the plot
auto ids = pl.get_map<IdData>();
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<double, 3> 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<int, 3> 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<IdData>();
// select only cell ID data and flip the y-axis
xt::xtensor<int32_t, 2> 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<double>(width);
double out_pixel = (plt->width_[1])/static_cast<double>(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<IdData>();
// 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<const PlotBase*>(plot);
if (!plt) {
set_errmsg("Invalid slice pointer passed to openmc_id_map");
return OPENMC_E_INVALID_ARGUMENT;
}
auto props = plt->get_map<PropertyData>();
// write id data to array
std::copy(props.data_.begin(), props.data_.end(), data_out);
return 0;
}
} // namespace openmc

2
tests/regression_tests/plot/results_true.dat
View file

@ -1 +1 @@
30f435795c755b6fb28565a33c68e95415e5c09bca59f5549a6cbd2dd165bac9d90f5a599e12a5041cb927698abde562b7bac14175441c57b5c0dccd00544506
6b3eb36488b995b42233e6b7c0164cf6c92a1823b3c91985b97fd076f86c0aad93fbdb74e70026f5f12c61a6aee52a09588171a7fd839511ee7d9efc3952beb2

2
tests/regression_tests/plot_voxel/results_true.dat
View file

@ -1 +1 @@
76274390783bf0fd54a63e3ee141e072e797a564935a16f22a97788bc051e33fd7160dd311f0c50e45b0b498701f6735a37fabacbf9f36903381d9e18759a684
ff596b17d2cd33f964a952279b824292272ab73f57ded94812286e6bce12b5c54852d5ccd1dd137fa6312101568c06510cd89f822469d6c904a17f96259e4d03

32
tests/unit_tests/test_capi.py
View file

@ -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)