diff --git a/CMakeLists.txt b/CMakeLists.txt index c1228c493..ef9daabd2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -439,6 +439,8 @@ set(LIBOPENMC_CXX_SRC src/cell.h src/error.h src/hdf5_interface.h + src/lattice.cpp + src/lattice.h src/random_lcg.cpp src/random_lcg.h src/surface.cpp diff --git a/src/cell.cpp b/src/cell.cpp index 83942cf73..427e30890 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -301,23 +301,25 @@ Cell::to_hdf5(hid_t cell_group) const //write_int(cell_group, "universe", universe); // Write the region specification. - std::stringstream region_spec {}; - for (int32_t token : region) { - if (token == OP_LEFT_PAREN) { - region_spec << " ("; - } else if (token == OP_RIGHT_PAREN) { - region_spec << " )"; - } else if (token == OP_COMPLEMENT) { - region_spec << " ~"; - } else if (token == OP_INTERSECTION) { - } else if (token == OP_UNION) { - region_spec << " |"; - } else { - // Note the off-by-one indexing - region_spec << " " << copysign(surfaces_c[abs(token)-1]->id, token); + if (!region.empty()) { + std::stringstream region_spec {}; + for (int32_t token : region) { + if (token == OP_LEFT_PAREN) { + region_spec << " ("; + } else if (token == OP_RIGHT_PAREN) { + region_spec << " )"; + } else if (token == OP_COMPLEMENT) { + region_spec << " ~"; + } else if (token == OP_INTERSECTION) { + } else if (token == OP_UNION) { + region_spec << " |"; + } else { + // Note the off-by-one indexing + region_spec << " " << copysign(surfaces_c[abs(token)-1]->id, token); + } } + write_string(cell_group, "region", region_spec.str()); } - write_string(cell_group, "region", region_spec.str()); // close_group(cell_group); } diff --git a/src/cell.h b/src/cell.h index 5f9369e55..59cec1729 100644 --- a/src/cell.h +++ b/src/cell.h @@ -25,6 +25,19 @@ std::vector cells_c; std::map cell_dict; +//============================================================================== +//! A geometry primitive that fills all space and contains cells. +//============================================================================== + +class Universe +{ + public: + int32_t id; //! Unique ID + int32_t type; + std::vector cells; //! Cells within this universe + double x0, y0, z0; //! Translation coordinates. +}; + //============================================================================== //! A geometry primitive that links surfaces, universes, and materials //============================================================================== diff --git a/src/geometry.F90 b/src/geometry.F90 index 74e9731ff..1f241f66d 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -260,7 +260,7 @@ contains ! Particle is outside the lattice. if (lat % outer == NO_OUTER_UNIVERSE) then call p % mark_as_lost("Particle " // trim(to_str(p %id)) & - // " is outside lattice " // trim(to_str(lat % id)) & + // " is outside lattice " // trim(to_str(lat % id())) & // " but the lattice has no defined outer universe.") return else @@ -299,7 +299,7 @@ contains lat => lattices(p % coord(j) % lattice) % obj if (verbosity >= 10 .or. trace) then - call write_message(" Crossing lattice " // trim(to_str(lat % id)) & + call write_message(" Crossing lattice " // trim(to_str(lat % id())) & &// ". Current position (" // trim(to_str(p % coord(j) % lattice_x)) & &// "," // trim(to_str(p % coord(j) % lattice_y)) // "," & &// trim(to_str(p % coord(j) % lattice_z)) // ")") diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 32511ea3d..ad34c735c 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -79,6 +79,29 @@ module geometry_header type(C_PTR), intent(in), value :: cell_ptr 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) + use ISO_C_BINDING + implicit none + integer(C_INT32_T), intent(in), value :: lat_ind + type(C_PTR) :: ptr + end function lattice_pointer_c + + function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT32_T) :: id + end function lattice_id_c + + subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5') + use ISO_C_BINDING + use hdf5 + implicit none + type(C_PTR), intent(in), value :: lat_ptr + integer(HID_T), intent(in), value :: group + end subroutine lattice_to_hdf5_c end interface !=============================================================================== @@ -99,8 +122,9 @@ module geometry_header !=============================================================================== type, abstract :: Lattice - integer :: id ! Universe number for lattice - character(len=104) :: name = "" ! User-defined name + type(C_PTR) :: ptr + + !integer :: id ! Universe number for lattice real(8), allocatable :: pitch(:) ! Pitch along each axis integer, allocatable :: universes(:,:,:) ! Specified universes integer :: outside ! Material to fill area outside @@ -108,6 +132,10 @@ module geometry_header logical :: is_3d ! Lattice has cells on z axis integer, allocatable :: offset(:,:,:,:) ! Distribcell offsets contains + + procedure :: id => lattice_id + procedure :: to_hdf5 => lattice_to_hdf5 + procedure(lattice_are_valid_indices_), deferred :: are_valid_indices procedure(lattice_get_indices_), deferred :: get_indices procedure(lattice_get_local_xyz_), deferred :: get_local_xyz @@ -247,6 +275,18 @@ module geometry_header contains + function lattice_id(this) result(id) + class(Lattice), intent(in) :: this + integer(C_INT32_T) :: id + id = lattice_id_c(this % ptr) + end function lattice_id + + subroutine lattice_to_hdf5(this, group) + class(Lattice), intent(in) :: this + integer(HID_T), intent(in) :: group + call lattice_to_hdf5_c(this % ptr, group) + end subroutine lattice_to_hdf5 + !=============================================================================== function valid_inds_rect(this, i_xyz) result(is_valid) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 71b3b8d64..78dd0bb87 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -60,6 +60,12 @@ module input_xml implicit none type(C_PTR) :: node_ptr end subroutine read_cells + + subroutine read_lattices(node_ptr) bind(C) + use ISO_C_BINDING + implicit none + type(C_PTR) :: node_ptr + end subroutine read_lattices end interface contains @@ -1240,6 +1246,8 @@ contains ! ========================================================================== ! READ LATTICES FROM GEOMETRY.XML + call read_lattices(root % ptr) + ! Get pointer to list of XML call get_node_list(root, "lattice", node_rlat_list) call get_node_list(root, "hex_lattice", node_hlat_list) @@ -1253,30 +1261,13 @@ contains RECT_LATTICES: do i = 1, n_rlats allocate(RectLattice::lattices(i) % obj) lat => lattices(i) % obj + lat % ptr = lattice_pointer_c(i - 1) select type(lat) type is (RectLattice) ! Get pointer to i-th lattice node_lat = node_rlat_list(i) - ! ID of lattice - if (check_for_node(node_lat, "id")) then - call get_node_value(node_lat, "id", lat % id) - else - call fatal_error("Must specify id of lattice in geometry XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (lattice_dict % has(lat % id)) then - call fatal_error("Two or more lattices use the same unique ID: " & - // to_str(lat % id)) - end if - - ! Copy lattice name - if (check_for_node(node_lat, "name")) then - call get_node_value(node_lat, "name", lat % name) - end if - ! Read number of lattice cells in each dimension n = node_word_count(node_lat, "dimension") if (n == 2) then @@ -1341,7 +1332,7 @@ contains n = node_word_count(node_lat, "universes") if (n /= n_x*n_y*n_z) then call fatal_error("Number of universes on does not match & - &size of lattice " // trim(to_str(lat % id)) // ".") + &size of lattice " // trim(to_str(lat % id())) // ".") end if allocate(temp_int_array(n)) @@ -1377,7 +1368,7 @@ contains end if ! Add lattice to dictionary - call lattice_dict % set(lat % id, i) + call lattice_dict % set(lat % id(), i) end select end do RECT_LATTICES @@ -1385,30 +1376,13 @@ 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) select type (lat) type is (HexLattice) ! Get pointer to i-th lattice node_lat = node_hlat_list(i) - ! ID of lattice - if (check_for_node(node_lat, "id")) then - call get_node_value(node_lat, "id", lat % id) - else - call fatal_error("Must specify id of lattice in geometry XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (lattice_dict % has(lat % id)) then - call fatal_error("Two or more lattices use the same unique ID: " & - // to_str(lat % id)) - end if - - ! Copy lattice name - if (check_for_node(node_lat, "name")) then - call get_node_value(node_lat, "name", lat % name) - end if - ! Read number of lattice cells in each dimension call get_node_value(node_lat, "n_rings", lat % n_rings) if (check_for_node(node_lat, "n_axial")) then @@ -1454,7 +1428,7 @@ contains n = node_word_count(node_lat, "universes") if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then call fatal_error("Number of universes on does not match & - &size of lattice " // trim(to_str(lat % id)) // ".") + &size of lattice " // trim(to_str(lat % id())) // ".") end if allocate(temp_int_array(n)) @@ -1564,7 +1538,7 @@ contains end if ! Add lattice to dictionary - call lattice_dict % set(lat % id, n_rlats + i) + call lattice_dict % set(lat % id(), n_rlats + i) end select end do HEX_LATTICES @@ -4353,7 +4327,7 @@ contains else call fatal_error("Invalid universe number " & &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) + &// trim(to_str(lat % id()))) end if end do end do @@ -4374,7 +4348,7 @@ contains else call fatal_error("Invalid universe number " & &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) + &// trim(to_str(lat % id()))) end if end do end do @@ -4388,7 +4362,7 @@ contains else call fatal_error("Invalid universe number " & &// trim(to_str(lat % outer)) & - &// " specified on lattice " // trim(to_str(lat % id))) + &// " specified on lattice " // trim(to_str(lat % id()))) end if end if diff --git a/src/lattice.cpp b/src/lattice.cpp new file mode 100644 index 000000000..73758b110 --- /dev/null +++ b/src/lattice.cpp @@ -0,0 +1,109 @@ +#include "lattice.h" + +#include +#include + +#include "error.h" +#include "hdf5_interface.h" +#include "xml_interface.h" + +//TODO: remove this include +#include + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +// Braces force n_lattices to be defined here, not just declared. +//extern "C" {int32_t n_lattices {0};} + +//Lattice **lattices_c; +std::vector lattices_c; + +std::map lattice_dict; + +//============================================================================== +// Lattice implementation +//============================================================================== + +Lattice::Lattice(pugi::xml_node lat_node) +{ + if (check_for_node(lat_node, "id")) { + id = stoi(get_node_value(lat_node, "id")); + } else { + fatal_error("Must specify id of lattice in geometry XML file."); + } + + if (check_for_node(lat_node, "name")) { + name = get_node_value(lat_node, "name"); + } +} + +void +Lattice::to_hdf5(hid_t lat_group) const +{ + if (!name.empty()) { + write_string(lat_group, "name", name); + } +} + +//============================================================================== +// RectLattice implementation +//============================================================================== + +RectLattice::RectLattice(pugi::xml_node lat_node) + : Lattice {lat_node} +{ +} + +//============================================================================== +// HexLattice implementation +//============================================================================== + +HexLattice::HexLattice(pugi::xml_node lat_node) + : Lattice {lat_node} +{ +} + +//============================================================================== + +extern "C" void +read_lattices(pugi::xml_node *node) +{ + for (pugi::xml_node lat_node: node->children("lattice")) { + lattices_c.push_back(new RectLattice(lat_node)); + } + for (pugi::xml_node lat_node: node->children("hex_lattice")) { + lattices_c.push_back(new HexLattice(lat_node)); + } + + // Fill the lattice dictionary. + for (int i_lat = 0; i_lat < lattices_c.size(); i_lat++) { + int id = lattices_c[i_lat]->id; + auto in_dict = lattice_dict.find(id); + if (in_dict == lattice_dict.end()) { + lattice_dict[id] = i_lat; + } else { + std::stringstream err_msg; + err_msg << "Two or more lattices use the same unique ID: " << id; + fatal_error(err_msg); + } + } +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + Lattice* lattice_pointer(int lat_ind) {return lattices_c[lat_ind];} + + int32_t lattice_id(Lattice *lat) {return lat->id;} + + void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);} +} + +} // namespace openmc diff --git a/src/lattice.h b/src/lattice.h new file mode 100644 index 000000000..d8e1ba39a --- /dev/null +++ b/src/lattice.h @@ -0,0 +1,77 @@ +#ifndef LATTICE_H +#define LATTICE_H + +#include +#include +#include +#include + +#include "hdf5.h" +#include "pugixml/pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +//extern "C" int32_t n_lattice; + +class Lattice; +//extern Lattice **lattices_c; +extern std::vector lattices_c; + +extern std::map lattice_dict; + +//============================================================================== +//! Abstract type for ordered array of universes +//============================================================================== + +class Lattice +{ +public: + int32_t id; //! Universe ID number + std::string name; //! User-defined name + //std::vector pitch; //! Pitch along each basis + //std::vector universes; //! Universes filling each lattice tile + //int32_t outer; //! Universe tiled outside the lattice + //std::vector offset; //! Distribcell offsets + + explicit Lattice(pugi::xml_node lat_node); + + virtual ~Lattice() {} + + //virtual bool are_valid_indices(const int i_xyz[3]) const = 0; + + //virtual void get_indices(const double global_xyz[3], int i_xyz[3]) const = 0; + + //virtual void get_local_xyz(const double global_xyz[3], const int i_xyz[3], + // double local_xyz[3]) const = 0; + + //! Write all information needed to reconstruct the lattice to an HDF5 group. + //! @param group_id An HDF5 group id. + void to_hdf5(hid_t group_id) const; +}; + +//============================================================================== +//============================================================================== + +class RectLattice : public Lattice +{ +public: + explicit RectLattice(pugi::xml_node lat_node); + + virtual ~RectLattice() {} +}; + +class HexLattice : public Lattice +{ +public: + explicit HexLattice(pugi::xml_node lat_node); + + virtual ~HexLattice() {} +}; + +} // namespace openmc +#endif // LATTICE_H diff --git a/src/output.F90 b/src/output.F90 index 483d171ca..81f6e1292 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -246,7 +246,7 @@ contains ! Print information on lattice if (p % coord(i) % lattice /= NONE) then l => lattices(p % coord(i) % lattice) % obj - write(ou,*) ' Lattice = ' // trim(to_str(l % id)) + write(ou,*) ' Lattice = ' // trim(to_str(l % id())) write(ou,*) ' Lattice position = (' // trim(to_str(& p % coord(i) % lattice_x)) // ',' // trim(to_str(& p % coord(i) % lattice_y)) // ')' diff --git a/src/summary.F90 b/src/summary.F90 index 5e5ab778a..59bb90447 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -200,7 +200,7 @@ contains case (FILL_LATTICE) call write_dataset(cell_group, "fill_type", "lattice") - call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id) + call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id()) end select call close_group(cell_group) @@ -258,10 +258,11 @@ contains ! Write information on each lattice LATTICE_LOOP: do i = 1, n_lattices lat => lattices(i)%obj - lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id))) + lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id()))) + + call lat % to_hdf5(lattice_group) ! Write name, pitch, and outer universe - call write_dataset(lattice_group, "name", lat%name) call write_dataset(lattice_group, "pitch", lat%pitch) if (lat % outer > 0) then call write_dataset(lattice_group, "outer", universes(lat % outer) % id) diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index 310af1cb6..e7a9645ae 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -270,7 +270,7 @@ contains type is (RectLattice) ! Write to the geometry stack - path = trim(path) // "->l" // to_str(lat%id) + path = trim(path) // "->l" // to_str(lat%id()) n_x = lat % n_cells(1) n_y = lat % n_cells(2) @@ -332,7 +332,7 @@ contains type is (HexLattice) ! Write to the geometry stack - path = trim(path) // "->l" // to_str(lat%id) + path = trim(path) // "->l" // to_str(lat%id()) n_z = lat % n_axial n_y = 2 * lat % n_rings - 1