crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1107 from pshriwise/plot_to_cpp

Browse source 
Plotting to C++
This commit is contained in:
Sterling Harper 2018-11-05 12:35:59 -05:00 • committed by GitHub
commit 56dc03c19c
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
20 changed files with 1321 additions and 1181 deletions
Download patch file Download diff file Expand all files Collapse all files

4
CMakeLists.txt
View file

@ -344,9 +344,6 @@ add_library(libopenmc SHARED
src/photon_physics.F90
src/physics_common.F90
src/physics.F90
src/plot.F90
src/plot_header.F90
src/progress_header.F90
src/pugixml/pugixml_f.F90
src/random_lcg.F90
src/reaction_header.F90
@ -414,6 +411,7 @@ add_library(libopenmc SHARED
src/physics_mg.cpp
src/plot.cpp
src/position.cpp
src/progress_bar.cpp
src/pugixml/pugixml_c.cpp
src/random_lcg.cpp
src/reaction.cpp

16
include/openmc/output.h
View file

@ -4,6 +4,7 @@
#ifndef OPENMC_OUTPUT_H
#define OPENMC_OUTPUT_H
#include <string>
namespace openmc {
@ -16,6 +17,21 @@ namespace openmc {
void header(const char* msg, int level);
//==============================================================================
//! Retrieve a time stamp
//!
//! \return current time stamp (format: "yyyy-mm-dd hh:mm:ss")
//==============================================================================
std::string time_stamp();
//==============================================================================
//! Display plot information
//==============================================================================
extern "C" void print_plot();
//==============================================================================
//! Display information regarding cell overlap checking.
//==============================================================================

172
include/openmc/plot.h
View file

@ -1,15 +1,181 @@
#ifndef OPENMC_PLOT_H
#define OPENMC_PLOT_H
#include <unordered_map>
#include <sstream>
#include "xtensor/xarray.hpp"
#include "hdf5.h"
#include "openmc/position.h"
#include "openmc/constants.h"
#include "openmc/particle.h"
#include "openmc/xml_interface.h"
namespace openmc {
extern "C" void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
//===============================================================================
// Global variables
//===============================================================================
extern int PLOT_LEVEL_LOWEST; //!< lower bound on plot universe level
extern std::unordered_map<int, int> plot_map; //!< map of plot ids to index
extern "C" int32_t n_plots; //!< number of plots in openmc run
class Plot;
extern std::vector<Plot> plots; //!< Plot instance container
//===============================================================================
// RGBColor holds color information for plotted objects
//===============================================================================
struct RGBColor {
//Constructors
RGBColor() : red(0), green(0), blue(0) { };
RGBColor(const int v[3]) : red(v[0]), green(v[1]), blue(v[2]) { };
RGBColor(int r, int g, int b) : red(r), green(g), blue(b) { };
RGBColor(const std::vector<int> &v) {
if (v.size() != 3) {
throw std::out_of_range("Incorrect vector size for RGBColor.");
}
red = v[0];
green = v[1];
blue = v[2];
}
// Members
uint8_t red, green, blue;
};
typedef xt::xtensor<RGBColor, 2> ImageData;
enum class PlotType {
slice = 1,
voxel = 2
};
enum class PlotBasis {
xy = 1,
xz = 2,
yz = 3
};
enum class PlotColorBy {
cells = 1,
mats = 2
};
//===============================================================================
// Plot class
//===============================================================================
class Plot
{
public:
// Constructor
Plot(pugi::xml_node plot);
// Methods
private:
void set_id(pugi::xml_node plot_node);
void set_type(pugi::xml_node plot_node);
void set_output_path(pugi::xml_node plot_node);
void set_bg_color(pugi::xml_node plot_node);
void set_basis(pugi::xml_node plot_node);
void set_origin(pugi::xml_node plot_node);
void set_width(pugi::xml_node plot_node);
void set_universe(pugi::xml_node plot_node);
void set_default_colors(pugi::xml_node plot_node);
void set_user_colors(pugi::xml_node plot_node);
void set_meshlines(pugi::xml_node plot_node);
void set_mask(pugi::xml_node plot_node);
// Members
public:
int id_; //!< Plot ID
PlotType type_; //!< Plot type (Slice/Voxel)
PlotColorBy color_by_; //!< Plot coloring (cell/material)
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
int meshlines_width_; //!< Width of lines added to the plot
int level_; //!< Plot universe level
int index_meshlines_mesh_; //!< Index of the mesh to draw on the plot
RGBColor meshlines_color_; //!< Color of meshlines on the plot
RGBColor not_found_; //!< Plot background color
std::vector<RGBColor> colors_; //!< Plot colors
std::string path_plot_; //!< Plot output filename
};
//===============================================================================
// Non-member functions
//===============================================================================
//! Add mesh lines to image data of a plot object
//! \param[in] plot object
//! \param[out] image data associated with the plot object
void draw_mesh_lines(Plot pl, ImageData& data);
//! Write a ppm image to file using a plot object's image data
//! \param[in] plot object
//! \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)
//! \param[out] dataspace pointer to voxel data
//! \param[out] dataset pointer to voxesl data
//! \param[out] pointer to memory space of voxel data
void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
hid_t* dset, hid_t* memspace);
extern "C" void voxel_write_slice(int x, hid_t dspace, hid_t dset,
//! Write a section of the voxel data to hdf5
//! \param[in] voxel slice
//! \param[out] dataspace pointer to voxel data
//! \param[out] dataset pointer to voxesl data
//! \param[out] pointer to data to write
void voxel_write_slice(int x, hid_t dspace, hid_t dset,
hid_t memspace, void* buf);
extern "C" void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
//! Close voxel file entities
//! \param[in] data space to close
//! \param[in] dataset to close
//! \param[in] memory space to close
void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
//===============================================================================
// External functions
//===============================================================================
//! Read plots from plots.xml node
//! \param[in] plot node of plots.xml
extern "C" void read_plots(pugi::xml_node* plot_node);
//! Create a ppm image for a plot object
//! \param[in] plot object
void create_ppm(Plot pl);
//! Create an hdf5 voxel file for a plot object
//! \param[in] plot object
void create_voxel(Plot pl);
//! Create a randomly generated RGB color
//! \return RGBColor with random value
RGBColor random_color();
} // namespace openmc
#endif // OPENMC_PLOT_H

22
include/openmc/progress_bar.h Normal file
View file

@ -0,0 +1,22 @@
#ifndef OPENMC_PROGRESSBAR_H
#define OPENMC_PROGRESSBAR_H
#include <string>
class ProgressBar {
public:
// Constructor
ProgressBar();
void set_value(double val);
private:
std::string bar;
char bar_old[72] = "???% | |";
};
#endif // OPENMC_PROGRESSBAR_H

3
include/openmc/settings.h
View file

@ -61,7 +61,8 @@ extern std::string path_statepoint; //!< path to a statepoint file
extern "C" int32_t index_entropy_mesh; //!< Index of entropy mesh in global mesh array
extern "C" int32_t index_ufs_mesh; //!< Index of UFS mesh in global mesh array
extern "C" int32_t index_cmfd_mesh; //!< Index of CMFD mesh in global mesh array
extern "C" int32_t n_batches; //!< number of (inactive+active) batches
extern "C" int32_t n_inactive; //!< number of inactive batches
extern "C" int32_t gen_per_batch; //!< number of generations per batch

2
include/openmc/string_functions.h
View file

@ -14,5 +14,7 @@ char* strtrim(char* c_str);
void to_lower(std::string& str);
int word_count(std::string const& str);
} // namespace openmc
#endif // STRING_FUNCTIONS_H

1
include/openmc/xml_interface.h
View file

@ -10,7 +10,6 @@
#include "xtensor/xarray.hpp"
#include "xtensor/xadapt.hpp"
namespace openmc {
inline bool

2
src/api.F90
View file

@ -63,7 +63,6 @@ contains
use geometry_header
use material_header
use photon_header
use plot_header
use sab_header
use settings
use simulation_header
@ -85,7 +84,6 @@ contains
call free_memory_geometry()
call free_memory_surfaces()
call free_memory_material()
call free_memory_plot()
call free_memory_volume()
call free_memory_simulation()
call free_memory_nuclide()

4
src/cmfd_execute.cpp
View file

@ -9,10 +9,10 @@
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/simulation.h"
#include "openmc/settings.h"
namespace openmc {
extern "C" int index_cmfd_mesh;
extern "C" void
cmfd_populate_sourcecounts(int n_energy, const double* energies,
@ -24,7 +24,7 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies,
openmc_source_bank(&source_bank, &n);
// Get source counts in each mesh bin / energy bin
auto& m = meshes.at(index_cmfd_mesh);
auto& m = meshes.at(settings::index_cmfd_mesh);
xt::xarray<double> counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside);
// Copy data from the xarray into the source counts array

447
src/input_xml.F90
View file

@ -23,9 +23,8 @@ module input_xml
use mgxs_interface
use nuclide_header
use multipole_header
use output, only: title, header, print_plot
use output, only: title, header
use photon_header
use plot_header
use random_lcg, only: prn
use surface_header
use set_header, only: SetChar, SetInt
@ -97,6 +96,14 @@ module input_xml
integer(C_INT) :: n
end function maximum_levels
subroutine read_plots(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_plots
subroutine print_plot() bind(C)
end subroutine print_plot
subroutine set_particle_energy_bounds(particle, E_min, E_max) bind(C)
import C_INT, C_DOUBLE
integer(C_INT), value :: particle
@ -2030,26 +2037,10 @@ contains
subroutine read_plots_xml()
integer :: i, j
integer :: n_cols, col_id, n_comp, n_masks, n_meshlines
integer :: meshid
integer(C_INT) :: err, idx
integer, allocatable :: iarray(:)
logical :: file_exists ! does plots.xml file exist?
character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
character(MAX_LINE_LEN) :: temp_str
character(MAX_WORD_LEN) :: meshtype
type(ObjectPlot), pointer :: pl => null()
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_plot
type(XMLNode) :: node_col
type(XMLNode) :: node_mask
type(XMLNode) :: node_meshlines
type(XMLNode), allocatable :: node_plot_list(:)
type(XMLNode), allocatable :: node_col_list(:)
type(XMLNode), allocatable :: node_mask_list(:)
type(XMLNode), allocatable :: node_meshline_list(:)
! Check if plots.xml exists
filename = trim(path_input) // "plots.xml"
@ -2066,425 +2057,7 @@ contains
call doc % load_file(filename)
root = doc % document_element()
! Get list pointer to XML <plot>
call get_node_list(root, "plot", node_plot_list)
! Allocate plots array
n_plots = size(node_plot_list)
allocate(plots(n_plots))
READ_PLOTS: do i = 1, n_plots
pl => plots(i)
! Get pointer to plot XML node
node_plot = node_plot_list(i)
! Copy data into plots
if (check_for_node(node_plot, "id")) then
call get_node_value(node_plot, "id", pl % id)
else
call fatal_error("Must specify plot id in plots XML file.")
end if
! Check to make sure 'id' hasn't been used
if (plot_dict % has(pl % id)) then
call fatal_error("Two or more plots use the same unique ID: " &
// to_str(pl % id))
end if
! Copy plot type
temp_str = 'slice'
if (check_for_node(node_plot, "type")) &
call get_node_value(node_plot, "type", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("slice")
pl % type = PLOT_TYPE_SLICE
case ("voxel")
pl % type = PLOT_TYPE_VOXEL
case default
call fatal_error("Unsupported plot type '" // trim(temp_str) &
// "' in plot " // trim(to_str(pl % id)))
end select
! Set output file path
filename = "plot_" // trim(to_str(pl % id))
if (check_for_node(node_plot, "filename")) &
call get_node_value(node_plot, "filename", filename)
select case (pl % type)
case (PLOT_TYPE_SLICE)
pl % path_plot = trim(path_input) // trim(filename) // ".ppm"
case (PLOT_TYPE_VOXEL)
pl % path_plot = trim(path_input) // trim(filename) // ".h5"
end select
! Copy plot pixel size
if (pl % type == PLOT_TYPE_SLICE) then
if (node_word_count(node_plot, "pixels") == 2) then
call get_node_array(node_plot, "pixels", pl % pixels(1:2))
else
call fatal_error("<pixels> must be length 2 in slice plot " &
// trim(to_str(pl % id)))
end if
else if (pl % type == PLOT_TYPE_VOXEL) then
if (node_word_count(node_plot, "pixels") == 3) then
call get_node_array(node_plot, "pixels", pl % pixels(1:3))
else
call fatal_error("<pixels> must be length 3 in voxel plot " &
// trim(to_str(pl % id)))
end if
end if
! Copy plot background color
if (check_for_node(node_plot, "background")) then
if (pl % type == PLOT_TYPE_VOXEL) then
if (master) call warning("Background color ignored in voxel plot " &
// trim(to_str(pl % id)))
end if
if (node_word_count(node_plot, "background") == 3) then
call get_node_array(node_plot, "background", pl % not_found % rgb)
else
call fatal_error("Bad background RGB in plot " &
// trim(to_str(pl % id)))
end if
else
pl % not_found % rgb = (/ 255, 255, 255 /)
end if
! Copy plot basis
if (pl % type == PLOT_TYPE_SLICE) then
temp_str = 'xy'
if (check_for_node(node_plot, "basis")) &
call get_node_value(node_plot, "basis", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("xy")
pl % basis = PLOT_BASIS_XY
case ("xz")
pl % basis = PLOT_BASIS_XZ
case ("yz")
pl % basis = PLOT_BASIS_YZ
case default
call fatal_error("Unsupported plot basis '" // trim(temp_str) &
// "' in plot " // trim(to_str(pl % id)))
end select
end if
! Copy plotting origin
if (node_word_count(node_plot, "origin") == 3) then
call get_node_array(node_plot, "origin", pl % origin)
else
call fatal_error("Origin must be length 3 in plot " &
// trim(to_str(pl % id)))
end if
! Copy plotting width
if (pl % type == PLOT_TYPE_SLICE) then
if (node_word_count(node_plot, "width") == 2) then
call get_node_array(node_plot, "width", pl % width(1:2))
else
call fatal_error("<width> must be length 2 in slice plot " &
// trim(to_str(pl % id)))
end if
else if (pl % type == PLOT_TYPE_VOXEL) then
if (node_word_count(node_plot, "width") == 3) then
call get_node_array(node_plot, "width", pl % width(1:3))
else
call fatal_error("<width> must be length 3 in voxel plot " &
// trim(to_str(pl % id)))
end if
end if
! Copy plot cell universe level
if (check_for_node(node_plot, "level")) then
call get_node_value(node_plot, "level", pl % level)
if (pl % level < 0) then
call fatal_error("Bad universe level in plot " &
// trim(to_str(pl % id)))
end if
else
pl % level = PLOT_LEVEL_LOWEST
end if
! Copy plot color type and initialize all colors randomly
temp_str = "cell"
if (check_for_node(node_plot, "color_by")) &
call get_node_value(node_plot, "color_by", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("cell")
pl % color_by = PLOT_COLOR_CELLS
allocate(pl % colors(n_cells))
do j = 1, n_cells
pl % colors(j) % rgb(1) = int(prn()*255)
pl % colors(j) % rgb(2) = int(prn()*255)
pl % colors(j) % rgb(3) = int(prn()*255)
end do
case ("material")
pl % color_by = PLOT_COLOR_MATS
allocate(pl % colors(n_materials))
do j = 1, n_materials
pl % colors(j) % rgb(1) = int(prn()*255)
pl % colors(j) % rgb(2) = int(prn()*255)
pl % colors(j) % rgb(3) = int(prn()*255)
end do
case default
call fatal_error("Unsupported plot color type '" // trim(temp_str) &
// "' in plot " // trim(to_str(pl % id)))
end select
! Get the number of <color> nodes and get a list of them
call get_node_list(node_plot, "color", node_col_list)
n_cols = size(node_col_list)
! Copy user specified colors
if (n_cols /= 0) then
if (pl % type == PLOT_TYPE_VOXEL) then
if (master) call warning("Color specifications ignored in voxel &
&plot " // trim(to_str(pl % id)))
end if
do j = 1, n_cols
! Get pointer to color spec XML node
node_col = node_col_list(j)
! Check and make sure 3 values are specified for RGB
if (node_word_count(node_col, "rgb") /= 3) then
call fatal_error("Bad RGB in plot " &
// trim(to_str(pl % id)))
end if
! Ensure that there is an id for this color specification
if (check_for_node(node_col, "id")) then
call get_node_value(node_col, "id", col_id)
else
call fatal_error("Must specify id for color specification in &
&plot " // trim(to_str(pl % id)))
end if
! Add RGB
if (pl % color_by == PLOT_COLOR_CELLS) then
if (cell_dict % has(col_id)) then
col_id = cell_dict % get(col_id)
call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
else
call fatal_error("Could not find cell " // trim(to_str(col_id)) &
// " specified in plot " // trim(to_str(pl % id)))
end if
else if (pl % color_by == PLOT_COLOR_MATS) then
if (material_dict % has(col_id)) then
col_id = material_dict % get(col_id)
call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
else
call fatal_error("Could not find material " &
// trim(to_str(col_id)) // " specified in plot " &
// trim(to_str(pl % id)))
end if
end if
end do
end if
! Deal with meshlines
call get_node_list(node_plot, "meshlines", node_meshline_list)
n_meshlines = size(node_meshline_list)
if (n_meshlines /= 0) then
if (pl % type == PLOT_TYPE_VOXEL) then
call warning("Meshlines ignored in voxel plot " &
// trim(to_str(pl % id)))
end if
select case(n_meshlines)
case (0)
! Skip if no meshlines are specified
case (1)
! Get pointer to meshlines
node_meshlines = node_meshline_list(1)
! Check mesh type
if (check_for_node(node_meshlines, "meshtype")) then
call get_node_value(node_meshlines, "meshtype", meshtype)
else
call fatal_error("Must specify a meshtype for meshlines &
&specification in plot " // trim(to_str(pl % id)))
end if
! Ensure that there is a linewidth for this meshlines specification
if (check_for_node(node_meshlines, "linewidth")) then
call get_node_value(node_meshlines, "linewidth", &
pl % meshlines_width)
else
call fatal_error("Must specify a linewidth for meshlines &
&specification in plot " // trim(to_str(pl % id)))
end if
! Check for color
if (check_for_node(node_meshlines, "color")) then
! Check and make sure 3 values are specified for RGB
if (node_word_count(node_meshlines, "color") /= 3) then
call fatal_error("Bad RGB for meshlines color in plot " &
// trim(to_str(pl % id)))
end if
call get_node_array(node_meshlines, "color", &
pl % meshlines_color % rgb)
else
pl % meshlines_color % rgb = (/ 0, 0, 0 /)
end if
! Set mesh based on type
select case (trim(meshtype))
case ('ufs')
if (index_ufs_mesh < 0) then
call fatal_error("No UFS mesh for meshlines on plot " &
// trim(to_str(pl % id)))
end if
pl % index_meshlines_mesh = index_ufs_mesh
case ('cmfd')
if (.not. cmfd_run) then
call fatal_error("Need CMFD run to plot CMFD mesh for &
&meshlines on plot " // trim(to_str(pl % id)))
end if
pl % index_meshlines_mesh = index_cmfd_mesh
case ('entropy')
if (index_entropy_mesh < 0) then
call fatal_error("No entropy mesh for meshlines on plot " &
// trim(to_str(pl % id)))
end if
pl % index_meshlines_mesh = index_entropy_mesh
case ('tally')
! Ensure that there is a mesh id if the type is tally
if (check_for_node(node_meshlines, "id")) then
call get_node_value(node_meshlines, "id", meshid)
else
call fatal_error("Must specify a mesh id for meshlines tally &
&mesh specification in plot " // trim(to_str(pl % id)))
end if
! Check if the specified tally mesh exists
err = openmc_get_mesh_index(meshid, idx)
if (err /= 0) then
call fatal_error("Could not find mesh " &
// trim(to_str(meshid)) // " specified in meshlines for &
&plot " // trim(to_str(pl % id)))
end if
pl % index_meshlines_mesh = idx
case default
call fatal_error("Invalid type for meshlines on plot " &
// trim(to_str(pl % id)) // ": " // trim(meshtype))
end select
case default
call fatal_error("Mutliple meshlines specified in plot " &
// trim(to_str(pl % id)))
end select
end if
! Deal with masks
call get_node_list(node_plot, "mask", node_mask_list)
n_masks = size(node_mask_list)
if (n_masks /= 0) then
if (pl % type == PLOT_TYPE_VOXEL) then
if (master) call warning("Mask ignored in voxel plot " &
// trim(to_str(pl % id)))
end if
select case(n_masks)
case default
call fatal_error("Mutliple masks specified in plot " &
// trim(to_str(pl % id)))
case (1)
! Get pointer to mask
node_mask = node_mask_list(1)
! Determine how many components there are and allocate
n_comp = 0
n_comp = node_word_count(node_mask, "components")
if (n_comp == 0) then
call fatal_error("Missing <components> in mask of plot " &
// trim(to_str(pl % id)))
end if
allocate(iarray(n_comp))
call get_node_array(node_mask, "components", iarray)
! First we need to change the user-specified identifiers to indices
! in the cell and material arrays
do j=1, n_comp
col_id = iarray(j)
if (pl % color_by == PLOT_COLOR_CELLS) then
if (cell_dict % has(col_id)) then
iarray(j) = cell_dict % get(col_id)
else
call fatal_error("Could not find cell " &
// trim(to_str(col_id)) // " specified in the mask in &
&plot " // trim(to_str(pl % id)))
end if
else if (pl % color_by == PLOT_COLOR_MATS) then
if (material_dict % has(col_id)) then
iarray(j) = material_dict % get(col_id)
else
call fatal_error("Could not find material " &
// trim(to_str(col_id)) // " specified in the mask in &
&plot " // trim(to_str(pl % id)))
end if
end if
end do
! Alter colors based on mask information
do j = 1, size(pl % colors)
if (.not. any(j == iarray)) then
if (check_for_node(node_mask, "background")) then
call get_node_array(node_mask, "background", pl % colors(j) % rgb)
else
pl % colors(j) % rgb(:) = [255, 255, 255]
end if
end if
end do
deallocate(iarray)
end select
end if
! Add plot to dictionary
call plot_dict % set(pl % id, i)
end do READ_PLOTS
call read_plots(root % ptr)
! Close plots XML file
call doc % clear()

74
src/output.F90
View file

@ -14,7 +14,6 @@ module output
use mgxs_interface
use nuclide_header
use particle_header, only: LocalCoord, Particle
use plot_header
use sab_header, only: SAlphaBeta
use settings
use simulation_header
@ -412,79 +411,6 @@ contains
end subroutine print_batch_keff
!===============================================================================
! PRINT_PLOT displays selected options for plotting
!===============================================================================
subroutine print_plot()
integer :: i ! loop index for plots
type(ObjectPlot), pointer :: pl
! Display header for plotting
call header("PLOTTING SUMMARY", 5)
do i = 1, n_plots
pl => plots(i)
! Plot id
write(ou,100) "Plot ID:", trim(to_str(pl % id))
! Plot filename
write(ou,100) "Plot file:", trim(pl % path_plot)
! Plot level
write(ou,100) "Universe depth:", trim(to_str(pl % level))
! Plot type
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Plot Type:", "Slice"
else if (pl % type == PLOT_TYPE_VOXEL) then
write(ou,100) "Plot Type:", "Voxel"
end if
! Plot parameters
write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2)))
else if (pl % type == PLOT_TYPE_VOXEL) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2))) // &
" " // trim(to_str(pl % width(3)))
end if
if (pl % color_by == PLOT_COLOR_CELLS) then
write(ou,100) "Coloring:", "Cells"
else if (pl % color_by == PLOT_COLOR_MATS) then
write(ou,100) "Coloring:", "Materials"
end if
if (pl % type == PLOT_TYPE_SLICE) then
select case (pl % basis)
case (PLOT_BASIS_XY)
write(ou,100) "Basis:", "xy"
case (PLOT_BASIS_XZ)
write(ou,100) "Basis:", "xz"
case (PLOT_BASIS_YZ)
write(ou,100) "Basis:", "yz"
end select
write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
trim(to_str(pl % pixels(2)))
else if (pl % type == PLOT_TYPE_VOXEL) then
write(ou,100) "Voxels:", trim(to_str(pl % pixels(1))) // " " // &
trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3)))
end if
write(ou,*)
end do
! Format descriptor for columns
100 format (1X,A,T25,A)
end subroutine print_plot
!===============================================================================
! PRINT_RUNTIME displays the total time elapsed for the entire run, for
! initialization, for computation, and for intergeneration synchronization.

92
src/output.cpp
View file

@ -5,13 +5,14 @@
#include <iomanip> // for setw
#include <iostream>
#include <sstream>
#include <ctime>
#include "openmc/cell.h"
#include "openmc/geometry.h"
#include "openmc/message_passing.h"
#include "openmc/capi.h"
#include "openmc/settings.h"
#include "openmc/plot.h"
namespace openmc {
@ -35,13 +36,98 @@ header(const char* msg, int level) {
// Print header based on verbosity level.
if (settings::verbosity >= level) {
std::cout << out.str() << std::endl << std::endl;
std::cout << out.str() << "\n\n";
}
}
std::string time_stamp()
{
int base_year = 1990;
std::stringstream ts;
std::time_t t = std::time(0); // get time now
std::tm* now = std::localtime(&t);
ts << now->tm_year + base_year << "-" << now->tm_mon
<< "-" << now->tm_mday << " " << now->tm_hour
<< ":" << now->tm_min << ":" << now->tm_sec;
return ts.str();
}
//==============================================================================
void print_overlap_check() {
//===============================================================================
// PRINT_PLOT displays selected options for plotting
//===============================================================================
void print_plot() {
header("PLOTTING SUMMARY", 5);
for (auto pl : plots) {
// Plot id
std::cout << "Plot ID: " << pl.id_ << "\n";
// Plot filename
std::cout << "Plot file: " << pl.path_plot_ << "\n";
// Plot level
std::cout << "Universe depth: " << pl.level_ << "\n";
// Plot type
if (PlotType::slice == pl.type_) {
std::cout << "Plot Type: Slice" << "\n";
} else if (PlotType::voxel == pl.type_) {
std::cout << "Plot Type: Voxel" << "\n";
}
// Plot parameters
std::cout << "Origin: " << pl.origin_[0] << " "
<< pl.origin_[1] << " "
<< pl.origin_[2] << "\n";
if (PlotType::slice == pl.type_) {
std::cout << std::setprecision(4)
<< "Width: "
<< pl.width_[0] << " "
<< pl.width_[1] << "\n";
} else if (PlotType::voxel == pl.type_) {
std::cout << std::setprecision(4)
<< "Width: "
<< pl.width_[0] << " "
<< pl.width_[1] << " "
<< pl.width_[2] << "\n";
}
if (PlotColorBy::cells == pl.color_by_) {
std::cout << "Coloring: Cells" << "\n";
} else if (PlotColorBy::mats == pl.color_by_) {
std::cout << "Coloring: Materials" << "\n";
}
if (PlotType::slice == pl.type_) {
switch(pl.basis_) {
case PlotBasis::xy:
std::cout << "Basis: XY" << "\n";
break;
case PlotBasis::xz:
std::cout << "Basis: XZ" << "\n";
break;
case PlotBasis::yz:
std::cout << "Basis: YZ" << "\n";
break;
}
std::cout << "Pixels: " << pl.pixels_[0] << " "
<< pl.pixels_[1] << " " << "\n";
} else if (PlotType::voxel == pl.type_) {
std::cout << "Voxels: " << pl.pixels_[0] << " "
<< pl.pixels_[1] << " "
<< pl.pixels_[2] << "\n";
}
std::cout << "\n";
}
}
void
print_overlap_check() {
#ifdef OPENMC_MPI
std::vector<int64_t> temp(overlap_check_count);
int err = MPI_Reduce(temp.data(), overlap_check_count.data(),

456
src/plot.F90
View file

@ -1,456 +0,0 @@
module plot
use, intrinsic :: ISO_C_BINDING
use constants
use error, only: fatal_error, write_message
use geometry, only: find_cell, check_cell_overlap
use geometry_header, only: Cell, root_universe, cells
use hdf5_interface
use output, only: time_stamp
use material_header, only: materials
use mesh_header, only: meshes, RegularMesh
use particle_header
use plot_header
use progress_header, only: ProgressBar
use settings, only: check_overlaps
use string, only: to_str
implicit none
private
public :: openmc_plot_geometry
integer, parameter :: RED = 1
integer, parameter :: GREEN = 2
integer, parameter :: BLUE = 3
contains
!===============================================================================
! RUN_PLOT controls the logic for making one or many plots
!===============================================================================
function openmc_plot_geometry() result(err) bind(C)
integer(C_INT) :: err
integer :: i ! loop index for plots
do i = 1, n_plots
associate (pl => plots(i))
! Display output message
call write_message("Processing plot " // trim(to_str(pl % id)) &
// ": " // trim(pl % path_plot) // " ...", 5)
if (pl % type == PLOT_TYPE_SLICE) then
! create 2d image
call create_ppm(pl)
else if (pl % type == PLOT_TYPE_VOXEL) then
! create dump for 3D silomesh utility script
call create_voxel(pl)
end if
end associate
end do
err = 0
end function openmc_plot_geometry
!===============================================================================
! POSITION_RGB computes the red/green/blue values for a given plot with the
! current particle's position
!===============================================================================
subroutine position_rgb(p, pl, rgb, id)
type(Particle), intent(inout) :: p
type(ObjectPlot), intent(in) :: pl
integer, intent(out) :: rgb(3)
integer, intent(out) :: id
integer :: j
logical :: found_cell
p % n_coord = 1
call find_cell(p, found_cell)
j = p % n_coord
if (check_overlaps) call check_cell_overlap(p)
! Set coordinate level if specified
if (pl % level >= 0) j = pl % level + 1
if (.not. found_cell) then
! If no cell, revert to default color
rgb = pl % not_found % rgb
id = -1
else
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
associate (c => cells(p % coord(j) % cell + 1))
if (c % type() == FILL_UNIVERSE) then
! If we stopped on a middle universe level, treat as if not found
rgb = pl % not_found % rgb
id = -1
else if (p % material == MATERIAL_VOID) then
! By default, color void cells white
rgb = 255
id = -1
else
rgb = pl % colors(p % material) % rgb
id = materials(p % material) % id()
end if
end associate
else if (pl % color_by == PLOT_COLOR_CELLS) then
! Assign color based on cell
rgb = pl % colors(p % coord(j) % cell + 1) % rgb
id = cells(p % coord(j) % cell + 1) % id()
else
rgb = 0
id = -1
end if
end if
end subroutine position_rgb
!===============================================================================
! CREATE_PPM creates an image based on user input from a plots.xml <plot>
! specification in the portable pixmap format (PPM)
!===============================================================================
subroutine create_ppm(pl)
type(ObjectPlot), intent(in) :: pl
integer :: in_i
integer :: out_i
integer :: x, y ! pixel location
integer :: rgb(3) ! colors (red, green, blue) from 0-255
integer :: id
integer :: height, width
real(8) :: in_pixel
real(8) :: out_pixel
real(8) :: xyz(3)
integer, allocatable :: data(:,:,:)
type(Particle) :: p
width = pl % pixels(1)
height = pl % pixels(2)
in_pixel = pl % width(1)/dble(width)
out_pixel = pl % width(2)/dble(height)
! Allocate and initialize results array
allocate(data(3, width, height))
data(:,:,:) = 0
if (pl % basis == PLOT_BASIS_XY) then
in_i = 1
out_i = 2
xyz(1) = pl % origin(1) - pl % width(1) / TWO
xyz(2) = pl % origin(2) + pl % width(2) / TWO
xyz(3) = pl % origin(3)
else if (pl % basis == PLOT_BASIS_XZ) then
in_i = 1
out_i = 3
xyz(1) = pl % origin(1) - pl % width(1) / TWO
xyz(2) = pl % origin(2)
xyz(3) = pl % origin(3) + pl % width(2) / TWO
else if (pl % basis == PLOT_BASIS_YZ) then
in_i = 2
out_i = 3
xyz(1) = pl % origin(1)
xyz(2) = pl % origin(2) - pl % width(1) / TWO
xyz(3) = pl % origin(3) + pl % width(2) / TWO
end if
! allocate and initialize particle
call particle_initialize(p)
p % coord(1) % xyz = xyz
p % coord(1) % uvw = [ HALF, HALF, HALF ]
p % coord(1) % universe = root_universe
!$omp parallel do firstprivate(p) private(x, rgb, id) reduction(+ : data)
do y = 1, height
! Set y coordinate
p % coord(1) % xyz(out_i) = xyz(out_i) - out_pixel*(y - 1)
do x = 1, width
! Set x coordinate
p % coord(1) % xyz(in_i) = xyz(in_i) + in_pixel*(x - 1)
! get pixel color
call position_rgb(p, pl, rgb, id)
! Create a pixel at (x,y) with color (r,g,b)
data(:, x, y) = rgb
end do
end do
!$omp end parallel do
! Draw tally mesh boundaries on the image if requested
if (pl % index_meshlines_mesh >= 0) call draw_mesh_lines(pl, data)
! Write out the ppm to a file
call output_ppm(pl, data)
end subroutine create_ppm
!===============================================================================
! DRAW_MESH_LINES draws mesh line boundaries on an image
!===============================================================================
subroutine draw_mesh_lines(pl, data)
type(ObjectPlot), intent(in) :: pl
integer, intent(inout) :: data(:,:,:)
logical :: in_mesh
integer :: out_, in_ ! pixel location
integer :: rgb(3) ! RGB color for meshlines pixels
integer :: outrange(2), inrange(2) ! range of pixel locations
integer :: i, j ! loop indices
integer :: plus
integer :: ijk_ll(3) ! mesh bin ijk indicies of plot lower left
integer :: ijk_ur(3) ! mesh bin ijk indicies of plot upper right
integer :: outer, inner
real(8) :: frac
real(8) :: width(3) ! real widths of the plot
real(8) :: xyz_ll_plot(3) ! lower left xyz of plot image
real(8) :: xyz_ur_plot(3) ! upper right xyz of plot image
real(8) :: xyz_ll(3) ! lower left xyz
real(8) :: xyz_ur(3) ! upper right xyz
type(RegularMesh) :: m
rgb(:) = pl % meshlines_color % rgb
select case (pl % basis)
case(PLOT_BASIS_XY)
outer = 1
inner = 2
case(PLOT_BASIS_XZ)
outer = 1
inner = 3
case(PLOT_BASIS_YZ)
outer = 2
inner = 3
end select
xyz_ll_plot = pl % origin
xyz_ur_plot = pl % origin
xyz_ll_plot(outer) = pl % origin(outer) - pl % width(1) / TWO
xyz_ll_plot(inner) = pl % origin(inner) - pl % width(2) / TWO
xyz_ur_plot(outer) = pl % origin(outer) + pl % width(1) / TWO
xyz_ur_plot(inner) = pl % origin(inner) + pl % width(2) / TWO
width = xyz_ur_plot - xyz_ll_plot
m = meshes(pl % index_meshlines_mesh)
call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension()), in_mesh)
call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension()), in_mesh)
! sweep through all meshbins on this plane and draw borders
do i = ijk_ll(outer), ijk_ur(outer)
do j = ijk_ll(inner), ijk_ur(inner)
! check if we're in the mesh for this ijk
if (i > 0 .and. i <= m % dimension(outer) .and. &
j > 0 .and. j <= m % dimension(inner)) then
! get xyz's of lower left and upper right of this mesh cell
xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)
xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1)
xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i
xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j
! map the xyz ranges to pixel ranges
frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
outrange(1) = int(frac * real(pl % pixels(1), 8))
frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
outrange(2) = int(frac * real(pl % pixels(1), 8))
frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8))
frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8))
! draw lines
do out_ = outrange(1), outrange(2)
do plus = 0, pl % meshlines_width
data(:, out_ + 1, inrange(1) + plus + 1) = rgb
data(:, out_ + 1, inrange(2) + plus + 1) = rgb
data(:, out_ + 1, inrange(1) - plus + 1) = rgb
data(:, out_ + 1, inrange(2) - plus + 1) = rgb
end do
end do
do in_ = inrange(1), inrange(2)
do plus = 0, pl % meshlines_width
data(:, outrange(1) + plus + 1, in_ + 1) = rgb
data(:, outrange(2) + plus + 1, in_ + 1) = rgb
data(:, outrange(1) - plus + 1, in_ + 1) = rgb
data(:, outrange(2) - plus + 1, in_ + 1) = rgb
end do
end do
end if
end do
end do
end subroutine draw_mesh_lines
!===============================================================================
! OUTPUT_PPM writes out a previously generated image to a PPM file
!===============================================================================
subroutine output_ppm(pl, data)
type(ObjectPlot), intent(in) :: pl
integer, intent(in) :: data(:,:,:)
integer :: y ! loop index for height
integer :: x ! loop index for width
integer :: unit_plot
! Open PPM file for writing
open(NEWUNIT=unit_plot, FILE=pl % path_plot)
! Write header
write(unit_plot, '(A2)') 'P6'
write(unit_plot, '(I0,'' '',I0)') pl % pixels(1), pl % pixels(2)
write(unit_plot, '(A)') '255'
! Write color for each pixel
do y = 1, pl % pixels(2)
do x = 1, pl % pixels(1)
write(unit_plot, '(3A1)', advance='no') achar(data(RED, x, y)), &
achar(data(GREEN, x, y)), achar(data(BLUE, x, y))
end do
end do
! Close plot file
close(UNIT=unit_plot)
end subroutine output_ppm
!===============================================================================
! 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 cell or material
! id. For 1 million voxels this produces a file of approximately 15MB.
!===============================================================================
subroutine create_voxel(pl)
type(ObjectPlot), intent(in) :: pl
integer(C_INT) :: x, y, z ! voxel location indices
integer :: rgb(3) ! colors (red, green, blue) from 0-255
integer :: id ! id of cell or material
integer(C_INT), target :: data(pl%pixels(3),pl%pixels(2))
integer(HID_T) :: file_id
integer(HID_T) :: dspace
integer(HID_T) :: memspace
integer(HID_T) :: dset
integer(HSIZE_T) :: dims(3)
real(8) :: vox(3) ! x, y, and z voxel widths
real(8) :: ll(3) ! lower left starting point for each sweep direction
type(Particle) :: p
type(ProgressBar) :: progress
interface
subroutine voxel_init(file_id, dims, dspace, dset, memspace) bind(C)
import HID_T, HSIZE_T
integer(HID_T), value :: file_id
integer(HSIZE_T), intent(in) :: dims(*)
integer(HID_T), intent(out) :: dspace
integer(HID_T), intent(out) :: dset
integer(HID_T), intent(out) :: memspace
end subroutine voxel_init
subroutine voxel_write_slice(x, dspace, dset, memspace, buf) bind(C)
import C_INT, HID_T, C_PTR
integer(C_INT), value :: x
integer(HID_T), value :: dspace
integer(HID_T), value :: dset
integer(HID_T), value :: memspace
type(C_PTR), value :: buf
end subroutine voxel_write_slice
subroutine voxel_finalize(dspace, dset, memspace) bind(C)
import HID_T
integer(HID_T), value :: dspace
integer(HID_T), value :: dset
integer(HID_T), value :: memspace
end subroutine voxel_finalize
end interface
! compute voxel widths in each direction
vox = pl % width/dble(pl % pixels)
! initial particle position
ll = pl % origin - pl % width / TWO
! allocate and initialize particle
call particle_initialize(p)
p % coord(1) % xyz = ll
p % coord(1) % uvw = [ HALF, HALF, HALF ]
p % coord(1) % universe = root_universe
! Open binary plot file for writing
file_id = file_open(pl%path_plot, 'w')
! write header info
call write_attribute(file_id, "filetype", 'voxel')
call write_attribute(file_id, "version", VERSION_VOXEL)
call write_attribute(file_id, "openmc_version", VERSION)
#ifdef GIT_SHA1
call write_attribute(file_id, "git_sha1", GIT_SHA1)
#endif
! Write current date and time
call write_attribute(file_id, "date_and_time", time_stamp())
call write_attribute(file_id, "num_voxels", pl%pixels)
call write_attribute(file_id, "voxel_width", vox)
call write_attribute(file_id, "lower_left", ll)
! Create dataset for voxel data -- note that the dimensions are reversed
! since we want the order in the file to be z, y, x
dims(:) = pl % pixels
call voxel_init(file_id, dims, dspace, dset, memspace)
! move to center of voxels
ll = ll + vox / TWO
do x = 1, pl % pixels(1)
call progress % set_value(dble(x)/dble(pl % pixels(1))*100)
do y = 1, pl % pixels(2)
do z = 1, pl % pixels(3)
! get voxel color
call position_rgb(p, pl, rgb, id)
! write to plot file
data(z,y) = id
! advance particle in z direction
p % coord(1) % xyz(3) = p % coord(1) % xyz(3) + vox(3)
end do
! advance particle in y direction
p % coord(1) % xyz(2) = p % coord(1) % xyz(2) + vox(2)
p % coord(1) % xyz(3) = ll(3)
end do
! advance particle in y direction
p % coord(1) % xyz(1) = p % coord(1) % xyz(1) + vox(1)
p % coord(1) % xyz(2) = ll(2)
p % coord(1) % xyz(3) = ll(3)
! Write to HDF5 dataset
call voxel_write_slice(x, dspace, dset, memspace, c_loc(data))
end do
call voxel_finalize(dspace, dset, memspace)
call file_close(file_id)
end subroutine create_voxel
end module plot

931
src/plot.cpp
View file

@ -1,10 +1,931 @@
#include <fstream>
#include <sstream>
#include "openmc/plot.h"
#include "openmc/constants.h"
#include "openmc/settings.h"
#include "openmc/error.h"
#include "openmc/particle.h"
#include "openmc/geometry.h"
#include "openmc/cell.h"
#include "openmc/material.h"
#include "openmc/string_functions.h"
#include "openmc/mesh.h"
#include "openmc/output.h"
#include "openmc/hdf5_interface.h"
#include "openmc/random_lcg.h"
#include "openmc/output.h"
#include "openmc/progress_bar.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
int PLOT_LEVEL_LOWEST = -1;
std::unordered_map<int, int> plot_map;
int n_plots;
std::vector<Plot> plots;
const RGBColor WHITE = {255, 255, 255};
//==============================================================================
// RUN_PLOT controls the logic for making one or many plots
//==============================================================================
extern "C"
int openmc_plot_geometry()
{
int err;
for (auto pl : plots) {
std::stringstream ss;
ss << "Processing plot " << pl.id_ << ": "
<< pl.path_plot_ << "...";
write_message(ss.str(), 5);
if (PlotType::slice == pl.type_) {
// create 2D image
create_ppm(pl);
} else if (PlotType::voxel == pl.type_) {
// create voxel file for 3D viewing
create_voxel(pl);
}
}
return 0;
}
void
voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
hid_t* memspace)
read_plots(pugi::xml_node* plots_node)
{
n_plots = 0;
for (auto node : plots_node->children("plot")) {
Plot pl(node);
plots.push_back(pl);
plot_map[pl.id_] = n_plots++;
}
}
//==============================================================================
// CREATE_PPM creates an image based on user input from a plots.xml <plot>
// specification in the portable pixmap format (PPM)
//==============================================================================
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;
data.resize({width, height});
int in_i, out_i;
double xyz[3];
switch(pl.basis_) {
case PlotBasis::xy :
in_i = 0;
out_i = 1;
xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
xyz[1] = pl.origin_[1] + pl.width_[1] / 2.;
xyz[2] = pl.origin_[2];
break;
case PlotBasis::xz :
in_i = 0;
out_i = 2;
xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
xyz[1] = pl.origin_[1];
xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
break;
case PlotBasis::yz :
in_i = 1;
out_i = 2;
xyz[0] = pl.origin_[0];
xyz[1] = pl.origin_[1] - pl.width_[0] / 2.;
xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
break;
}
double dir[3] = {0.5, 0.5, 0.5};
#pragma omp parallel
{
Particle p;
p.initialize();
std::copy(xyz, xyz+3, p.coord[0].xyz);
std::copy(dir, dir+3, p.coord[0].uvw);
p.coord[0].universe = openmc_root_universe;
#pragma omp for
for (int y = 0; y < height; y++) {
p.coord[0].xyz[out_i] = xyz[out_i] - out_pixel * y;
for (int x = 0; x < width; x++) {
// local variables
RGBColor rgb;
int id;
p.coord[0].xyz[in_i] = xyz[in_i] + in_pixel * x;
position_rgb(p, pl, rgb, id);
data(x,y) = rgb;
}
}
}
// draw mesh lines if present
if (pl.index_meshlines_mesh_ >= 0) {draw_mesh_lines(pl, data);}
// write ppm data to file
output_ppm(pl, data);
}
void
Plot::set_id(pugi::xml_node plot_node)
{
// Copy data into plots
if (check_for_node(plot_node, "id")) {
id_ = std::stoi(get_node_value(plot_node, "id"));
} else {
fatal_error("Must specify plot id in plots XML file.");
}
// Check to make sure 'id' hasn't been used
if (plot_map.find(id_) != plot_map.end()) {
std::stringstream err_msg;
err_msg << "Two or more plots use the same unique ID: " << id_;
fatal_error(err_msg.str());
}
}
void
Plot::set_type(pugi::xml_node plot_node)
{
// Copy plot type
// Default is slice
type_ = PlotType::slice;
// check type specified on plot node
if (check_for_node(plot_node, "type")) {
std::string type_str = get_node_value(plot_node, "type", true);
// set type using node value
if (type_str == "slice") {
type_ = PlotType::slice;
}
else if (type_str == "voxel") {
type_ = PlotType::voxel;
} else {
// if we're here, something is wrong
std::stringstream err_msg;
err_msg << "Unsupported plot type '" << type_str
<< "' in plot " << id_;
fatal_error(err_msg.str());
}
}
}
void
Plot::set_output_path(pugi::xml_node plot_node)
{
// Set output file path
std::stringstream filename;
if (check_for_node(plot_node, "filename")) {
filename << get_node_value(plot_node, "filename");
} else {
filename << "plot_" << id_;
}
// add appropriate file extension to name
switch(type_) {
case PlotType::slice:
filename << ".ppm";
break;
case PlotType::voxel:
filename << ".h5";
break;
}
path_plot_ = filename.str();
// Copy plot pixel size
std::vector<int> pxls = get_node_array<int>(plot_node, "pixels");
if (PlotType::slice == type_) {
if (pxls.size() == 2) {
pixels_[0] = pxls[0];
pixels_[1] = pxls[1];
} else {
std::stringstream err_msg;
err_msg << "<pixels> must be length 2 in slice plot "
<< id_;
fatal_error(err_msg.str());
}
} else if (PlotType::voxel == type_) {
if (pxls.size() == 3) {
pixels_[0] = pxls[0];
pixels_[1] = pxls[1];
pixels_[2] = pxls[2];
} else {
std::stringstream err_msg;
err_msg << "<pixels> must be length 3 in voxel plot "
<< id_;
fatal_error(err_msg.str());
}
}
}
void
Plot::set_bg_color(pugi::xml_node plot_node)
{
// Copy plot background color
if (check_for_node(plot_node, "background")) {
std::vector<int> bg_rgb = get_node_array<int>(plot_node, "background");
if (PlotType::voxel == type_) {
if (openmc_master) {
std::stringstream err_msg;
err_msg << "Background color ignored in voxel plot "
<< id_;
warning(err_msg.str());
}
}
if (bg_rgb.size() == 3) {
not_found_ = bg_rgb;
} else {
std::stringstream err_msg;
err_msg << "Bad background RGB in plot "
<< id_;
fatal_error(err_msg);
}
} else {
// default to a white background
not_found_ = WHITE;
}
}
void
Plot::set_basis(pugi::xml_node plot_node)
{
// Copy plot basis
if (PlotType::slice == type_) {
std::string pl_basis = "xy";
if (check_for_node(plot_node, "basis")) {
pl_basis = get_node_value(plot_node, "basis", true);
}
if ("xy" == pl_basis) {
basis_ = PlotBasis::xy;
} else if ("xz" == pl_basis) {
basis_ = PlotBasis::xz;
} else if ("yz" == pl_basis) {
basis_ = PlotBasis::yz;
} else {
std::stringstream err_msg;
err_msg << "Unsupported plot basis '" << pl_basis
<< "' in plot " << id_;
fatal_error(err_msg);
}
}
}
void
Plot::set_origin(pugi::xml_node plot_node)
{
// Copy plotting origin
std::vector<double> pl_origin = get_node_array<double>(plot_node, "origin");
if (pl_origin.size() == 3) {
origin_[0] = pl_origin[0];
origin_[1] = pl_origin[1];
origin_[2] = pl_origin[2];
} else {
std::stringstream err_msg;
err_msg << "Origin must be length 3 in plot "
<< id_;
fatal_error(err_msg);
}
}
void
Plot::set_width(pugi::xml_node plot_node)
{
// Copy plotting width
std::vector<double> pl_width = get_node_array<double>(plot_node, "width");
if (PlotType::slice == type_) {
if (pl_width.size() == 2) {
width_[0] = pl_width[0];
width_[1] = pl_width[1];
} else {
std::stringstream err_msg;
err_msg << "<width> must be length 2 in slice plot "
<< id_;
fatal_error(err_msg);
}
} else if (PlotType::voxel == type_) {
if (pl_width.size() == 3) {
pl_width = get_node_array<double>(plot_node, "width");
width_[0] = pl_width[0];
width_[1] = pl_width[1];
width_[2] = pl_width[2];
} else {
std::stringstream err_msg;
err_msg << "<width> must be length 3 in voxel plot "
<< id_;
fatal_error(err_msg);
}
}
}
void
Plot::set_universe(pugi::xml_node plot_node)
{
// Copy plot universe level
if (check_for_node(plot_node, "level")) {
level_ = std::stoi(get_node_value(plot_node, "level"));
if (level_ < 0) {
std::stringstream err_msg;
err_msg << "Bad universe level in plot " << id_;
fatal_error(err_msg);
}
} else {
level_ = PLOT_LEVEL_LOWEST;
}
}
void
Plot::set_default_colors(pugi::xml_node plot_node)
{
// Copy plot color type and initialize all colors randomly
std::string pl_color_by = "cell";
if (check_for_node(plot_node, "color_by")) {
pl_color_by = get_node_value(plot_node, "color_by", true);
}
if ("cell" == pl_color_by) {
color_by_ = PlotColorBy::cells;
colors_.resize(n_cells);
for (int i = 0; i < n_cells; i++) {
colors_[i] = random_color();
}
} else if("material" == pl_color_by) {
color_by_ = PlotColorBy::mats;
colors_.resize(n_materials);
for (int i = 0; i < materials.size(); i++) {
colors_[i] = random_color();
}
} else {
std::stringstream err_msg;
err_msg << "Unsupported plot color type '" << pl_color_by
<< "' in plot " << id_;
fatal_error(err_msg);
}
}
void
Plot::set_user_colors(pugi::xml_node plot_node)
{
if (!plot_node.select_nodes("color").empty() && PlotType::voxel == type_) {
if (openmc_master) {
std::stringstream err_msg;
err_msg << "Color specifications ignored in voxel plot "
<< id_;
warning(err_msg);
}
}
for (auto cn : plot_node.children("color")) {
// Make sure 3 values are specified for RGB
std::vector<int> user_rgb = get_node_array<int>(cn, "rgb");
if (user_rgb.size() != 3) {
std::stringstream err_msg;
err_msg << "Bad RGB in plot " << id_;
fatal_error(err_msg);
}
// Ensure that there is an id for this color specification
int col_id;
if (check_for_node(cn, "id")) {
col_id = std::stoi(get_node_value(cn, "id"));
} else {
std::stringstream err_msg;
err_msg << "Must specify id for color specification in plot "
<< id_;
fatal_error(err_msg);
}
// Add RGB
if (PlotColorBy::cells == color_by_) {
if (cell_map.find(col_id) != cell_map.end()) {
col_id = cell_map[col_id];
colors_[col_id] = user_rgb;
} else {
std::stringstream err_msg;
err_msg << "Could not find cell " << col_id
<< " specified in plot " << id_;
fatal_error(err_msg);
}
} else if (PlotColorBy::mats == color_by_) {
if (material_map.find(col_id) != material_map.end()) {
col_id = material_map[col_id];
colors_[col_id] = user_rgb;
} else {
std::stringstream err_msg;
err_msg << "Could not find material " << col_id
<< " specified in plot " << id_;
fatal_error(err_msg);
}
}
} // color node loop
}
void
Plot::set_meshlines(pugi::xml_node plot_node)
{
// Deal with meshlines
pugi::xpath_node_set mesh_line_nodes = plot_node.select_nodes("meshlines");
if (!mesh_line_nodes.empty()) {
if (PlotType::voxel == type_) {
std::stringstream msg;
msg << "Meshlines ignored in voxel plot " << id_;
warning(msg);
}
if (mesh_line_nodes.size() == 1) {
// Get first meshline node
pugi::xml_node meshlines_node = mesh_line_nodes[0].node();
// Check mesh type
std::string meshtype;
if (check_for_node(meshlines_node, "meshtype")) {
meshtype = get_node_value(meshlines_node, "meshtype");
} else {
std::stringstream err_msg;
err_msg << "Must specify a meshtype for meshlines specification in plot " << id_;
fatal_error(err_msg);
}
// Ensure that there is a linewidth for this meshlines specification
std::string meshline_width;
if (check_for_node(meshlines_node, "linewidth")) {
meshline_width = get_node_value(meshlines_node, "linewidth");
meshlines_width_ = std::stoi(meshline_width);
} else {
std::stringstream err_msg;
err_msg << "Must specify a linewidth for meshlines specification in plot " << id_;
fatal_error(err_msg);
}
// Check for color
if (check_for_node(meshlines_node, "color")) {
// Check and make sure 3 values are specified for RGB
std::vector<int> ml_rgb = get_node_array<int>(meshlines_node, "color");
if (ml_rgb.size() != 3) {
std::stringstream err_msg;
err_msg << "Bad RGB for meshlines color in plot " << id_;
fatal_error(err_msg);
}
meshlines_color_ = ml_rgb;
}
// Set mesh based on type
if ("ufs" == meshtype) {
if (settings::index_ufs_mesh < 0) {
std::stringstream err_msg;
err_msg << "No UFS mesh for meshlines on plot " << id_;
fatal_error(err_msg);
} else {
index_meshlines_mesh_ = settings::index_ufs_mesh;
}
} else if ("cmfd" == meshtype) {
if (!settings::cmfd_run) {
std::stringstream err_msg;
err_msg << "Need CMFD run to plot CMFD mesh for meshlines on plot " << id_;
fatal_error(err_msg);
} else {
index_meshlines_mesh_ = settings::index_cmfd_mesh;
}
} else if ("entropy" == meshtype) {
if (settings::index_entropy_mesh < 0) {
std::stringstream err_msg;
err_msg <<"No entropy mesh for meshlines on plot " << id_;
fatal_error(err_msg);
} else {
index_meshlines_mesh_ = settings::index_entropy_mesh;
}
} else if ("tally" == meshtype) {
// Ensure that there is a mesh id if the type is tally
int tally_mesh_id;
if (check_for_node(meshlines_node, "id")) {
tally_mesh_id = std::stoi(get_node_value(meshlines_node, "id"));
} else {
std::stringstream err_msg;
err_msg << "Must specify a mesh id for meshlines tally "
<< "mesh specification in plot " << id_;
fatal_error(err_msg);
}
// find the tally index
int idx;
int err = openmc_get_mesh_index(tally_mesh_id, &idx);
if (err != 0) {
std::stringstream err_msg;
err_msg << "Could not find mesh " << tally_mesh_id
<< " specified in meshlines for plot " << id_;
fatal_error(err_msg);
}
index_meshlines_mesh_ = idx;
} else {
std::stringstream err_msg;
err_msg << "Invalid type for meshlines on plot " << id_ ;
fatal_error(err_msg);
}
} else {
std::stringstream err_msg;
err_msg << "Mutliple meshlines specified in plot " << id_;
fatal_error(err_msg);
}
}
}
void
Plot::set_mask(pugi::xml_node plot_node)
{
// Deal with masks
pugi::xpath_node_set mask_nodes = plot_node.select_nodes("mask");
if (!mask_nodes.empty()) {
if (PlotType::voxel == type_) {
if (openmc_master) {
std::stringstream wrn_msg;
wrn_msg << "Mask ignored in voxel plot " << id_;
warning(wrn_msg);
}
}
if (mask_nodes.size() == 1) {
// Get pointer to mask
pugi::xml_node mask_node = mask_nodes[0].node();
// Determine how many components there are and allocate
std::vector<int> iarray = get_node_array<int>(mask_node, "components");
if (iarray.size() == 0) {
std::stringstream err_msg;
err_msg << "Missing <components> in mask of plot " << id_;
fatal_error(err_msg);
}
// First we need to change the user-specified identifiers to indices
// in the cell and material arrays
for (auto& col_id : iarray) {
if (PlotColorBy::cells == color_by_) {
if (cell_map.find(col_id) != cell_map.end()) {
col_id = cell_map[col_id];
}
else {
std::stringstream err_msg;
err_msg << "Could not find cell " << col_id
<< " specified in the mask in plot " << id_;
fatal_error(err_msg);
}
} else if (PlotColorBy::mats == color_by_) {
if (material_map.find(col_id) != material_map.end()) {
col_id = material_map[col_id];
}
else {
std::stringstream err_msg;
err_msg << "Could not find material " << col_id
<< " specified in the mask in plot " << id_;
fatal_error(err_msg);
}
}
}
// Alter colors based on mask information
for (int j = 0; j < colors_.size(); j++) {
if (std::find(iarray.begin(), iarray.end(), j) == iarray.end()) {
if (check_for_node(mask_node, "background")) {
std::vector<int> bg_rgb = get_node_array<int>(mask_node, "background");
colors_[j] = bg_rgb;
} else {
colors_[j] = WHITE;
}
}
}
} else {
std::stringstream err_msg;
err_msg << "Mutliple masks specified in plot " << id_;
fatal_error(err_msg);
}
}
}
Plot::Plot(pugi::xml_node plot_node):
index_meshlines_mesh_(-1)
{
set_id(plot_node);
set_type(plot_node);
set_output_path(plot_node);
set_bg_color(plot_node);
set_basis(plot_node);
set_origin(plot_node);
set_width(plot_node);
set_universe(plot_node);
set_default_colors(plot_node);
set_user_colors(plot_node);
set_meshlines(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 = -1;
} else {
if (PlotColorBy::mats == pl.color_by_) {
// Assign color based on material
Cell* c = 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 = -1;
} else if (p.material == MATERIAL_VOID) {
// By default, color void cells white
rgb = WHITE;
id = -1;
} else {
rgb = pl.colors_[p.material - 1];
id = materials[p.material - 1]->id_;
}
} else if (PlotColorBy::cells == pl.color_by_) {
// Assign color based on cell
rgb = pl.colors_[p.coord[j].cell];
id = cells[p.coord[j].cell]->id_;
}
} // endif found_cell
}
//==============================================================================
// OUTPUT_PPM writes out a previously generated image to a PPM file
//==============================================================================
void output_ppm(Plot pl, const ImageData& data)
{
// Open PPM file for writing
std::string fname = pl.path_plot_;
fname = strtrim(fname);
std::ofstream of;
of.open(fname);
// Write header
of << "P6" << "\n";
of << pl.pixels_[0] << " " << pl.pixels_[1] << "\n";
of << "255" << "\n";
of.close();
of.open(fname, std::ios::binary | std::ios::app);
// Write color for each pixel
for (int y = 0; y < pl.pixels_[1]; y++) {
for (int x = 0; x < pl.pixels_[0]; x++) {
RGBColor rgb = data(x,y);
of << rgb.red << rgb.green << rgb.blue;
}
}
// Close file
// THIS IS HERE TO MATCH FORTRAN VERSION, NOT TECHNICALLY NECESSARY
of << "\n";
of.close();
}
//==============================================================================
// DRAW_MESH_LINES draws mesh line boundaries on an image
//==============================================================================
void draw_mesh_lines(Plot pl, ImageData& data)
{
RGBColor rgb;
rgb = pl.meshlines_color_;
int outer, inner;
switch(pl.basis_) {
case PlotBasis::xy :
outer = 0;
inner = 1;
break;
case PlotBasis::xz :
outer = 0;
inner = 2;
break;
case PlotBasis::yz :
outer = 1;
inner = 2;
break;
}
double xyz_ll_plot[3], xyz_ur_plot[3];
xyz_ll_plot[0] = pl.origin_[0];
xyz_ll_plot[1] = pl.origin_[1];
xyz_ll_plot[2] = pl.origin_[2];
xyz_ur_plot[0] = pl.origin_[0];
xyz_ur_plot[1] = pl.origin_[1];
xyz_ur_plot[2] = pl.origin_[2];
xyz_ll_plot[outer] = pl.origin_[outer] - pl.width_[0] / 2.;
xyz_ll_plot[inner] = pl.origin_[inner] - pl.width_[1] / 2.;
xyz_ur_plot[outer] = pl.origin_[outer] + pl.width_[0] / 2.;
xyz_ur_plot[inner] = pl.origin_[inner] + pl.width_[1] / 2.;
int width[3];
width[0] = xyz_ur_plot[0] - xyz_ll_plot[0];
width[1] = xyz_ur_plot[1] - xyz_ll_plot[1];
width[2] = xyz_ur_plot[2] - xyz_ll_plot[2];
auto& m = meshes[pl.index_meshlines_mesh_];
int ijk_ll[3], ijk_ur[3];
bool in_mesh;
m->get_indices(Position(xyz_ll_plot), &(ijk_ll[0]), &in_mesh);
m->get_indices(Position(xyz_ur_plot), &(ijk_ur[0]), &in_mesh);
// Fortran/C++ index correction
ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++;
double xyz_ll[3], xyz_ur[3];
// sweep through all meshbins on this plane and draw borders
for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) {
for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) {
// check if we're in the mesh for this ijk
if (i > 0 && i <= m->shape_[outer] && j >0 && j <= m->shape_[inner] ) {
int outrange[3], inrange[3];
// get xyz's of lower left and upper right of this mesh cell
xyz_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1);
xyz_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1);
xyz_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i;
xyz_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j;
// map the xyz ranges to pixel ranges
double frac = (xyz_ll[outer] - xyz_ll_plot[outer]) / width[outer];
outrange[0] = int(frac * double(pl.pixels_[0]));
frac = (xyz_ur[outer] - xyz_ll_plot[outer]) / width[outer];
outrange[1] = int(frac * double(pl.pixels_[0]));
frac = (xyz_ur[inner] - xyz_ll_plot[inner]) / width[inner];
inrange[0] = int((1. - frac) * (double)pl.pixels_[1]);
frac = (xyz_ll[inner] - xyz_ll_plot[inner]) / width[inner];
inrange[1] = int((1. - frac) * (double)pl.pixels_[1]);
// draw lines
for (int out_ = outrange[0]; out_ <= outrange[1]; out_++) {
for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
data(out_, inrange[0] + plus) = rgb;
data(out_, inrange[1] + plus) = rgb;
data(out_, inrange[0] - plus) = rgb;
data(out_, inrange[1] - plus) = rgb;
}
}
for (int in_ = inrange[0]; in_ <= inrange[1]; in_++) {
for (int plus = 0; plus <= pl.meshlines_width_; plus++) {
data(outrange[0] + plus, in_) = rgb;
data(outrange[1] + plus, in_) = rgb;
data(outrange[0] - plus, in_) = rgb;
data(outrange[1] - plus, in_) = rgb;
}
}
} // end if(in mesh)
}
} // end outer loops
} // end draw_mesh_lines
//==============================================================================
// 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
// 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];
vox[1] = pl.width_[1]/(double)pl.pixels_[1];
vox[2] = pl.width_[2]/(double)pl.pixels_[2];
// initial particle position
std::array<double, 3> ll;
ll[0] = pl.origin_[0] - pl.width_[0] / 2.;
ll[1] = pl.origin_[1] - pl.width_[1] / 2.;
ll[2] = pl.origin_[2] - pl.width_[2] / 2.;
// allocate and initialize particle
double dir[3] = {0.5, 0.5, 0.5};
Particle p;
p.initialize();
std::copy(ll.begin(), ll.begin()+ll.size(), p.coord[0].xyz);
std::copy(dir, dir+3, p.coord[0].uvw);
p.coord[0].universe = openmc_root_universe;
// Open binary plot file for writing
std::ofstream of;
std::string fname = std::string(pl.path_plot_);
fname = strtrim(fname);
hid_t file_id = file_open(fname, 'w');
// write header info
write_attribute(file_id, "filetype", "voxel");
write_attribute(file_id, "version", VERSION_VOXEL);
write_attribute(file_id, "openmc_version", VERSION);
#ifdef GIT_SHA1
write_attribute(file_id, "git_sha1", GIT_SHA1);
#endif
// 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_);
write_attribute(file_id, "voxel_width", vox);
write_attribute(file_id, "lower_left", ll);
// Create dataset for voxel data -- note that the dimensions are reversed
// since we want the order in the file to be z, y, x
hsize_t dims[3];
dims[0] = pl.pixels_[0];
dims[1] = pl.pixels_[1];
dims[2] = pl.pixels_[2];
hid_t dspace, dset, memspace;
voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
// move to center of voxels
ll[0] = ll[0] + vox[0] / 2.;
ll[1] = ll[1] + vox[1] / 2.;
ll[2] = ll[2] + vox[2] / 2.;
int data[pl.pixels_[1]][pl.pixels_[2]];
ProgressBar pb;
RGBColor rgb;
int id;
for (int x = 0; x < pl.pixels_[0]; x++) {
pb.set_value(100.*(double)x/(double)(pl.pixels_[0]-1));
for (int y = 0; y < pl.pixels_[1]; y++) {
for (int z = 0; z < pl.pixels_[2]; z++) {
// get voxel color
position_rgb(p, pl, rgb, id);
// write to plot data
data[y][z] = id;
// advance particle in z direction
p.coord[0].xyz[2] = p.coord[0].xyz[2] + vox[2];
}
// advance particle in y direction
p.coord[0].xyz[1] = p.coord[0].xyz[1] + vox[1];
p.coord[0].xyz[2] = ll[2];
}
// advance particle in x direction
p.coord[0].xyz[0] = p.coord[0].xyz[0] + vox[0];
p.coord[0].xyz[1] = ll[1];
p.coord[0].xyz[2] = ll[2];
// Write to HDF5 dataset
voxel_write_slice(x, dspace, dset, memspace, &(data[0]));
}
voxel_finalize(dspace, dset, memspace);
file_close(file_id);
}
void
voxel_init(hid_t file_id, const hsize_t* dims,
hid_t* dspace, hid_t* dset, hid_t* memspace)
{
// Create dataspace/dataset for voxel data
*dspace = H5Screate_simple(3, dims, nullptr);
@ -25,7 +946,7 @@ voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
void
voxel_write_slice(int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf)
{
hssize_t offset[3] {x - 1, 0, 0};
hssize_t offset[3] {x, 0, 0};
H5Soffset_simple(dspace, offset);
H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf);
}
@ -39,4 +960,8 @@ voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
H5Sclose(memspace);
}
RGBColor random_color() {
return {int(prn()*255), int(prn()*255), int(prn()*255)};
}
} // namespace openmc

74
src/plot_header.F90
View file

@ -1,74 +0,0 @@
module plot_header
use, intrinsic :: ISO_C_BINDING
use constants
use dict_header, only: DictIntInt
implicit none
!===============================================================================
! ObjectColor holds color information for plotted objects
!===============================================================================
type ObjectColor
integer :: rgb(3)
end type ObjectColor
!===============================================================================
! PLOTSLICE holds plot information
!===============================================================================
type ObjectPlot
integer :: id ! Unique ID
character(MAX_LINE_LEN) :: path_plot ! path for plot file
integer :: type ! Type
integer :: color_by ! quantity to color regions by
real(8) :: origin(3) ! xyz center of plot location
real(8) :: width(3) ! xyz widths of plot
integer :: basis ! direction of plot slice
integer :: pixels(3) ! pixel width/height of plot slice
integer :: meshlines_width ! pixel width of meshlines
integer :: level ! universe depth to plot the cells of
integer :: index_meshlines_mesh = -1 ! index of mesh to plot
type(ObjectColor) :: meshlines_color ! Color for meshlines
type(ObjectColor) :: not_found ! color for positions where no cell found
type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats
end type ObjectPlot
! Plot type
integer, parameter :: PLOT_TYPE_SLICE = 1
integer, parameter :: PLOT_TYPE_VOXEL = 2
! Plot level
integer, parameter :: PLOT_LEVEL_LOWEST = -1
! Plot basis plane
integer, parameter :: PLOT_BASIS_XY = 1
integer, parameter :: PLOT_BASIS_XZ = 2
integer, parameter :: PLOT_BASIS_YZ = 3
! Indicate whether color refers to unique cell or unique material
integer, parameter :: PLOT_COLOR_CELLS = 1
integer, parameter :: PLOT_COLOR_MATS = 2
integer(C_INT32_T), bind(C) :: n_plots ! # of plots
type(ObjectPlot), allocatable, target :: plots(:)
! Dictionary that maps user IDs to indices in 'plots'
type(DictIntInt) :: plot_dict
contains
!===============================================================================
! FREE_MEMORY_PLOT deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_plot()
n_plots = 0
if (allocated(plots)) deallocate(plots)
call plot_dict % clear()
end subroutine free_memory_plot
end module plot_header

67
src/progress_bar.cpp Normal file
View file

@ -0,0 +1,67 @@
#include "openmc/progress_bar.h"
#include <sstream>
#include <iostream>
#include <iomanip>
#ifdef UNIX
#include <unistd.h>
#endif
#define BAR_WIDTH 72
bool is_terminal() {
#ifdef UNIX
return isatty(STDOUT_FILENO) != 0;
#else
return false;
#endif
}
ProgressBar::ProgressBar() {
// initialize bar
set_value(0.0);
}
void
ProgressBar::set_value(double val) {
if (!is_terminal()) return;
// set the bar percentage
if (val >= 100.0) {
bar.append("100");
} else if (val <= 0.0) {
bar.append(" 0");
} else {
std::stringstream ss;
ss << std::setfill(' ') << std::setw(3) << (int)val;
bar.append(ss.str());
}
bar.append("% |");
// remaining width of the bar
int remaining_width = BAR_WIDTH - bar.size() - 2;
// set the bar width
if (val >= 100.0) {
bar.append(remaining_width, '=');
} else if (val < 0.0) {
bar.append(remaining_width, ' ');
} else {
int width = (int)((double)remaining_width*val/100);
bar.append(width, '=');
bar.append(1, '>');
bar.append(remaining_width-width-1, ' ');
}
bar.append("|+");
// write the bar
std::cout << '\r' << bar << std::flush;
if (val >= 100.0) { std::cout << "\n"; }
// reset the bar value
bar = "";
}

102
src/progress_header.F90
View file

@ -1,102 +0,0 @@
module progress_header
use, intrinsic :: ISO_FORTRAN_ENV, only: OUTPUT_UNIT
implicit none
#ifdef UNIX
interface
function check_isatty(fd) bind(C, name = 'isatty')
use, intrinsic :: ISO_C_BINDING, only: c_int
integer(c_int) :: check_isatty
integer(c_int), value :: fd
end function check_isatty
end interface
#endif
!===============================================================================
! PROGRESSBAR
!===============================================================================
type ProgressBar
private
character(len=72) :: bar="???% | " // &
" |"
contains
procedure :: set_value => bar_set_value
end type ProgressBar
contains
!===============================================================================
! IS_TERMINAL checks if output is currently being output to a terminal. Relies
! on a POSIX implementation of isatty, and defaults to false if that is not
! available
!===============================================================================
function is_terminal() result(istty)
logical :: istty
istty = .true.
#ifdef UNIX
if (check_isatty(1) == 0) istty = .false.
#else
istty = .false.
#endif
end function is_terminal
!===============================================================================
! BAR_SET_VALUE prints the progress bar without advancing. The value is
! specified as percent completion, from 0 to 100. If the value is ever set to
! 100 or above, the bar is set to 100 and a newline is written.
!===============================================================================
subroutine bar_set_value(self, val)
class(ProgressBar), intent(inout) :: self
real(8), intent(in) :: val
integer :: i
if (.not. is_terminal()) return
! set the percentage
if (val >= 100.) then
write(self % bar(1:3), "(I3)") 100
else if (val <= 0.) then
write(self % bar(1:3), "(I3)") 0
else
write(self % bar(1:3), "(I3)") int(val)
end if
! set the bar width
if (val >= 100.) then
do i=1,65
write(self % bar(i+6:i+6), '(A)') '='
end do
else
do i=1,int(dble(65)*val/100.)
write(self % bar(i+6:i+6), '(A)') '='
end do
end if
write(OUTPUT_UNIT, '(A1,A1,A72)', ADVANCE='no') '+', char(13), self % bar
flush(OUTPUT_UNIT)
if (val >= 100.) then
! make new line
write(OUTPUT_UNIT, "(A)") ""
flush(OUTPUT_UNIT)
! reset the bar in case we want to use this instance again
self % bar = "???% | " // &
" |"
end if
end subroutine bar_set_value
end module progress_header

3
src/settings.cpp
View file

@ -72,7 +72,8 @@ std::string path_statepoint;
int32_t index_entropy_mesh {-1};
int32_t index_ufs_mesh {-1};
int32_t index_cmfd_mesh {-1};
int32_t n_batches;
int32_t n_inactive {0};
int32_t gen_per_batch {1};

10
src/string_functions.cpp
View file

@ -1,4 +1,5 @@
#include "openmc/string_functions.h"
#include <sstream>
namespace openmc {
@ -27,4 +28,13 @@ void to_lower(std::string& str)
for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
}
int word_count(std::string const& str)
{
std::stringstream stream(str);
std::string dum;
int count = 0;
while (stream >> dum) {count++;}
return count;
}
} // namespace openmc

20
src/surface.cpp
View file

@ -8,7 +8,7 @@
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/xml_interface.h"
#include "openmc/string_functions.h"
namespace openmc {
@ -35,24 +35,6 @@ std::map<int, int> surface_map;
// Helper functions for reading the "coeffs" node of an XML surface element
//==============================================================================
int word_count(const std::string &text)
{
bool in_word = false;
int count {0};
for (auto c = text.begin(); c != text.end(); c++) {
if (std::isspace(*c)) {
if (in_word) {
in_word = false;
count++;
}
} else {
in_word = true;
}
}
if (in_word) count++;
return count;
}
void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1)
{
// Check the given number of coefficients.