Move cell materials array to C++

This commit is contained in:
Sterling Harper 2018-08-15 19:50:05 -04:00
parent e5d154b955
commit 1f11eabc71
17 changed files with 197 additions and 230 deletions

View file

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

View file

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

View file

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

View file

@ -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<std::string> mats
{get_node_array<std::string>(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)

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

@ -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 <material>
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.")

View file

@ -1,5 +1,8 @@
#include "material.h"
#include <string>
#include <sstream>
#include "error.h"
#include "xml_interface.h"
@ -14,6 +17,7 @@ namespace openmc {
//==============================================================================
std::vector<Material*> global_materials;
std::unordered_map<int32_t, int32_t> 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

View file

@ -1,6 +1,7 @@
#ifndef OPENMC_CELL_H
#define OPENMC_CELL_H
#ifndef OPENMC_MATERIAL_H
#define OPENMC_MATERIAL_H
#include <unordered_map>
#include <vector>
#include "pugixml.hpp"
@ -14,6 +15,7 @@ namespace openmc {
class Material;
extern std::vector<Material*> global_materials;
extern std::unordered_map<int32_t, int32_t> 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

View file

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

View file

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

View file

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

View file

@ -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 <typename T>
std::vector<T> get_node_array(pugi::xml_node node, const char* name)