From 1bc0ca274ca7dd148293467059f4dc115d0742ea Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 15 Aug 2018 16:06:05 -0400 Subject: [PATCH 1/6] Fix pointer * placement style --- src/cell.cpp | 32 ++++++++++++++++---------------- src/surface.cpp | 42 +++++++++++++++++++++--------------------- src/surface.h | 38 +++++++++++++++++++------------------- 3 files changed, 56 insertions(+), 56 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 798935358f..554ca74f1b 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -429,7 +429,7 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const //============================================================================== extern "C" void -read_cells(pugi::xml_node *node) +read_cells(pugi::xml_node* node) { // Count the number of cells. for (pugi::xml_node cell_node: node->children("cell")) {n_cells++;} @@ -467,35 +467,35 @@ read_cells(pugi::xml_node *node) extern "C" { Cell* cell_pointer(int32_t cell_ind) {return global_cells[cell_ind];} - int32_t cell_id(Cell *c) {return c->id;} + int32_t cell_id(Cell* c) {return c->id;} - void cell_set_id(Cell *c, int32_t id) {c->id = id;} + void cell_set_id(Cell* c, int32_t id) {c->id = id;} - int cell_type(Cell *c) {return c->type;} + int cell_type(Cell* c) {return c->type;} - void cell_set_type(Cell *c, int type) {c->type = type;} + void cell_set_type(Cell* c, int type) {c->type = type;} - int32_t cell_universe(Cell *c) {return c->universe;} + int32_t cell_universe(Cell* c) {return c->universe;} - void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;} + void cell_set_universe(Cell* c, int32_t universe) {c->universe = universe;} - int32_t cell_fill(Cell *c) {return c->fill;} + int32_t cell_fill(Cell* c) {return c->fill;} - int32_t* cell_fill_ptr(Cell *c) {return &c->fill;} + int32_t* cell_fill_ptr(Cell* c) {return &c->fill;} - int32_t cell_n_instances(Cell *c) {return c->n_instances;} + int32_t cell_n_instances(Cell* c) {return c->n_instances;} - bool cell_simple(Cell *c) {return c->simple;} + bool cell_simple(Cell* c) {return c->simple;} - bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface) + bool cell_contains(Cell* c, double xyz[3], double uvw[3], int32_t on_surface) { Position r {xyz}; Direction u {uvw}; return c->contains(r, u, on_surface); } - void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface, - double *min_dist, int32_t *i_surf) + void cell_distance(Cell* c, double xyz[3], double uvw[3], int32_t on_surface, + double* min_dist, int32_t* i_surf) { Position r {xyz}; Direction u {uvw}; @@ -504,9 +504,9 @@ extern "C" { *i_surf = out.second; } - int32_t cell_offset(Cell *c, int map) {return c->offset[map];} + int32_t cell_offset(Cell* c, int map) {return c->offset[map];} - void cell_to_hdf5(Cell *c, hid_t group) {c->to_hdf5(group);} + void cell_to_hdf5(Cell* c, hid_t group) {c->to_hdf5(group);} void extend_cells_c(int32_t n) { diff --git a/src/surface.cpp b/src/surface.cpp index 261a8d7d24..99cfd82b05 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -27,7 +27,7 @@ extern "C" const int BC_PERIODIC {3}; int32_t n_surfaces; -Surface **surfaces_c; +Surface** surfaces_c; std::map surface_map; @@ -298,8 +298,8 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r, - Direction& u) const +bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other, + Position& r, Direction& u) const { Direction other_n = other->normal(r); if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) { @@ -359,8 +359,8 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r, - Direction& u) const +bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, + Position& r, Direction& u) const { Direction other_n = other->normal(r); if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) { @@ -421,8 +421,8 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r, - Direction& u) const +bool SurfaceZPlane::periodic_translate(const PeriodicSurface* other, + Position& r, Direction& u) const { // Assume the other plane is aligned along z. Just change the z coord. r.z = z0; @@ -478,7 +478,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r, +bool SurfacePlane::periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const { // This function assumes the other plane shares this plane's normal direction. @@ -1023,7 +1023,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const //============================================================================== extern "C" void -read_surfaces(pugi::xml_node *node) +read_surfaces(pugi::xml_node* node) { // Count the number of surfaces. for (pugi::xml_node surf_node: node->children("surface")) {n_surfaces++;} @@ -1106,8 +1106,8 @@ read_surfaces(pugi::xml_node *node) for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { if (surfaces_c[i_surf]->bc == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface *surf_base = surfaces_c[i_surf]; - PeriodicSurface *surf = dynamic_cast(surf_base); + Surface* surf_base = surfaces_c[i_surf]; + PeriodicSurface* surf = dynamic_cast(surf_base); // Make sure this surface inherits from PeriodicSurface. if (!surf) { @@ -1151,12 +1151,12 @@ read_surfaces(pugi::xml_node *node) for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { if (surfaces_c[i_surf]->bc == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface *surf_base = surfaces_c[i_surf]; - PeriodicSurface *surf = dynamic_cast(surf_base); + Surface* surf_base = surfaces_c[i_surf]; + PeriodicSurface* surf = dynamic_cast(surf_base); // Also try downcasting to the SurfacePlane type (which must be handled // differently). - SurfacePlane *surf_p = dynamic_cast(surf); + SurfacePlane* surf_p = dynamic_cast(surf); if (!surf_p) { // This is not a SurfacePlane. @@ -1215,11 +1215,11 @@ read_surfaces(pugi::xml_node *node) extern "C" { Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];} - int surface_id(Surface *surf) {return surf->id;} + int surface_id(Surface* surf) {return surf->id;} - int surface_bc(Surface *surf) {return surf->bc;} + int surface_bc(Surface* surf) {return surf->bc;} - void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) + void surface_reflect(Surface* surf, double xyz[3], double uvw[3]) { Position r {xyz}; Direction u {uvw}; @@ -1230,7 +1230,7 @@ extern "C" { uvw[2] = u.z; } - void surface_normal(Surface *surf, double xyz[3], double uvw[3]) + void surface_normal(Surface* surf, double xyz[3], double uvw[3]) { Position r {xyz}; Direction u = surf->normal(r); @@ -1239,12 +1239,12 @@ extern "C" { uvw[2] = u.z; } - void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);} + void surface_to_hdf5(Surface* surf, hid_t group) {surf->to_hdf5(group);} - int surface_i_periodic(PeriodicSurface *surf) {return surf->i_periodic;} + int surface_i_periodic(PeriodicSurface* surf) {return surf->i_periodic;} bool - surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], + surface_periodic(PeriodicSurface* surf, PeriodicSurface* other, double xyz[3], double uvw[3]) { Position r {xyz}; diff --git a/src/surface.h b/src/surface.h index d4e6c0230a..d43b24ae06 100644 --- a/src/surface.h +++ b/src/surface.h @@ -31,7 +31,7 @@ extern "C" const int BC_PERIODIC; extern "C" int32_t n_surfaces; class Surface; -extern Surface **surfaces_c; +extern Surface** surfaces_c; extern std::map surface_map; @@ -131,7 +131,7 @@ public: //! periodicity. //! \return true if this surface and its partner make a rotationally-periodic //! boundary condition. - virtual bool periodic_translate(const PeriodicSurface *other, Position& r, + virtual bool periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const = 0; //! Get the bounding box for this surface. @@ -153,8 +153,8 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) - const; + bool periodic_translate(const PeriodicSurface* other, Position& r, + Direction& u) const; BoundingBox bounding_box() const; }; @@ -173,8 +173,8 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) - const; + bool periodic_translate(const PeriodicSurface* other, Position& r, + Direction& u) const; BoundingBox bounding_box() const; }; @@ -193,8 +193,8 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) - const; + bool periodic_translate(const PeriodicSurface* other, Position& r, + Direction& u) const; BoundingBox bounding_box() const; }; @@ -213,8 +213,8 @@ public: double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) - const; + bool periodic_translate(const PeriodicSurface* other, Position& r, + Direction& u) const; BoundingBox bounding_box() const; }; @@ -368,16 +368,16 @@ public: extern "C" { Surface* surface_pointer(int surf_ind); - int surface_id(Surface *surf); - int surface_bc(Surface *surf); - bool surface_sense(Surface *surf, double xyz[3], double uvw[3]); - void surface_reflect(Surface *surf, double xyz[3], double uvw[3]); - double surface_distance(Surface *surf, double xyz[3], double uvw[3], + int surface_id(Surface* surf); + int surface_bc(Surface* surf); + bool surface_sense(Surface* surf, double xyz[3], double uvw[3]); + void surface_reflect(Surface* surf, double xyz[3], double uvw[3]); + double surface_distance(Surface* surf, double xyz[3], double uvw[3], bool coincident); - void surface_normal(Surface *surf, double xyz[3], double uvw[3]); - void surface_to_hdf5(Surface *surf, hid_t group); - int surface_i_periodic(PeriodicSurface *surf); - bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, + void surface_normal(Surface* surf, double xyz[3], double uvw[3]); + void surface_to_hdf5(Surface* surf, hid_t group); + int surface_i_periodic(PeriodicSurface* surf); + bool surface_periodic(PeriodicSurface* surf, PeriodicSurface* other, double xyz[3], double uvw[3]); void free_memory_surfaces_c(); } From e5d154b95595c1506983170977acf7289f909dd0 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 15 Aug 2018 16:08:32 -0400 Subject: [PATCH 2/6] Add basic C++ Material class --- CMakeLists.txt | 1 + src/api.F90 | 2 +- src/geometry_header.F90 | 11 ++--- src/input_xml.F90 | 48 +++++++++---------- src/material.cpp | 66 +++++++++++++++++++++++++++ src/material.h | 33 ++++++++++++++ src/material_header.F90 | 57 +++++++++++++++++++++-- src/plot.F90 | 2 +- src/summary.F90 | 6 +-- src/surface_header.F90 | 5 +- src/tallies/tally.F90 | 60 ++++++++++++------------ src/tallies/tally_filter_material.F90 | 4 +- src/volume_calc.F90 | 2 +- 13 files changed, 222 insertions(+), 75 deletions(-) create mode 100644 src/material.cpp create mode 100644 src/material.h diff --git a/CMakeLists.txt b/CMakeLists.txt index dc9a1231ba..6a2e5aee9b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -397,6 +397,7 @@ add_library(libopenmc SHARED src/geometry_aux.cpp src/hdf5_interface.cpp src/lattice.cpp + src/material.cpp src/math_functions.cpp src/message_passing.cpp src/mgxs.cpp diff --git a/src/api.F90 b/src/api.F90 index cc9120c7b9..6a6c10671a 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -214,7 +214,7 @@ contains if (p % material == MATERIAL_VOID) then id = 0 else - id = materials(p % material) % id + id = materials(p % material) % id() end if end if instance = p % cell_instance - 1 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 71c4581055..5a2dca282f 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -16,11 +16,11 @@ module geometry_header implicit none interface - function cell_pointer_c(cell_ind) bind(C, name='cell_pointer') result(ptr) + function cell_pointer(cell_ind) bind(C) result(ptr) import C_PTR, C_INT32_T integer(C_INT32_T), intent(in), value :: cell_ind type(C_PTR) :: ptr - end function cell_pointer_c + end function cell_pointer function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id) import C_PTR, C_INT32_T @@ -110,12 +110,11 @@ module geometry_header integer(HID_T), intent(in), value :: group end subroutine cell_to_hdf5_c - function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') & - result(ptr) + function lattice_pointer(lat_ind) bind(C) result(ptr) import C_PTR, C_INT32_T integer(C_INT32_T), intent(in), value :: lat_ind type(C_PTR) :: ptr - end function lattice_pointer_c + end function lattice_pointer function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id) import C_PTR, C_INT32_T @@ -571,7 +570,7 @@ contains ! Extend the C++ cells array and get pointers to the C++ objects call extend_cells_c(n) do i = n_cells - n, n_cells - cells(i) % ptr = cell_pointer_c(i - 1) + cells(i) % ptr = cell_pointer(i - 1) end do err = 0 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index da33a6627f..aef9a14c07 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -87,6 +87,11 @@ module input_xml type(C_PTR) :: node_ptr end subroutine read_settings + subroutine read_materials(node_ptr) bind(C) + import C_PTR + type(C_PTR) :: node_ptr + end subroutine read_materials + function find_root_universe() bind(C) result(root) import C_INT32_T integer(C_INT32_T) :: root @@ -1051,7 +1056,7 @@ contains allocate(surfaces(n_surfaces)) do i = 1, n_surfaces - surfaces(i) % ptr = surface_pointer_c(i - 1); + surfaces(i) % ptr = surface_pointer(i - 1); if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true. @@ -1095,7 +1100,7 @@ contains do i = 1, n_cells c => cells(i) - c % ptr = cell_pointer_c(i - 1) + c % ptr = cell_pointer(i - 1) ! Initialize distribcell instances and distribcell index c % distribcell_index = NONE @@ -1309,7 +1314,7 @@ contains RECT_LATTICES: do i = 1, n_rlats allocate(RectLattice::lattices(i) % obj) lat => lattices(i) % obj - lat % ptr = lattice_pointer_c(i - 1) + lat % ptr = lattice_pointer(i - 1) select type(lat) type is (RectLattice) @@ -1325,7 +1330,7 @@ contains HEX_LATTICES: do i = 1, n_hlats allocate(HexLattice::lattices(n_rlats + i) % obj) lat => lattices(n_rlats + i) % obj - lat % ptr = lattice_pointer_c(n_rlats + i - 1) + lat % ptr = lattice_pointer(n_rlats + i - 1) select type (lat) type is (HexLattice) @@ -1540,6 +1545,8 @@ contains call doc % load_file(filename) root = doc % document_element() + call read_materials(root % ptr) + ! Get pointer to list of XML call get_node_list(root, "material", node_mat_list) @@ -1556,16 +1563,11 @@ contains do i = 1, n_materials mat => materials(i) + mat % ptr = material_pointer(i - 1) + ! Get pointer to i-th material node node_mat = node_mat_list(i) - ! Copy material id - if (check_for_node(node_mat, "id")) then - call get_node_value(node_mat, "id", mat % id) - else - call fatal_error("Must specify id of material in materials XML file") - end if - ! Check if material is depletable if (check_for_node(node_mat, "depletable")) then call get_node_value(node_mat, "depletable", temp_str) @@ -1574,9 +1576,9 @@ contains end if ! Check to make sure 'id' hasn't been used - if (material_dict % has(mat % id)) then + if (material_dict % has(mat % id())) then call fatal_error("Two or more materials use the same unique ID: " & - // to_str(mat % id)) + // to_str(mat % id())) end if ! Copy material name @@ -1596,7 +1598,7 @@ contains node_dens = node_mat % child("density") else call fatal_error("Must specify density element in material " & - // trim(to_str(mat % id))) + // trim(to_str(mat % id()))) end if ! Copy units @@ -1626,7 +1628,7 @@ contains sum_density = .false. if (val <= ZERO) then call fatal_error("Need to specify a positive density on material " & - // trim(to_str(mat % id)) // ".") + // trim(to_str(mat % id())) // ".") end if ! Adjust material density based on specified units @@ -1641,7 +1643,7 @@ contains mat % density = 1.0e-24_8 * val case default call fatal_error("Unkwown units '" // trim(units) & - // "' specified on material " // trim(to_str(mat % id))) + // "' specified on material " // trim(to_str(mat % id()))) end select end if @@ -1650,7 +1652,7 @@ contains if (size(node_ele_list) > 0) then call fatal_error("Unable to add an element to material " & - // trim(to_str(mat % id)) // " since the element option has & + // trim(to_str(mat % id())) // " since the element option has & &been removed from the xml input. Elements can only be added via & &the Python API, which will expand elements into their natural & &nuclides.") @@ -1663,7 +1665,7 @@ contains if (.not. check_for_node(node_mat, "nuclide") .and. & .not. check_for_node(node_mat, "macroscopic")) then call fatal_error("No macroscopic data or nuclides specified on & - &material " // trim(to_str(mat % id))) + &material " // trim(to_str(mat % id()))) end if ! Create list of macroscopic x/s based on those specified, just treat @@ -1677,7 +1679,7 @@ contains & mode!") else if (size(node_macro_list) > 1) then call fatal_error("Only one macroscopic object permitted per material, " & - // trim(to_str(mat % id))) + // trim(to_str(mat % id()))) else if (size(node_macro_list) == 1) then node_nuc = node_macro_list(1) @@ -1685,7 +1687,7 @@ contains ! Check for empty name on nuclide if (.not. check_for_node(node_nuc, "name")) then call fatal_error("No name specified on macroscopic data in material " & - // trim(to_str(mat % id))) + // trim(to_str(mat % id()))) end if ! store nuclide name @@ -1715,7 +1717,7 @@ contains ! Check for empty name on nuclide if (.not. check_for_node(node_nuc, "name")) then call fatal_error("No name specified on nuclide in material " & - // trim(to_str(mat % id))) + // trim(to_str(mat % id()))) end if ! store nuclide name @@ -1854,7 +1856,7 @@ contains if (.not. (all(mat % atom_density >= ZERO) .or. & all(mat % atom_density <= ZERO))) then call fatal_error("Cannot mix atom and weight percents in material " & - // to_str(mat % id)) + // to_str(mat % id())) end if ! Determine density if it is a sum value @@ -1935,7 +1937,7 @@ contains end if ! Add material to dictionary - call material_dict % set(mat % id, i) + call material_dict % set(mat % id(), i) end do ! Set total number of nuclides and S(a,b) tables diff --git a/src/material.cpp b/src/material.cpp new file mode 100644 index 0000000000..49f0e24158 --- /dev/null +++ b/src/material.cpp @@ -0,0 +1,66 @@ +#include "material.h" + +#include "error.h" +#include "xml_interface.h" + +//TODO: remove this include +#include + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +std::vector global_materials; + +//============================================================================== +// Material implementation +//============================================================================== + +Material::Material(pugi::xml_node material_node) +{ + if (check_for_node(material_node, "id")) { + id = stoi(get_node_value(material_node, "id")); + } else { + fatal_error("Must specify id of material in materials XML file."); + } + + std::cout << "Reading material with id " << id << std::endl; +} + +//============================================================================== +// Non-method functions +//============================================================================== + +extern "C" void +read_materials(pugi::xml_node* node) +{ + // Loop over XML material elements and populate the array. + for (pugi::xml_node material_node: node->children("material")) { + global_materials.push_back(new Material(material_node)); + } +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + Material* material_pointer(int32_t indx) {return global_materials[indx];} + + int32_t material_id(Material* mat) {return mat->id;} + + void material_set_id(Material* mat, int32_t id) {mat->id = id;} + + void extend_materials_c(int32_t n) + { + global_materials.reserve(global_materials.size() + n); + for (int32_t i = 0; i < n; i++) { + global_materials.push_back(new Material()); + } + } +} + +} // namespace openmc diff --git a/src/material.h b/src/material.h new file mode 100644 index 0000000000..f2e3df2fb4 --- /dev/null +++ b/src/material.h @@ -0,0 +1,33 @@ +#ifndef OPENMC_CELL_H +#define OPENMC_CELL_H + +#include + +#include "pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +class Material; +extern std::vector global_materials; + +//============================================================================== +//! A substance with constituent nuclides and thermal scattering data +//============================================================================== + +class Material +{ +public: + int32_t id; //!< Unique ID + + Material() {}; + + explicit Material(pugi::xml_node material_node); +}; + +} // namespace openmc +#endif // OPENMC_CELL_H diff --git a/src/material_header.F90 b/src/material_header.F90 index 03ecf69d9f..5bfbadfcf9 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -27,13 +27,39 @@ module material_header public :: openmc_material_set_density public :: openmc_material_set_densities public :: openmc_material_set_id + public :: material_pointer + + interface + function material_pointer(mat_ind) bind(C) result(ptr) + import C_PTR, C_INT32_T + integer(C_INT32_T), intent(in), value :: mat_ind + type(C_PTR) :: ptr + end function material_pointer + + function material_id_c(mat_ptr) bind(C, name='material_id') result(id) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: mat_ptr + integer(C_INT32_T) :: id + end function material_id_c + + subroutine material_set_id_c(mat_ptr, id) bind(C, name='material_set_id') + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: mat_ptr + integer(C_INT32_T), intent(in), value :: id + end subroutine material_set_id_c + + subroutine extend_materials_c(n) bind(C) + import C_INT32_t + integer(C_INT32_T), intent(in), value :: n + end subroutine extend_materials_c + end interface !=============================================================================== ! MATERIAL describes a material by its constituent nuclides !=============================================================================== type, public :: Material - integer :: id ! unique identifier + type(C_PTR) :: ptr character(len=104) :: name = "" ! User-defined name integer :: n_nuclides = 0 ! number of nuclides integer, allocatable :: nuclide(:) ! index in nuclides array @@ -67,6 +93,8 @@ module material_header logical, allocatable :: p0(:) contains + procedure :: id => material_id + procedure :: set_id => material_set_id procedure :: set_density => material_set_density procedure :: init_nuclide_index => material_init_nuclide_index procedure :: assign_sab_tables => material_assign_sab_tables @@ -88,6 +116,18 @@ contains ! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm. !=============================================================================== + function material_id(this) result(id) + class(Material), intent(in) :: this + integer(C_INT32_T) :: id + id = material_id_c(this % ptr) + end function material_id + + subroutine material_set_id(this, id) + class(Material), intent(in) :: this + integer(C_INT32_T), intent(in) :: id + call material_set_id_c(this % ptr, id) + end subroutine material_set_id + function material_set_density(this, density) result(err) class(Material), intent(inout) :: this real(8), intent(in) :: density @@ -185,7 +225,7 @@ contains if (.not. found) then call fatal_error("S(a,b) table " // trim(this % & sab_names(k)) // " did not match any nuclide on material " & - // trim(to_str(this % id))) + // trim(to_str(this % id()))) end if end do ASSIGN_SAB @@ -195,7 +235,7 @@ contains if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then call fatal_error(trim( & nuclides(this % nuclide(i_sab_nuclides % data(j))) % name) & - // " in material " // trim(to_str(this % id)) // " was found & + // " in material " // trim(to_str(this % id())) // " was found & &in multiple S(a,b) tables. Each nuclide can only appear in & &one S(a,b) table per material.") end if @@ -460,6 +500,7 @@ contains integer(C_INT32_T), optional, intent(out) :: index_end integer(C_INT) :: err + integer :: i type(Material), allocatable :: temp(:) ! temporary materials array if (n_materials == 0) then @@ -481,6 +522,12 @@ contains if (present(index_end)) index_end = n_materials + n n_materials = n_materials + n + ! Extend the C++ materials array and get pointers to the C++ objects + call extend_materials_c(n) + do i = n_materials - n, n_materials + materials(i) % ptr = material_pointer(i - 1) + end do + err = 0 end function openmc_extend_materials @@ -606,7 +653,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= size(materials)) then - id = materials(index) % id + id = materials(index) % id() err = 0 else err = E_OUT_OF_BOUNDS @@ -622,7 +669,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= n_materials) then - materials(index) % id = id + call materials(index) % set_id(id) call material_dict % set(id, index) err = 0 else diff --git a/src/plot.F90 b/src/plot.F90 index d8db97ec6f..2fd898206c 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -95,7 +95,7 @@ contains id = -1 else rgb = pl % colors(p % material) % rgb - id = materials(p % material) % id + id = materials(p % material) % id() end if end associate else if (pl % color_by == PLOT_COLOR_CELLS) then diff --git a/src/summary.F90 b/src/summary.F90 index 062b4efc24..507642f713 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -203,7 +203,7 @@ contains call write_dataset(cell_group, "material", MATERIAL_VOID) else call write_dataset(cell_group, "material", & - materials(c % material(1)) % id) + materials(c % material(1)) % id()) end if else allocate(cell_materials(size(c % material))) @@ -211,7 +211,7 @@ contains if (c % material(j) == MATERIAL_VOID) then cell_materials(j) = MATERIAL_VOID else - cell_materials(j) = materials(c % material(j)) % id + cell_materials(j) = materials(c % material(j)) % id() end if end do call write_dataset(cell_group, "material", cell_materials) @@ -333,7 +333,7 @@ contains do i = 1, n_materials m => materials(i) material_group = create_group(materials_group, "material " // & - trim(to_str(m%id))) + trim(to_str(m%id()))) if (m % depletable) then call write_attribute(material_group, "depletable", 1) diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 916f9c11c9..4a5834ec2c 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -8,13 +8,12 @@ module surface_header implicit none interface - pure function surface_pointer_c(surf_ind) & - bind(C, name='surface_pointer') result(ptr) + pure function surface_pointer(surf_ind) bind(C) result(ptr) use ISO_C_BINDING implicit none integer(C_INT), intent(in), value :: surf_ind type(C_PTR) :: ptr - end function surface_pointer_c + end function surface_pointer pure function surface_id_c(surf_ptr) bind(C, name='surface_id') result(id) use ISO_C_BINDING diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 2c660ce87b..94abb42c05 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -3043,7 +3043,7 @@ contains case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & SCORE_NU_FISSION) - if (materials(p % material) % id == deriv % diff_material) then + if (materials(p % material) % id() == deriv % diff_material) then score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else @@ -3064,7 +3064,7 @@ contains case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & SCORE_NU_FISSION) - if (materials(p % material) % id == deriv % diff_material) then + if (materials(p % material) % id() == deriv % diff_material) then score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else @@ -3106,7 +3106,7 @@ contains case (SCORE_TOTAL, SCORE_SCATTER, SCORE_ABSORPTION, SCORE_FISSION, & SCORE_NU_FISSION) - if (materials(p % material) % id == deriv % diff_material & + if (materials(p % material) % id() == deriv % diff_material & .and. p % event_nuclide == deriv % diff_nuclide) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. @@ -3128,7 +3128,7 @@ contains case (ESTIMATOR_COLLISION) scoring_diff_nuclide = & - (materials(p % material) % id == deriv % diff_material) & + (materials(p % material) % id() == deriv % diff_material) & .and. (i_nuclide == deriv % diff_nuclide) select case (score_bin) @@ -3138,7 +3138,7 @@ contains case (SCORE_TOTAL) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % total /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % total & @@ -3152,7 +3152,7 @@ contains case (SCORE_SCATTER) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % total - material_xs % absorption /= ZERO) then score = score * (flux_deriv & + (micro_xs(deriv % diff_nuclide) % total & @@ -3168,7 +3168,7 @@ contains case (SCORE_ABSORPTION) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % absorption /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % absorption & @@ -3182,7 +3182,7 @@ contains case (SCORE_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % fission /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % fission & @@ -3196,7 +3196,7 @@ contains case (SCORE_NU_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % nu_fission /= ZERO) then score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % nu_fission & @@ -3242,7 +3242,7 @@ contains score = score * flux_deriv case (SCORE_TOTAL) - if (materials(p % material) % id == deriv % diff_material .and. & + if (materials(p % material) % id() == deriv % diff_material .and. & micro_xs(p % event_nuclide) % total > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. @@ -3267,7 +3267,7 @@ contains end if case (SCORE_SCATTER) - if (materials(p % material) % id == deriv % diff_material .and. & + if (materials(p % material) % id() == deriv % diff_material .and. & (micro_xs(p % event_nuclide) % total & - micro_xs(p % event_nuclide) % absorption) > ZERO) then associate(mat => materials(p % material)) @@ -3295,7 +3295,7 @@ contains end if case (SCORE_ABSORPTION) - if (materials(p % material) % id == deriv % diff_material .and. & + if (materials(p % material) % id() == deriv % diff_material .and. & micro_xs(p % event_nuclide) % absorption > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. @@ -3320,7 +3320,7 @@ contains end if case (SCORE_FISSION) - if (materials(p % material) % id == deriv % diff_material .and. & + if (materials(p % material) % id() == deriv % diff_material .and. & micro_xs(p % event_nuclide) % fission > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. @@ -3345,7 +3345,7 @@ contains end if case (SCORE_NU_FISSION) - if (materials(p % material) % id == deriv % diff_material .and. & + if (materials(p % material) % id() == deriv % diff_material .and. & micro_xs(p % event_nuclide) % nu_fission > ZERO) then associate(mat => materials(p % material)) ! Search for the index of the perturbed nuclide. @@ -3385,7 +3385,7 @@ contains case (SCORE_TOTAL) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % total > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) @@ -3404,7 +3404,7 @@ contains end associate score = score * (flux_deriv & + cum_dsig / material_xs % total) - else if (materials(p % material) % id == deriv % diff_material & + else if (materials(p % material) % id() == deriv % diff_material & .and. material_xs % total > ZERO) then dsig_t = ZERO associate (nuc => nuclides(i_nuclide)) @@ -3423,7 +3423,7 @@ contains case (SCORE_SCATTER) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & (material_xs % total - material_xs % absorption) > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) @@ -3444,7 +3444,7 @@ contains end associate score = score * (flux_deriv + cum_dsig & / (material_xs % total - material_xs % absorption)) - else if ( materials(p % material) % id == deriv % diff_material & + else if ( materials(p % material) % id() == deriv % diff_material & .and. (material_xs % total - material_xs % absorption) > ZERO)& then dsig_t = ZERO @@ -3466,7 +3466,7 @@ contains case (SCORE_ABSORPTION) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % absorption > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) @@ -3485,7 +3485,7 @@ contains end associate score = score * (flux_deriv & + cum_dsig / material_xs % absorption) - else if (materials(p % material) % id == deriv % diff_material & + else if (materials(p % material) % id() == deriv % diff_material & .and. material_xs % absorption > ZERO) then dsig_a = ZERO associate (nuc => nuclides(i_nuclide)) @@ -3504,7 +3504,7 @@ contains case (SCORE_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % fission > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) @@ -3523,7 +3523,7 @@ contains end associate score = score * (flux_deriv & + cum_dsig / material_xs % fission) - else if (materials(p % material) % id == deriv % diff_material & + else if (materials(p % material) % id() == deriv % diff_material & .and. material_xs % fission > ZERO) then dsig_f = ZERO associate (nuc => nuclides(i_nuclide)) @@ -3542,7 +3542,7 @@ contains case (SCORE_NU_FISSION) if (i_nuclide == -1 .and. & - materials(p % material) % id == deriv % diff_material .and. & + materials(p % material) % id() == deriv % diff_material .and. & material_xs % nu_fission > ZERO) then cum_dsig = ZERO associate(mat => materials(p % material)) @@ -3563,7 +3563,7 @@ contains end associate score = score * (flux_deriv & + cum_dsig / material_xs % nu_fission) - else if (materials(p % material) % id == deriv % diff_material & + else if (materials(p % material) % id() == deriv % diff_material & .and. material_xs % nu_fission > ZERO) then dsig_f = ZERO associate (nuc => nuclides(i_nuclide)) @@ -3614,7 +3614,7 @@ contains case (DIFF_DENSITY) associate (mat => materials(p % material)) - if (mat % id == deriv % diff_material) then + if (mat % id() == deriv % diff_material) then ! phi is proportional to e^(-Sigma_tot * dist) ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist @@ -3625,7 +3625,7 @@ contains case (DIFF_NUCLIDE_DENSITY) associate (mat => materials(p % material)) - if (mat % id == deriv % diff_material) then + if (mat % id() == deriv % diff_material) then ! phi is proportional to e^(-Sigma_tot * dist) ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist @@ -3636,7 +3636,7 @@ contains case (DIFF_TEMPERATURE) associate (mat => materials(p % material)) - if (mat % id == deriv % diff_material) then + if (mat % id() == deriv % diff_material) then do l=1, mat % n_nuclides associate (nuc => nuclides(mat % nuclide(l))) if (nuc % mp_present .and. & @@ -3690,7 +3690,7 @@ contains case (DIFF_DENSITY) associate (mat => materials(p % material)) - if (mat % id == deriv % diff_material) then + if (mat % id() == deriv % diff_material) then ! phi is proportional to Sigma_s ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s ! (1 / phi) * (d_phi / d_rho) = 1 / rho @@ -3701,7 +3701,7 @@ contains case (DIFF_NUCLIDE_DENSITY) associate (mat => materials(p % material)) - if (mat % id == deriv % diff_material & + if (mat % id() == deriv % diff_material & .and. p % event_nuclide == deriv % diff_nuclide) then ! Find the index in this material for the diff_nuclide. do j = 1, mat % n_nuclides @@ -3722,7 +3722,7 @@ contains case (DIFF_TEMPERATURE) associate (mat => materials(p % material)) - if (mat % id == deriv % diff_material) then + if (mat % id() == deriv % diff_material) then do l=1, mat % n_nuclides associate (nuc => nuclides(mat % nuclide(l))) if (mat % nuclide(l) == p % event_nuclide .and. & diff --git a/src/tallies/tally_filter_material.F90 b/src/tallies/tally_filter_material.F90 index 443dc62855..51ab26594e 100644 --- a/src/tallies/tally_filter_material.F90 +++ b/src/tallies/tally_filter_material.F90 @@ -76,7 +76,7 @@ contains allocate(material_ids(size(this % materials))) do i = 1, size(this % materials) - material_ids(i) = materials(this % materials(i)) % id + material_ids(i) = materials(this % materials(i)) % id() end do call write_dataset(filter_group, "bins", material_ids) end subroutine to_statepoint_material @@ -110,7 +110,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Material " // to_str(materials(this % materials(bin)) % id) + label = "Material " // to_str(materials(this % materials(bin)) % id()) end function text_label_material !=============================================================================== diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 93087df899..e0c3e88594 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -196,7 +196,7 @@ contains i_material = p % material if (i_material /= MATERIAL_VOID) then do i_domain = 1, size(this % domain_id) - if (materials(i_material) % id == this % domain_id(i_domain)) then + if (materials(i_material) % id() == this % domain_id(i_domain)) then call check_hit(i_domain, i_material, indices, hits, n_mat) end if end do From 1f11eabc71264003e740e285809a317c34de5f04 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 15 Aug 2018 19:50:05 -0400 Subject: [PATCH 3/6] Move cell materials array to C++ --- docs/source/io_formats/summary.rst | 2 +- openmc/capi/cell.py | 9 +- src/api.F90 | 2 - src/cell.cpp | 99 +++++++++++++++++++-- src/constants.F90 | 2 +- src/constants.h | 3 +- src/geometry.F90 | 4 +- src/geometry_aux.cpp | 18 +++- src/geometry_aux.h | 2 +- src/geometry_header.F90 | 138 +++++++---------------------- src/input_xml.F90 | 99 ++------------------- src/material.cpp | 26 +++++- src/material.h | 8 +- src/material_header.F90 | 5 ++ src/mgxs_data.F90 | 2 +- src/summary.F90 | 6 +- src/xml_interface.h | 2 +- 17 files changed, 197 insertions(+), 230 deletions(-) diff --git a/docs/source/io_formats/summary.rst b/docs/source/io_formats/summary.rst index 1f30c69bb1..05c6fff9c0 100644 --- a/docs/source/io_formats/summary.rst +++ b/docs/source/io_formats/summary.rst @@ -34,7 +34,7 @@ The current version of the summary file format is 6.0. 'material', 'universe', or 'lattice'. - **material** (*int* or *int[]*) -- Unique ID of the material(s) assigned to the cell. This dataset is present only if fill_type is - set to 'normal'. The value '-1' signifies void material. The data + set to 'normal'. The value '-2' signifies void material. The data is an array if the cell uses distributed materials, otherwise it is a scalar. - **temperature** (*double[]*) -- Temperature of the cell in Kelvin. diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 258e773dce..75986d2988 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -106,9 +106,10 @@ class Cell(_FortranObjectWithID): if fill_type.value == 1: if n.value > 1: - return [Material(index=i) for i in indices[:n.value]] + #TODO: off-by-one + return [Material(index=i+1) for i in indices[:n.value]] else: - return Material(index=indices[0]) + return Material(index=indices[0]+1) else: raise NotImplementedError @@ -116,14 +117,14 @@ class Cell(_FortranObjectWithID): def fill(self, fill): if isinstance(fill, Iterable): n = len(fill) - indices = (c_int32*n)(*(m._index if m is not None else -1 + indices = (c_int32*n)(*(m._index if m is not None else -2 for m in fill)) _dll.openmc_cell_set_fill(self._index, 1, n, indices) elif isinstance(fill, Material): indices = (c_int32*1)(fill._index) _dll.openmc_cell_set_fill(self._index, 1, 1, indices) elif fill is None: - indices = (c_int32*1)(-1) + indices = (c_int32*1)(-2) _dll.openmc_cell_set_fill(self._index, 1, 1, indices) def set_temperature(self, T, instance=None): diff --git a/src/api.F90 b/src/api.F90 index 6a6c10671a..60d94784de 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -37,8 +37,6 @@ module openmc_api public :: openmc_calculate_volumes public :: openmc_cell_filter_get_bins public :: openmc_cell_get_id - public :: openmc_cell_get_fill - public :: openmc_cell_set_fill public :: openmc_cell_set_id public :: openmc_cell_set_temperature public :: openmc_energy_filter_get_bins diff --git a/src/cell.cpp b/src/cell.cpp index 554ca74f1b..52b9a2e44f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -8,6 +8,8 @@ #include "error.h" #include "hdf5_interface.h" #include "lattice.h" +#include "material.h" +#include "openmc.h" #include "surface.h" #include "xml_interface.h" @@ -219,10 +221,25 @@ Cell::Cell(pugi::xml_node cell_node) fill = C_NONE; } + // Read the material element. There can be zero materials (filled with a + // universe), more than one material (distribmats), and some materials may + // be "void". if (check_for_node(cell_node, "material")) { - //TODO: read material ids. - material.push_back(C_NONE+1); - material.shrink_to_fit(); + std::vector mats + {get_node_array(cell_node, "material")}; + if (mats.size() > 0) { + material.reserve(mats.size()); + for (std::string mat : mats) { + if (mat.compare("void") == 0) { + material.push_back(MATERIAL_VOID); + } else { + material.push_back(std::stoi(mat)); + } + } + } else { + material.push_back(C_NONE); + material.shrink_to_fit(); + } } else { material.push_back(C_NONE); material.shrink_to_fit(); @@ -460,6 +477,65 @@ read_cells(pugi::xml_node* node) } } +//============================================================================== +// C-API functions +//============================================================================== + +extern "C" int +openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n) +{ + if (index >= 1 && index <= global_cells.size()) { + //TODO: off-by-one + Cell& c {*global_cells[index - 1]}; + *type = c.type; + if (c.type == FILL_MATERIAL) { + *indices = c.material.data(); + *n = c.material.size(); + } else { + *indices = &c.fill; + *n = 1; + } + } else { + strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + return 0; +} + +extern "C" int +openmc_cell_set_fill(int32_t index, int type, int32_t n, + const int32_t* indices) +{ + if (index >= 1 && index <= global_cells.size()) { + //TODO: off-by-one + Cell& c {*global_cells[index - 1]}; + if (type == FILL_MATERIAL) { + c.type = FILL_MATERIAL; + c.material.clear(); + for (int i = 0; i < n; i++) { + int i_mat = indices[i]; + if (i_mat == MATERIAL_VOID) { + c.material.push_back(MATERIAL_VOID); + } else if (i_mat >= 1 && i_mat <= global_materials.size()) { + //TODO: off-by-one + c.material.push_back(i_mat - 1); + } else { + strcpy(openmc_err_msg, "Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + } + } else if (type == FILL_UNIVERSE) { + c.type = FILL_UNIVERSE; + } else { + c.type = FILL_LATTICE; + } + } else { + strcpy(openmc_err_msg, "Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + return 0; +} + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -473,18 +549,23 @@ extern "C" { int cell_type(Cell* c) {return c->type;} - void cell_set_type(Cell* c, int type) {c->type = type;} - int32_t cell_universe(Cell* c) {return c->universe;} - void cell_set_universe(Cell* c, int32_t universe) {c->universe = universe;} - int32_t cell_fill(Cell* c) {return c->fill;} - int32_t* cell_fill_ptr(Cell* c) {return &c->fill;} - int32_t cell_n_instances(Cell* c) {return c->n_instances;} + int cell_material_size(Cell* c) {return c->material.size();} + + //TODO: off-by-one + int32_t cell_material(Cell* c, int i) + { + int32_t mat = c->material[i-1]; + if (mat == C_NONE) return F90_NONE; + if (mat == MATERIAL_VOID) return MATERIAL_VOID; + return mat + 1; + } + bool cell_simple(Cell* c) {return c->simple;} bool cell_contains(Cell* c, double xyz[3], double uvw[3], int32_t on_surface) diff --git a/src/constants.F90 b/src/constants.F90 index ab85b277a6..4b0f392fd7 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -129,7 +129,7 @@ module constants integer(C_INT), bind(C, name='FILL_LATTICE') :: FILL_LATTICE_C = FILL_LATTICE ! Void material - integer, parameter :: MATERIAL_VOID = -1 + integer, parameter :: MATERIAL_VOID = -2 ! Flag to say that the outside of a lattice is not defined integer, parameter :: NO_OUTER_UNIVERSE = -1 diff --git a/src/constants.h b/src/constants.h index b6bb82d375..72c72bad3f 100644 --- a/src/constants.h +++ b/src/constants.h @@ -112,7 +112,7 @@ constexpr char SUBSHELLS[][4] { // Void material // TODO: refactor and remove -constexpr int MATERIAL_VOID {-1}; +constexpr int MATERIAL_VOID {-2}; // ============================================================================ // CROSS SECTION RELATED CONSTANTS @@ -434,6 +434,7 @@ constexpr int DIFF_NUCLIDE_DENSITY {2}; constexpr int DIFF_TEMPERATURE {3}; constexpr int C_NONE {-1}; +constexpr int F90_NONE {0}; // Interpolation rules enum class Interpolation { diff --git a/src/geometry.F90 b/src/geometry.F90 index 5bf0012a38..81e123d174 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -172,7 +172,7 @@ contains p % last_sqrtkT = p % sqrtkT ! Get distributed offset - if (size(c % material) > 1 .or. size(c % sqrtkT) > 1) then + if (c % material_size() > 1 .or. size(c % sqrtkT) > 1) then ! Distributed instances of this cell have different ! materials/temperatures. Determine which instance this is for ! assigning the matching material/temperature. @@ -204,7 +204,7 @@ contains end if ! Save the material - if (size(c % material) > 1) then + if (c % material_size() > 1) then p % material = c % material(offset + 1) else p % material = c % material(1) diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 3f36cf299a..8399817fa4 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -8,6 +8,7 @@ #include "constants.h" #include "error.h" #include "lattice.h" +#include "material.h" namespace openmc { @@ -15,7 +16,7 @@ namespace openmc { //============================================================================== void -adjust_indices_c() +adjust_indices() { // Adjust material/fill idices. for (Cell *c : global_cells) { @@ -36,8 +37,21 @@ adjust_indices_c() fatal_error(err_msg); } } else { - //TODO: materials c->type = FILL_MATERIAL; + for (auto it = c->material.begin(); it != c->material.end(); it++) { + int32_t mid = *it; + if (mid != MATERIAL_VOID) { + auto search = material_map.find(mid); + if (search != material_map.end()) { + *it = search->second; + } else { + std::stringstream err_msg; + err_msg << "Could not find material " << mid + << " specified on cell " << c->id; + fatal_error(err_msg); + } + } + } } } diff --git a/src/geometry_aux.h b/src/geometry_aux.h index b998ad32c9..fb5940d1c9 100644 --- a/src/geometry_aux.h +++ b/src/geometry_aux.h @@ -13,7 +13,7 @@ namespace openmc { //! Replace Universe, Lattice, and Material IDs with indices. //============================================================================== -extern "C" void adjust_indices_c(); +extern "C" void adjust_indices(); //============================================================================== //! Figure out which Universe is the root universe. diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 5a2dca282f..93bc9918b2 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -40,12 +40,6 @@ module geometry_header integer(C_INT) :: type end function cell_type_c - subroutine cell_set_type_c(cell_ptr, type) bind(C, name='cell_set_type') - import C_PTR, C_INT - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT), intent(in), value :: type - end subroutine cell_set_type_c - function cell_universe_c(cell_ptr) bind(C, name='cell_universe') & result(universe) import C_PTR, C_INT32_T @@ -53,25 +47,12 @@ module geometry_header integer(C_INT32_T) :: universe end function cell_universe_c - subroutine cell_set_universe_c(cell_ptr, universe) & - bind(C, name='cell_set_universe') - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT32_T), intent(in), value :: universe - end subroutine cell_set_universe_c - function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill) import C_PTR, C_INT32_T type(C_PTR), intent(in), value :: cell_ptr integer(C_INT32_T) :: fill end function cell_fill_c - function cell_fill_ptr(cell_ptr) bind(C) result(fill_ptr) - import C_PTR - type(C_PTR), intent(in), value :: cell_ptr - type(C_PTR) :: fill_ptr - end function cell_fill_ptr - function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') & result(n_instances) import C_PTR, C_INT32_T @@ -79,6 +60,21 @@ module geometry_header integer(C_INT32_T) :: n_instances end function cell_n_instances_c + function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') & + result(n) + import C_PTR, C_INT + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT) :: n + end function cell_material_size_c + + function cell_material_c(cell_ptr, i) bind(C, name='cell_material') & + result(mat) + import C_PTR, C_INT, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT), intent(in), value :: i + integer(C_INT32_T) :: mat + end function cell_material_c + function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple) import C_PTR, C_BOOL type(C_PTR), intent(in), value :: cell_ptr @@ -253,9 +249,6 @@ module geometry_header type Cell type(C_PTR) :: ptr - integer, allocatable :: material(:) ! Material within cell. Multiple - ! materials for distribcell - ! instances. 0 signifies a universe integer, allocatable :: region(:) ! Definition of spatial region as ! Boolean expression of half-spaces integer :: distribcell_index ! Index corresponding to this cell in @@ -274,11 +267,11 @@ module geometry_header procedure :: id => cell_id procedure :: set_id => cell_set_id procedure :: type => cell_type - procedure :: set_type => cell_set_type procedure :: universe => cell_universe - procedure :: set_universe => cell_set_universe procedure :: fill => cell_fill procedure :: n_instances => cell_n_instances + procedure :: material_size => cell_material_size + procedure :: material => cell_material procedure :: simple => cell_simple procedure :: distance => cell_distance procedure :: offset => cell_offset @@ -389,24 +382,12 @@ contains type = cell_type_c(this % ptr) end function cell_type - subroutine cell_set_type(this, type) - class(Cell), intent(in) :: this - integer(C_INT), intent(in) :: type - call cell_set_type_c(this % ptr, type) - end subroutine cell_set_type - function cell_universe(this) result(universe) class(Cell), intent(in) :: this integer(C_INT32_T) :: universe universe = cell_universe_c(this % ptr) end function cell_universe - subroutine cell_set_universe(this, universe) - class(Cell), intent(in) :: this - integer(C_INT32_T), intent(in) :: universe - call cell_set_universe_c(this % ptr, universe) - end subroutine cell_set_universe - function cell_fill(this) result(fill) class(Cell), intent(in) :: this integer(C_INT32_T) :: fill @@ -419,6 +400,19 @@ contains n_instances = cell_n_instances_c(this % ptr) end function cell_n_instances + function cell_material_size(this) result(n) + class(Cell), intent(in) :: this + integer(C_INT) :: n + n = cell_material_size_c(this % ptr) + end function cell_material_size + + function cell_material(this, i) result(mat) + class(Cell), intent(in) :: this + integer, intent(in) :: i + integer(C_INT32_T) :: mat + mat = cell_material_c(this % ptr, i) + end function cell_material + function cell_simple(this) result(simple) class(Cell), intent(in) :: this logical(C_BOOL) :: simple @@ -468,7 +462,7 @@ contains if (present(sab_temps)) allocate(sab_temps(n_sab_tables)) do i = 1, size(cells) - do j = 1, size(cells(i) % material) + do j = 1, cells(i) % material_size() ! Skip any non-material cells and void materials if (cells(i) % material(j) == NONE .or. & cells(i) % material(j) == MATERIAL_VOID) cycle @@ -598,33 +592,6 @@ contains end function openmc_get_cell_index - function openmc_cell_get_fill(index, type, indices, n) result(err) bind(C) - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT), intent(out) :: type - integer(C_INT32_T), intent(out) :: n - type(C_PTR), intent(out) :: indices - integer(C_INT) :: err - - err = 0 - if (index >= 1 .and. index <= size(cells)) then - associate (c => cells(index)) - type = c % type() - select case (type) - case (FILL_MATERIAL) - n = size(c % material) - indices = C_LOC(c % material(1)) - case (FILL_UNIVERSE, FILL_LATTICE) - n = 1 - indices = cell_fill_ptr(c % ptr) - end select - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in cells array is out of bounds.") - end if - end function openmc_cell_get_fill - - function openmc_cell_get_id(index, id) result(err) bind(C) ! Return the ID of a cell integer(C_INT32_T), value :: index @@ -641,49 +608,6 @@ contains end function openmc_cell_get_id - function openmc_cell_set_fill(index, type, n, indices) result(err) bind(C) - ! Set the fill for a cell - integer(C_INT32_T), value, intent(in) :: index ! index in cells - integer(C_INT), value, intent(in) :: type - integer(c_INT32_T), value, intent(in) :: n - integer(C_INT32_T), intent(in) :: indices(n) - integer(C_INT) :: err - - integer :: i, j - - err = 0 - if (index >= 1 .and. index <= size(cells)) then - associate (c => cells(index)) - select case (type) - case (FILL_MATERIAL) - if (allocated(c % material)) deallocate(c % material) - allocate(c % material(n)) - - call c % set_type(FILL_MATERIAL) - do i = 1, n - j = indices(i) - if ((j >= 1 .and. j <= n_materials) .or. j == MATERIAL_VOID) then - c % material(i) = j - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index " // trim(to_str(j)) // " in the & - &materials array is out of bounds.") - end if - end do - case (FILL_UNIVERSE) - call c % set_type(FILL_UNIVERSE) - case (FILL_LATTICE) - call c % set_type(FILL_LATTICE) - end select - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in cells array is out of bounds.") - end if - - end function openmc_cell_set_fill - - function openmc_cell_set_id(index, id) result(err) bind(C) ! Set the ID of a cell integer(C_INT32_T), value, intent(in) :: index diff --git a/src/input_xml.F90 b/src/input_xml.F90 index aef9a14c07..c10d31deab 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -48,8 +48,8 @@ module input_xml save interface - subroutine adjust_indices_c() bind(C) - end subroutine adjust_indices_c + subroutine adjust_indices() bind(C) + end subroutine adjust_indices subroutine allocate_offset_tables(n_maps) bind(C) import C_INT @@ -1114,42 +1114,6 @@ contains // to_str(c % id())) end if - ! Read material - if (check_for_node(node_cell, "material")) then - n_mats = node_word_count(node_cell, "material") - - if (n_mats > 0) then - allocate(sarray(n_mats)) - call get_node_array(node_cell, "material", sarray) - - allocate(c % material(n_mats)) - do j = 1, n_mats - select case(trim(to_lower(sarray(j)))) - case ('void') - c % material(j) = MATERIAL_VOID - case default - c % material(j) = int(str_to_int(sarray(j)), 4) - - ! Check for error - if (c % material(j) == ERROR_INT) then - call fatal_error("Invalid material specified on cell " & - // to_str(c % id())) - end if - end select - end do - - deallocate(sarray) - - else - allocate(c % material(1)) - c % material(1) = NONE - end if - - else - allocate(c % material(1)) - c % material(1) = NONE - end if - ! Check for region specification (also under deprecated name surfaces) if (check_for_node(node_cell, "surfaces")) then call warning("The use of 'surfaces' is deprecated and will be & @@ -1550,7 +1514,7 @@ contains ! Get pointer to list of XML call get_node_list(root, "material", node_mat_list) - ! Allocate cells array + ! Allocate materials array n_materials = size(node_mat_list) allocate(materials(n_materials)) allocate(material_temps(n_materials)) @@ -1575,12 +1539,6 @@ contains mat % depletable = .true. end if - ! Check to make sure 'id' hasn't been used - if (material_dict % has(mat % id())) then - call fatal_error("Two or more materials use the same unique ID: " & - // to_str(mat % id())) - end if - ! Copy material name if (check_for_node(node_mat, "name")) then call get_node_value(node_mat, "name", mat % name) @@ -3861,10 +3819,10 @@ contains ! Set the number of temperatures equal to the number of materials. deallocate(cells(i) % sqrtkT) - allocate(cells(i) % sqrtkT(size(cells(i) % material))) + allocate(cells(i) % sqrtkT(cells(i) % material_size())) ! Check each of the cell materials for temperature data. - do j = 1, size(cells(i) % material) + do j = 1, cells(i) % material_size() ! Arbitrarily set void regions to 0K. if (cells(i) % material(j) == MATERIAL_VOID) then cells(i) % sqrtkT(j) = ZERO @@ -3932,45 +3890,6 @@ contains end subroutine read_multipole_data -!=============================================================================== -! ADJUST_INDICES changes the values for 'surfaces' for each cell and the -! material index assigned to each to the indices in the surfaces and material -! array rather than the unique IDs assigned to each surface and material. Also -! assigns boundary conditions to surfaces based on those read into the bc_dict -! dictionary -!=============================================================================== - - subroutine adjust_indices() - - integer :: i ! index for various purposes - integer :: j ! index for various purposes - integer :: id ! user-specified id - - call adjust_indices_c() - - do i = 1, n_cells - associate (c => cells(i)) - - ! ======================================================================= - ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - - if (c % material(1) /= NONE) then - do j = 1, size(c % material) - id = c % material(j) - if (id == MATERIAL_VOID) then - else if (material_dict % has(id)) then - c % material(j) = material_dict % get(id) - else - call fatal_error("Could not find material " // trim(to_str(id)) & - // " specified on cell " // trim(to_str(c % id()))) - end if - end do - end if - end associate - end do - - end subroutine adjust_indices - !=============================================================================== ! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the ! offsets for distribcells @@ -3994,7 +3913,7 @@ contains ! Find all cells with multiple (distributed) materials or temperatures. do i = 1, n_cells - if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then + if (cells(i) % material_size() > 1 .or. size(cells(i) % sqrtkT) > 1) then call cell_list % add(i) end if end do @@ -4003,10 +3922,10 @@ contains ! number of respective cell instances. do i = 1, n_cells associate (c => cells(i)) - if (size(c % material) > 1) then - if (size(c % material) /= c % n_instances()) then + if (c % material_size() > 1) then + if (c % material_size() /= c % n_instances()) then call fatal_error("Cell " // trim(to_str(c % id())) // " was & - &specified with " // trim(to_str(size(c % material))) & + &specified with " // trim(to_str(c % material_size())) & // " materials but has " // trim(to_str(c % n_instances())) & // " distributed instances. The number of materials must & &equal one or the number of instances.") diff --git a/src/material.cpp b/src/material.cpp index 49f0e24158..3bc2452d39 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1,5 +1,8 @@ #include "material.h" +#include +#include + #include "error.h" #include "xml_interface.h" @@ -14,6 +17,7 @@ namespace openmc { //============================================================================== std::vector global_materials; +std::unordered_map material_map; //============================================================================== // Material implementation @@ -26,8 +30,6 @@ Material::Material(pugi::xml_node material_node) } else { fatal_error("Must specify id of material in materials XML file."); } - - std::cout << "Reading material with id " << id << std::endl; } //============================================================================== @@ -41,6 +43,19 @@ read_materials(pugi::xml_node* node) for (pugi::xml_node material_node: node->children("material")) { global_materials.push_back(new Material(material_node)); } + + // Populate the material map. + for (int i = 0; i < global_materials.size(); i++) { + int32_t mid = global_materials[i]->id; + auto search = material_map.find(mid); + if (search == material_map.end()) { + material_map[mid] = i; + } else { + std::stringstream err_msg; + err_msg << "Two or more materials use the same unique ID: " << mid; + fatal_error(err_msg); + } + } } //============================================================================== @@ -61,6 +76,13 @@ extern "C" { global_materials.push_back(new Material()); } } + + void free_memory_material_c() + { + for (Material *mat : global_materials) {delete mat;} + global_materials.clear(); + material_map.clear(); + } } } // namespace openmc diff --git a/src/material.h b/src/material.h index f2e3df2fb4..0a9a099a44 100644 --- a/src/material.h +++ b/src/material.h @@ -1,6 +1,7 @@ -#ifndef OPENMC_CELL_H -#define OPENMC_CELL_H +#ifndef OPENMC_MATERIAL_H +#define OPENMC_MATERIAL_H +#include #include #include "pugixml.hpp" @@ -14,6 +15,7 @@ namespace openmc { class Material; extern std::vector global_materials; +extern std::unordered_map material_map; //============================================================================== //! A substance with constituent nuclides and thermal scattering data @@ -30,4 +32,4 @@ public: }; } // namespace openmc -#endif // OPENMC_CELL_H +#endif // OPENMC_MATERIAL_H diff --git a/src/material_header.F90 b/src/material_header.F90 index 5bfbadfcf9..3a68a48694 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -484,6 +484,11 @@ contains !=============================================================================== subroutine free_memory_material() + interface + subroutine free_memory_material_c() bind(C) + end subroutine free_memory_material_c + end interface + call free_memory_material_c() n_materials = 0 if (allocated(materials)) deallocate(materials) call material_dict % clear() diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 1b3ed35113..f6289518d2 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -151,7 +151,7 @@ contains allocate(kTs(size(materials))) do i = 1, size(cells) - do j = 1, size(cells(i) % material) + do j = 1, cells(i) % material_size() ! Skip any non-material cells and void materials if (cells(i) % material(j) == NONE .or. & diff --git a/src/summary.F90 b/src/summary.F90 index 507642f713..13bc426637 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -198,7 +198,7 @@ contains case (FILL_MATERIAL) call write_dataset(cell_group, "fill_type", "material") - if (size(c % material) == 1) then + if (c % material_size() == 1) then if (c % material(1) == MATERIAL_VOID) then call write_dataset(cell_group, "material", MATERIAL_VOID) else @@ -206,8 +206,8 @@ contains materials(c % material(1)) % id()) end if else - allocate(cell_materials(size(c % material))) - do j = 1, size(c % material) + allocate(cell_materials(c % material_size())) + do j = 1, c % material_size() if (c % material(j) == MATERIAL_VOID) then cell_materials(j) = MATERIAL_VOID else diff --git a/src/xml_interface.h b/src/xml_interface.h index f10143697a..823948de2f 100644 --- a/src/xml_interface.h +++ b/src/xml_interface.h @@ -17,7 +17,7 @@ check_for_node(pugi::xml_node node, const char *name) } std::string get_node_value(pugi::xml_node node, const char *name, - bool lowercase=false, bool strip=false); + bool lowercase=true, bool strip=true); template std::vector get_node_array(pugi::xml_node node, const char* name) From 20eeb75c90b3dc74c0a309a35b873548e9d9b511 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 15 Aug 2018 21:41:09 -0400 Subject: [PATCH 4/6] Cleanup --- src/cell.cpp | 21 ++++++++-------- src/geometry_aux.cpp | 56 ++++++++++++++++++++--------------------- src/material.cpp | 8 +++++- src/material_header.F90 | 11 +++++--- src/surface.cpp | 52 ++++++++++++++++++-------------------- src/surface.h | 3 ++- 6 files changed, 80 insertions(+), 71 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 52b9a2e44f..e48a3903a2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -204,9 +204,8 @@ Cell::Cell(pugi::xml_node cell_node) fatal_error("Must specify id of cell in geometry XML file."); } - //TODO: don't automatically lowercase cell and surface names if (check_for_node(cell_node, "name")) { - name = get_node_value(cell_node, "name"); + name = get_node_value(cell_node, "name", false); } if (check_for_node(cell_node, "universe")) { @@ -238,8 +237,8 @@ Cell::Cell(pugi::xml_node cell_node) } } else { material.push_back(C_NONE); - material.shrink_to_fit(); } + material.shrink_to_fit(); } else { material.push_back(C_NONE); material.shrink_to_fit(); @@ -318,7 +317,7 @@ Cell::distance(Position r, Direction u, int32_t on_surface) const // Calculate the distance to this surface. // Note the off-by-one indexing bool coincident {token == on_surface}; - double d {surfaces_c[abs(token)-1]->distance(r, u, coincident)}; + double d {global_surfaces[abs(token)-1]->distance(r, u, coincident)}; // Check if this distance is the new minimum. if (d < min_dist) { @@ -359,7 +358,8 @@ Cell::to_hdf5(hid_t cell_group) const region_spec << " |"; } else { // Note the off-by-one indexing - region_spec << " " << copysign(surfaces_c[abs(token)-1]->id, token); + region_spec << " " + << copysign(global_surfaces[abs(token)-1]->id, token); } } write_string(cell_group, "region", region_spec.str(), false); @@ -382,7 +382,7 @@ Cell::contains_simple(Position r, Direction u, int32_t on_surface) const return false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(r, u); + bool sense = global_surfaces[abs(token)-1]->sense(r, u); if (sense != (token > 0)) {return false;} } } @@ -424,7 +424,7 @@ Cell::contains_complex(Position r, Direction u, int32_t on_surface) const stack[i_stack] = false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(r, u);; + bool sense = global_surfaces[abs(token)-1]->sense(r, u); stack[i_stack] = (sense == (token > 0)); } } @@ -454,13 +454,12 @@ read_cells(pugi::xml_node* node) fatal_error("No cells found in geometry.xml!"); } - // Allocate the vector of Cells. - global_cells.reserve(n_cells); - // Loop over XML cell elements and populate the array. + global_cells.reserve(n_cells); for (pugi::xml_node cell_node: node->children("cell")) { global_cells.push_back(new Cell(cell_node)); } + global_cells.shrink_to_fit(); // Populate the Universe vector and map. for (int i = 0; i < global_cells.size(); i++) { @@ -475,6 +474,7 @@ read_cells(pugi::xml_node* node) global_universes[it->second]->cells.push_back(i); } } + global_universes.shrink_to_fit(); } //============================================================================== @@ -524,6 +524,7 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n, return OPENMC_E_OUT_OF_BOUNDS; } } + c.material.shrink_to_fit(); } else if (type == FILL_UNIVERSE) { c.type = FILL_UNIVERSE; } else { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 8399817fa4..7167f24a85 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -19,7 +19,7 @@ void adjust_indices() { // Adjust material/fill idices. - for (Cell *c : global_cells) { + for (Cell* c : global_cells) { if (c->material[0] == C_NONE) { int32_t id = c->fill; auto search_univ = universe_map.find(id); @@ -56,7 +56,7 @@ adjust_indices() } // Change cell.universe values from IDs to indices. - for (Cell *c : global_cells) { + for (Cell* c : global_cells) { auto search = universe_map.find(c->universe); if (search != universe_map.end()) { //TODO: Remove this off-by-one indexing. @@ -70,7 +70,7 @@ adjust_indices() } // Change all lattice universe values from IDs to indices. - for (Lattice *l : lattices_c) { + for (Lattice* l : lattices_c) { l->adjust_indices(); } } @@ -82,12 +82,12 @@ find_root_universe() { // Find all the universes listed as a cell fill. std::unordered_set fill_univ_ids; - for (Cell *c : global_cells) { + for (Cell* c : global_cells) { fill_univ_ids.insert(c->fill); } // Find all the universes contained in a lattice. - for (Lattice *lat : lattices_c) { + for (Lattice* lat : lattices_c) { for (auto it = lat->begin(); it != lat->end(); ++it) { fill_univ_ids.insert(*it); } @@ -123,13 +123,13 @@ find_root_universe() void allocate_offset_tables(int n_maps) { - for (Cell *c : global_cells) { + for (Cell* c : global_cells) { if (c->type != FILL_MATERIAL) { c->offset.resize(n_maps, C_NONE); } } - for (Lattice *lat : lattices_c) { + for (Lattice* lat : lattices_c) { lat->allocate_offset_table(n_maps); } } @@ -140,7 +140,7 @@ void count_cell_instances(int32_t univ_indx) { for (int32_t cell_indx : global_universes[univ_indx]->cells) { - Cell &c = *global_cells[cell_indx]; + Cell& c = *global_cells[cell_indx]; ++c.n_instances; if (c.type == FILL_UNIVERSE) { @@ -149,7 +149,7 @@ count_cell_instances(int32_t univ_indx) } else if (c.type == FILL_LATTICE) { // This cell contains a lattice. Recurse into the lattice universes. - Lattice &lat = *lattices_c[c.fill]; + Lattice& lat = *lattices_c[c.fill]; for (auto it = lat.begin(); it != lat.end(); ++it) { count_cell_instances(*it); } @@ -169,14 +169,14 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id) int count {0}; for (int32_t cell_indx : global_universes[search_univ]->cells) { - Cell &c = *global_cells[cell_indx]; + Cell& c = *global_cells[cell_indx]; if (c.type == FILL_UNIVERSE) { int32_t next_univ = c.fill; count += count_universe_instances(next_univ, target_univ_id); } else if (c.type == FILL_LATTICE) { - Lattice &lat = *lattices_c[c.fill]; + Lattice& lat = *lattices_c[c.fill]; for (auto it = lat.begin(); it != lat.end(); ++it) { int32_t next_univ = *it; count += count_universe_instances(next_univ, target_univ_id); @@ -192,10 +192,10 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id) void fill_offset_tables(int32_t target_univ_id, int map) { - for (Universe *univ : global_universes) { + for (Universe* univ : global_universes) { int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation. for (int32_t cell_indx : univ->cells) { - Cell &c = *global_cells[cell_indx]; + Cell& c = *global_cells[cell_indx]; if (c.type == FILL_UNIVERSE) { c.offset[map] = offset; @@ -203,7 +203,7 @@ fill_offset_tables(int32_t target_univ_id, int map) offset += count_universe_instances(search_univ, target_univ_id); } else if (c.type == FILL_LATTICE) { - Lattice &lat = *lattices_c[c.fill]; + Lattice& lat = *lattices_c[c.fill]; offset = lat.fill_offset_table(offset, target_univ_id, map); } } @@ -214,7 +214,7 @@ fill_offset_tables(int32_t target_univ_id, int map) std::string distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, - const Universe &search_univ, int32_t offset) + const Universe& search_univ, int32_t offset) { std::stringstream path; @@ -224,7 +224,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, // write to the path and return. for (int32_t cell_indx : search_univ.cells) { if ((cell_indx == target_cell) && (offset == target_offset)) { - Cell &c = *global_cells[cell_indx]; + Cell& c = *global_cells[cell_indx]; path << "c" << c.id; return path.str(); } @@ -236,7 +236,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, std::vector::const_reverse_iterator cell_it {search_univ.cells.crbegin()}; for (; cell_it != search_univ.cells.crend(); ++cell_it) { - Cell &c = *global_cells[*cell_it]; + Cell& c = *global_cells[*cell_it]; // Material cells don't contain other cells so ignore them. if (c.type != FILL_MATERIAL) { @@ -244,7 +244,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, if (c.type == FILL_UNIVERSE) { temp_offset = offset + c.offset[map]; } else { - Lattice &lat = *lattices_c[c.fill]; + Lattice& lat = *lattices_c[c.fill]; int32_t indx = lat.universes.size()*map + lat.begin().indx; temp_offset = offset + lat.offsets[indx]; } @@ -256,7 +256,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, } // Add the cell to the path string. - Cell &c = *global_cells[*cell_it]; + Cell& c = *global_cells[*cell_it]; path << "c" << c.id << "->"; if (c.type == FILL_UNIVERSE) { @@ -267,7 +267,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, return path.str(); } else { // Recurse into the lattice cell. - Lattice &lat = *lattices_c[c.fill]; + Lattice& lat = *lattices_c[c.fill]; path << "l" << lat.id; for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) { int32_t indx = lat.universes.size()*map + it.indx; @@ -289,7 +289,7 @@ int distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, int32_t root_univ) { - Universe &root = *global_universes[root_univ]; + Universe& root = *global_universes[root_univ]; std::string path_ {distribcell_path_inner(target_cell, map, target_offset, root, 0)}; return path_.size() + 1; @@ -299,9 +299,9 @@ distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, void distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ, char *path) + int32_t root_univ, char* path) { - Universe &root = *global_universes[root_univ]; + Universe& root = *global_universes[root_univ]; std::string path_ {distribcell_path_inner(target_cell, map, target_offset, root, 0)}; path_.copy(path, path_.size()); @@ -316,12 +316,12 @@ maximum_levels(int32_t univ) int levels_below {0}; for (int32_t cell_indx : global_universes[univ]->cells) { - Cell &c = *global_cells[cell_indx]; + Cell& c = *global_cells[cell_indx]; if (c.type == FILL_UNIVERSE) { int32_t next_univ = c.fill; levels_below = std::max(levels_below, maximum_levels(next_univ)); } else if (c.type == FILL_LATTICE) { - Lattice &lat = *lattices_c[c.fill]; + Lattice& lat = *lattices_c[c.fill]; for (auto it = lat.begin(); it != lat.end(); ++it) { int32_t next_univ = *it; levels_below = std::max(levels_below, maximum_levels(next_univ)); @@ -338,16 +338,16 @@ maximum_levels(int32_t univ) void free_memory_geometry_c() { - for (Cell *c : global_cells) {delete c;} + for (Cell* c : global_cells) {delete c;} global_cells.clear(); cell_map.clear(); n_cells = 0; - for (Universe *u : global_universes) {delete u;} + for (Universe* u : global_universes) {delete u;} global_universes.clear(); universe_map.clear(); - for (Lattice *lat : lattices_c) {delete lat;} + for (Lattice* lat : lattices_c) {delete lat;} lattices_c.clear(); lattice_map.clear(); } diff --git a/src/material.cpp b/src/material.cpp index 3bc2452d39..7979fc5f09 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -43,6 +43,7 @@ read_materials(pugi::xml_node* node) for (pugi::xml_node material_node: node->children("material")) { global_materials.push_back(new Material(material_node)); } + global_materials.shrink_to_fit(); // Populate the material map. for (int i = 0; i < global_materials.size(); i++) { @@ -67,7 +68,12 @@ extern "C" { int32_t material_id(Material* mat) {return mat->id;} - void material_set_id(Material* mat, int32_t id) {mat->id = id;} + void material_set_id(Material* mat, int32_t id, int32_t index) + { + mat->id = id; + //TODO: off-by-one + material_map[id] = index - 1; + } void extend_materials_c(int32_t n) { diff --git a/src/material_header.F90 b/src/material_header.F90 index 3a68a48694..a9de08ec4c 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -42,10 +42,12 @@ module material_header integer(C_INT32_T) :: id end function material_id_c - subroutine material_set_id_c(mat_ptr, id) bind(C, name='material_set_id') + subroutine material_set_id_c(mat_ptr, id, index) & + bind(C, name='material_set_id') import C_PTR, C_INT32_T type(C_PTR), intent(in), value :: mat_ptr integer(C_INT32_T), intent(in), value :: id + integer(C_INT32_T), intent(in), value :: index end subroutine material_set_id_c subroutine extend_materials_c(n) bind(C) @@ -122,10 +124,11 @@ contains id = material_id_c(this % ptr) end function material_id - subroutine material_set_id(this, id) + subroutine material_set_id(this, id, index) class(Material), intent(in) :: this integer(C_INT32_T), intent(in) :: id - call material_set_id_c(this % ptr, id) + integer(C_INT32_T), intent(in) :: index + call material_set_id_c(this % ptr, id, index) end subroutine material_set_id function material_set_density(this, density) result(err) @@ -674,7 +677,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= n_materials) then - call materials(index) % set_id(id) + call materials(index) % set_id(id, index) call material_dict % set(id, index) err = 0 else diff --git a/src/surface.cpp b/src/surface.cpp index 99cfd82b05..b18d7c9195 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -27,7 +27,7 @@ extern "C" const int BC_PERIODIC {3}; int32_t n_surfaces; -Surface** surfaces_c; +std::vector global_surfaces; std::map surface_map; @@ -147,7 +147,7 @@ Surface::Surface(pugi::xml_node surf_node) } if (check_for_node(surf_node, "name")) { - name = get_node_value(surf_node, "name"); + name = get_node_value(surf_node, "name", false); } if (check_for_node(surf_node, "boundary")) { @@ -1031,10 +1031,8 @@ read_surfaces(pugi::xml_node* node) fatal_error("No surfaces found in geometry.xml!"); } - // Allocate the array of Surface pointers. - surfaces_c = new Surface* [n_surfaces]; - // Loop over XML surface elements and populate the array. + global_surfaces.reserve(n_surfaces); { pugi::xml_node surf_node; int i_surf; @@ -1043,40 +1041,40 @@ read_surfaces(pugi::xml_node* node) std::string surf_type = get_node_value(surf_node, "type", true, true); if (surf_type == "x-plane") { - surfaces_c[i_surf] = new SurfaceXPlane(surf_node); + global_surfaces.push_back(new SurfaceXPlane(surf_node)); } else if (surf_type == "y-plane") { - surfaces_c[i_surf] = new SurfaceYPlane(surf_node); + global_surfaces.push_back(new SurfaceYPlane(surf_node)); } else if (surf_type == "z-plane") { - surfaces_c[i_surf] = new SurfaceZPlane(surf_node); + global_surfaces.push_back(new SurfaceZPlane(surf_node)); } else if (surf_type == "plane") { - surfaces_c[i_surf] = new SurfacePlane(surf_node); + global_surfaces.push_back(new SurfacePlane(surf_node)); } else if (surf_type == "x-cylinder") { - surfaces_c[i_surf] = new SurfaceXCylinder(surf_node); + global_surfaces.push_back(new SurfaceXCylinder(surf_node)); } else if (surf_type == "y-cylinder") { - surfaces_c[i_surf] = new SurfaceYCylinder(surf_node); + global_surfaces.push_back(new SurfaceYCylinder(surf_node)); } else if (surf_type == "z-cylinder") { - surfaces_c[i_surf] = new SurfaceZCylinder(surf_node); + global_surfaces.push_back(new SurfaceZCylinder(surf_node)); } else if (surf_type == "sphere") { - surfaces_c[i_surf] = new SurfaceSphere(surf_node); + global_surfaces.push_back(new SurfaceSphere(surf_node)); } else if (surf_type == "x-cone") { - surfaces_c[i_surf] = new SurfaceXCone(surf_node); + global_surfaces.push_back(new SurfaceXCone(surf_node)); } else if (surf_type == "y-cone") { - surfaces_c[i_surf] = new SurfaceYCone(surf_node); + global_surfaces.push_back(new SurfaceYCone(surf_node)); } else if (surf_type == "z-cone") { - surfaces_c[i_surf] = new SurfaceZCone(surf_node); + global_surfaces.push_back(new SurfaceZCone(surf_node)); } else if (surf_type == "quadric") { - surfaces_c[i_surf] = new SurfaceQuadric(surf_node); + global_surfaces.push_back(new SurfaceQuadric(surf_node)); } else { std::stringstream err_msg; @@ -1085,10 +1083,11 @@ read_surfaces(pugi::xml_node* node) } } } + global_surfaces.shrink_to_fit(); // Fill the surface map. for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { - int id = surfaces_c[i_surf]->id; + int id = global_surfaces[i_surf]->id; auto in_map = surface_map.find(id); if (in_map == surface_map.end()) { surface_map[id] = i_surf; @@ -1104,9 +1103,9 @@ read_surfaces(pugi::xml_node* node) zmin {INFTY}, zmax {-INFTY}; int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax; for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { - if (surfaces_c[i_surf]->bc == BC_PERIODIC) { + if (global_surfaces[i_surf]->bc == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface* surf_base = surfaces_c[i_surf]; + Surface* surf_base = global_surfaces[i_surf]; PeriodicSurface* surf = dynamic_cast(surf_base); // Make sure this surface inherits from PeriodicSurface. @@ -1149,9 +1148,9 @@ read_surfaces(pugi::xml_node* node) // Set i_periodic for periodic BC surfaces. for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { - if (surfaces_c[i_surf]->bc == BC_PERIODIC) { + if (global_surfaces[i_surf]->bc == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface* surf_base = surfaces_c[i_surf]; + Surface* surf_base = global_surfaces[i_surf]; PeriodicSurface* surf = dynamic_cast(surf_base); // Also try downcasting to the SurfacePlane type (which must be handled @@ -1198,7 +1197,7 @@ read_surfaces(pugi::xml_node* node) } // Make sure the opposite surface is also periodic. - if (surfaces_c[surf->i_periodic]->bc != BC_PERIODIC) { + if (global_surfaces[surf->i_periodic]->bc != BC_PERIODIC) { std::stringstream err_msg; err_msg << "Could not find matching surface for periodic boundary " "condition on surface " << surf->id; @@ -1213,7 +1212,7 @@ read_surfaces(pugi::xml_node* node) //============================================================================== extern "C" { - Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];} + Surface* surface_pointer(int surf_ind) {return global_surfaces[surf_ind];} int surface_id(Surface* surf) {return surf->id;} @@ -1264,9 +1263,8 @@ extern "C" { void free_memory_surfaces_c() { - for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];} - delete surfaces_c; - surfaces_c = nullptr; + for (Surface* surf : global_surfaces) {delete surf;} + global_surfaces.clear(); n_surfaces = 0; surface_map.clear(); } diff --git a/src/surface.h b/src/surface.h index d43b24ae06..746db5c3b5 100644 --- a/src/surface.h +++ b/src/surface.h @@ -4,6 +4,7 @@ #include #include // For numeric_limits #include +#include #include "hdf5.h" #include "pugixml.hpp" @@ -31,7 +32,7 @@ extern "C" const int BC_PERIODIC; extern "C" int32_t n_surfaces; class Surface; -extern Surface** surfaces_c; +extern std::vector global_surfaces; extern std::map surface_map; From f70cd29dede7eea37d37d590159672255a6783d1 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 16 Aug 2018 14:28:28 -0400 Subject: [PATCH 5/6] Address #1045 comments --- docs/source/io_formats/summary.rst | 2 +- openmc/capi/cell.py | 11 +++--- src/cell.cpp | 54 ++++++++++++++---------------- src/constants.F90 | 2 +- src/constants.h | 3 +- src/geometry_aux.cpp | 2 +- src/geometry_header.F90 | 10 +++--- src/input_xml.F90 | 4 +-- src/lattice.cpp | 26 +++++++------- src/material.cpp | 7 ++-- src/mgxs_data.F90 | 8 +++-- src/surface.cpp | 4 +-- src/xml_interface.h | 7 ++-- 13 files changed, 71 insertions(+), 69 deletions(-) diff --git a/docs/source/io_formats/summary.rst b/docs/source/io_formats/summary.rst index 05c6fff9c0..1f30c69bb1 100644 --- a/docs/source/io_formats/summary.rst +++ b/docs/source/io_formats/summary.rst @@ -34,7 +34,7 @@ The current version of the summary file format is 6.0. 'material', 'universe', or 'lattice'. - **material** (*int* or *int[]*) -- Unique ID of the material(s) assigned to the cell. This dataset is present only if fill_type is - set to 'normal'. The value '-2' signifies void material. The data + set to 'normal'. The value '-1' signifies void material. The data is an array if the cell uses distributed materials, otherwise it is a scalar. - **temperature** (*double[]*) -- Temperature of the cell in Kelvin. diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 75986d2988..4b994ac168 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -107,9 +107,12 @@ class Cell(_FortranObjectWithID): if fill_type.value == 1: if n.value > 1: #TODO: off-by-one - return [Material(index=i+1) for i in indices[:n.value]] + return [Material(index=i+1 if i >= 0 else i) + for i in indices[:n.value]] else: - return Material(index=indices[0]+1) + #TODO: off-by-one + index = indices[0] + 1 if indices[0] >= 0 else indices[0] + return Material(index=index) else: raise NotImplementedError @@ -117,14 +120,14 @@ class Cell(_FortranObjectWithID): def fill(self, fill): if isinstance(fill, Iterable): n = len(fill) - indices = (c_int32*n)(*(m._index if m is not None else -2 + indices = (c_int32*n)(*(m._index if m is not None else -1 for m in fill)) _dll.openmc_cell_set_fill(self._index, 1, n, indices) elif isinstance(fill, Material): indices = (c_int32*1)(fill._index) _dll.openmc_cell_set_fill(self._index, 1, 1, indices) elif fill is None: - indices = (c_int32*1)(-2) + indices = (c_int32*1)(-1) _dll.openmc_cell_set_fill(self._index, 1, 1, indices) def set_temperature(self, T, instance=None): diff --git a/src/cell.cpp b/src/cell.cpp index e48a3903a2..f19930351b 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -199,23 +199,38 @@ generate_rpn(int32_t cell_id, std::vector infix) Cell::Cell(pugi::xml_node cell_node) { if (check_for_node(cell_node, "id")) { - id = stoi(get_node_value(cell_node, "id")); + id = std::stoi(get_node_value(cell_node, "id")); } else { fatal_error("Must specify id of cell in geometry XML file."); } if (check_for_node(cell_node, "name")) { - name = get_node_value(cell_node, "name", false); + name = get_node_value(cell_node, "name"); } if (check_for_node(cell_node, "universe")) { - universe = stoi(get_node_value(cell_node, "universe")); + universe = std::stoi(get_node_value(cell_node, "universe")); } else { universe = 0; } - if (check_for_node(cell_node, "fill")) { - fill = stoi(get_node_value(cell_node, "fill")); + // Make sure that either material or fill was specified, but not both. + bool fill_present = check_for_node(cell_node, "fill"); + bool material_present = check_for_node(cell_node, "material"); + if (!(fill_present || material_present)) { + std::stringstream err_msg; + err_msg << "Neither material nor fill was specified for cell " << id; + fatal_error(err_msg); + } + if (fill_present && material_present) { + std::stringstream err_msg; + err_msg << "Cell " << id << " has both a material and a fill specified; " + << "only one can be specified per cell"; + fatal_error(err_msg); + } + + if (fill_present) { + fill = std::stoi(get_node_value(cell_node, "fill")); } else { fill = C_NONE; } @@ -223,9 +238,9 @@ Cell::Cell(pugi::xml_node cell_node) // Read the material element. There can be zero materials (filled with a // universe), more than one material (distribmats), and some materials may // be "void". - if (check_for_node(cell_node, "material")) { + if (material_present) { std::vector mats - {get_node_array(cell_node, "material")}; + {get_node_array(cell_node, "material", true)}; if (mats.size() > 0) { material.reserve(mats.size()); for (std::string mat : mats) { @@ -235,28 +250,12 @@ Cell::Cell(pugi::xml_node cell_node) material.push_back(std::stoi(mat)); } } + material.shrink_to_fit(); } else { - material.push_back(C_NONE); + std::stringstream err_msg; + err_msg << "An empty material element was specified for cell " << id; + fatal_error(err_msg); } - material.shrink_to_fit(); - } else { - material.push_back(C_NONE); - material.shrink_to_fit(); - } - - // Make sure that either material or fill was specified. - if ((material[0] == C_NONE) && (fill == C_NONE)) { - std::stringstream err_msg; - err_msg << "Neither material nor fill was specified for cell " << id; - fatal_error(err_msg); - } - - // Make sure that material and fill haven't been specified simultaneously. - if ((material[0] != C_NONE) && (fill != C_NONE)) { - std::stringstream err_msg; - err_msg << "Cell " << id << " has both a material and a fill specified; " - << "only one can be specified per cell"; - fatal_error(err_msg); } // Read the region specification. @@ -562,7 +561,6 @@ extern "C" { int32_t cell_material(Cell* c, int i) { int32_t mat = c->material[i-1]; - if (mat == C_NONE) return F90_NONE; if (mat == MATERIAL_VOID) return MATERIAL_VOID; return mat + 1; } diff --git a/src/constants.F90 b/src/constants.F90 index 4b0f392fd7..ab85b277a6 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -129,7 +129,7 @@ module constants integer(C_INT), bind(C, name='FILL_LATTICE') :: FILL_LATTICE_C = FILL_LATTICE ! Void material - integer, parameter :: MATERIAL_VOID = -2 + integer, parameter :: MATERIAL_VOID = -1 ! Flag to say that the outside of a lattice is not defined integer, parameter :: NO_OUTER_UNIVERSE = -1 diff --git a/src/constants.h b/src/constants.h index 72c72bad3f..b6bb82d375 100644 --- a/src/constants.h +++ b/src/constants.h @@ -112,7 +112,7 @@ constexpr char SUBSHELLS[][4] { // Void material // TODO: refactor and remove -constexpr int MATERIAL_VOID {-2}; +constexpr int MATERIAL_VOID {-1}; // ============================================================================ // CROSS SECTION RELATED CONSTANTS @@ -434,7 +434,6 @@ constexpr int DIFF_NUCLIDE_DENSITY {2}; constexpr int DIFF_TEMPERATURE {3}; constexpr int C_NONE {-1}; -constexpr int F90_NONE {0}; // Interpolation rules enum class Interpolation { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 7167f24a85..b202620715 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -20,7 +20,7 @@ adjust_indices() { // Adjust material/fill idices. for (Cell* c : global_cells) { - if (c->material[0] == C_NONE) { + if (c->fill != C_NONE) { int32_t id = c->fill; auto search_univ = universe_map.find(id); auto search_lat = lattice_map.find(id); diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 93bc9918b2..7714601a42 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -4,7 +4,7 @@ module geometry_header use algorithm, only: find use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, & - MATERIAL_VOID, NONE + MATERIAL_VOID use dict_header, only: DictCharInt, DictIntInt use hdf5_interface, only: HID_T use material_header, only: Material, materials, material_dict, n_materials @@ -462,10 +462,12 @@ contains if (present(sab_temps)) allocate(sab_temps(n_sab_tables)) do i = 1, size(cells) + ! Skip non-material cells. + if (cells(i) % fill() /= C_NONE) cycle + do j = 1, cells(i) % material_size() - ! Skip any non-material cells and void materials - if (cells(i) % material(j) == NONE .or. & - cells(i) % material(j) == MATERIAL_VOID) cycle + ! Skip void materials + if (cells(i) % material(j) == MATERIAL_VOID) cycle ! Get temperature of cell (rounding to nearest integer) if (size(cells(i) % sqrtkT) > 1) then diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c10d31deab..23472e21af 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1215,7 +1215,7 @@ contains n = node_word_count(node_cell, "temperature") if (n > 0) then ! Make sure this is a "normal" cell. - if (c % material(1) == NONE) call fatal_error("Cell " & + if (c % fill() /= C_NONE) call fatal_error("Cell " & // trim(to_str(c % id())) // " was specified with a temperature & &but no material. Temperature specification is only valid for & &cells filled with a material.") @@ -3814,7 +3814,7 @@ contains do i = 1, n_cells ! Ignore non-normal cells and cells with defined temperature. - if (cells(i) % material(1) == NONE) cycle + if (cells(i) % fill() /= C_NONE) cycle if (cells(i) % sqrtkT(1) >= ZERO) cycle ! Set the number of temperatures equal to the number of materials. diff --git a/src/lattice.cpp b/src/lattice.cpp index 73dcb98667..e2871d9e1a 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -29,7 +29,7 @@ std::unordered_map lattice_map; Lattice::Lattice(pugi::xml_node lat_node) { if (check_for_node(lat_node, "id")) { - id = stoi(get_node_value(lat_node, "id")); + id = std::stoi(get_node_value(lat_node, "id")); } else { fatal_error("Must specify id of lattice in geometry XML file."); } @@ -39,7 +39,7 @@ Lattice::Lattice(pugi::xml_node lat_node) } if (check_for_node(lat_node, "outer")) { - outer = stoi(get_node_value(lat_node, "outer")); + outer = std::stoi(get_node_value(lat_node, "outer")); } } @@ -141,14 +141,14 @@ RectLattice::RectLattice(pugi::xml_node lat_node) std::string dimension_str {get_node_value(lat_node, "dimension")}; std::vector dimension_words {split(dimension_str)}; if (dimension_words.size() == 2) { - n_cells[0] = stoi(dimension_words[0]); - n_cells[1] = stoi(dimension_words[1]); + n_cells[0] = std::stoi(dimension_words[0]); + n_cells[1] = std::stoi(dimension_words[1]); n_cells[2] = 1; is_3d = false; } else if (dimension_words.size() == 3) { - n_cells[0] = stoi(dimension_words[0]); - n_cells[1] = stoi(dimension_words[1]); - n_cells[2] = stoi(dimension_words[2]); + n_cells[0] = std::stoi(dimension_words[0]); + n_cells[1] = std::stoi(dimension_words[1]); + n_cells[2] = std::stoi(dimension_words[2]); is_3d = true; } else { fatal_error("Rectangular lattice must be two or three dimensions."); @@ -195,7 +195,7 @@ RectLattice::RectLattice(pugi::xml_node lat_node) for (int ix = 0; ix < nx; ix++) { int indx1 = nx*ny*iz + nx*(ny-iy-1) + ix; int indx2 = nx*ny*iz + nx*iy + ix; - universes[indx1] = stoi(univ_words[indx2]); + universes[indx1] = std::stoi(univ_words[indx2]); } } } @@ -400,9 +400,9 @@ HexLattice::HexLattice(pugi::xml_node lat_node) : Lattice {lat_node} { // Read the number of lattice cells in each dimension. - n_rings = stoi(get_node_value(lat_node, "n_rings")); + n_rings = std::stoi(get_node_value(lat_node, "n_rings")); if (check_for_node(lat_node, "n_axial")) { - n_axial = stoi(get_node_value(lat_node, "n_axial")); + n_axial = std::stoi(get_node_value(lat_node, "n_axial")); is_3d = true; } else { n_axial = 1; @@ -476,7 +476,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node) int indx = (2*n_rings-1)*(2*n_rings-1) * m + (2*n_rings-1) * (i_a+n_rings-1) + (i_x+n_rings-1); - universes[indx] = stoi(univ_words[input_index]); + universes[indx] = std::stoi(univ_words[input_index]); input_index++; // Walk the index to the right neighbor (which is not adjacent). i_x += 2; @@ -505,7 +505,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node) int indx = (2*n_rings-1)*(2*n_rings-1) * m + (2*n_rings-1) * (i_a+n_rings-1) + (i_x+n_rings-1); - universes[indx] = stoi(univ_words[input_index]); + universes[indx] = std::stoi(univ_words[input_index]); input_index++; // Walk the index to the right neighbor (which is not adjacent). i_x += 2; @@ -528,7 +528,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node) int indx = (2*n_rings-1)*(2*n_rings-1) * m + (2*n_rings-1) * (i_a+n_rings-1) + (i_x+n_rings-1); - universes[indx] = stoi(univ_words[input_index]); + universes[indx] = std::stoi(univ_words[input_index]); input_index++; // Walk the index to the right neighbor (which is not adjacent). i_x += 2; diff --git a/src/material.cpp b/src/material.cpp index 7979fc5f09..c9ef47f074 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -6,9 +6,6 @@ #include "error.h" #include "xml_interface.h" -//TODO: remove this include -#include - namespace openmc { @@ -26,7 +23,7 @@ std::unordered_map material_map; Material::Material(pugi::xml_node material_node) { if (check_for_node(material_node, "id")) { - id = stoi(get_node_value(material_node, "id")); + id = std::stoi(get_node_value(material_node, "id")); } else { fatal_error("Must specify id of material in materials XML file."); } @@ -40,7 +37,7 @@ extern "C" void read_materials(pugi::xml_node* node) { // Loop over XML material elements and populate the array. - for (pugi::xml_node material_node: node->children("material")) { + for (pugi::xml_node material_node : node->children("material")) { global_materials.push_back(new Material(material_node)); } global_materials.shrink_to_fit(); diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index f6289518d2..b95c3a877a 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -151,11 +151,13 @@ contains allocate(kTs(size(materials))) do i = 1, size(cells) + ! Skip non-material cells + if (cells(i) % fill() /= C_NONE) cycle + do j = 1, cells(i) % material_size() - ! Skip any non-material cells and void materials - if (cells(i) % material(j) == NONE .or. & - cells(i) % material(j) == MATERIAL_VOID) cycle + ! Skip void materials + if (cells(i) % material(j) == MATERIAL_VOID) cycle ! Get temperature of cell (rounding to nearest integer) if (size(cells(i) % sqrtkT) > 1) then diff --git a/src/surface.cpp b/src/surface.cpp index b18d7c9195..77d19673b3 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -141,7 +141,7 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2, Surface::Surface(pugi::xml_node surf_node) { if (check_for_node(surf_node, "id")) { - id = stoi(get_node_value(surf_node, "id")); + id = std::stoi(get_node_value(surf_node, "id")); } else { fatal_error("Must specify id of surface in geometry XML file."); } @@ -247,7 +247,7 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node) : Surface {surf_node} { if (check_for_node(surf_node, "periodic_surface_id")) { - i_periodic = stoi(get_node_value(surf_node, "periodic_surface_id")); + i_periodic = std::stoi(get_node_value(surf_node, "periodic_surface_id")); } } diff --git a/src/xml_interface.h b/src/xml_interface.h index 823948de2f..49addba757 100644 --- a/src/xml_interface.h +++ b/src/xml_interface.h @@ -17,13 +17,14 @@ check_for_node(pugi::xml_node node, const char *name) } std::string get_node_value(pugi::xml_node node, const char *name, - bool lowercase=true, bool strip=true); + bool lowercase=false, bool strip=false); template -std::vector get_node_array(pugi::xml_node node, const char* name) +std::vector get_node_array(pugi::xml_node node, const char* name, + bool lowercase=false) { // Get value of node attribute/child - std::string s {get_node_value(node, name)}; + std::string s {get_node_value(node, name, lowercase)}; // Read values one by one into vector std::stringstream iss {s}; From 6cb9a06f5ad27489af81ee6621cb653d6bb13c8c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 16 Aug 2018 16:44:51 -0400 Subject: [PATCH 6/6] Remove unnecessary shrink_to_fit --- src/cell.cpp | 2 -- src/surface.cpp | 1 - 2 files changed, 3 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index f19930351b..d26fa84147 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -250,7 +250,6 @@ Cell::Cell(pugi::xml_node cell_node) material.push_back(std::stoi(mat)); } } - material.shrink_to_fit(); } else { std::stringstream err_msg; err_msg << "An empty material element was specified for cell " << id; @@ -458,7 +457,6 @@ read_cells(pugi::xml_node* node) for (pugi::xml_node cell_node: node->children("cell")) { global_cells.push_back(new Cell(cell_node)); } - global_cells.shrink_to_fit(); // Populate the Universe vector and map. for (int i = 0; i < global_cells.size(); i++) { diff --git a/src/surface.cpp b/src/surface.cpp index 77d19673b3..421c43d485 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1083,7 +1083,6 @@ read_surfaces(pugi::xml_node* node) } } } - global_surfaces.shrink_to_fit(); // Fill the surface map. for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {