diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 81bdfdc40..11b0e2199 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -1,12 +1,14 @@ -from collections import Mapping +from collections import Mapping, Iterable from ctypes import c_int, c_int32, c_double, c_char_p, POINTER from weakref import WeakValueDictionary import numpy as np +from numpy.ctypeslib import as_array from . import _dll from .core import _View from .error import _error_handler +from .material import MaterialView __all__ = ['CellView', 'cells'] @@ -14,6 +16,12 @@ __all__ = ['CellView', 'cells'] _dll.openmc_cell_get_id.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_cell_get_id.restype = c_int _dll.openmc_cell_get_id.errcheck = _error_handler +_dll.openmc_cell_get_fill.argtypes = [ + c_int32, POINTER(c_int), POINTER(POINTER(c_int32)), POINTER(c_int32)] +_dll.openmc_cell_set_fill.argtypes = [ + c_int32, c_int, c_int32, POINTER(c_int32)] +_dll.openmc_cell_set_fill.restype = c_int +_dll.openmc_cell_set_fill.errcheck = _error_handler _dll.openmc_cell_set_temperature.argtypes = [ c_int32, c_double, POINTER(c_int32)] _dll.openmc_cell_set_temperature.restype = c_int @@ -55,6 +63,34 @@ class CellView(_View): _dll.openmc_cell_get_id(self._index, cell_id) return cell_id.value + @property + def fill(self): + fill_type = c_int() + indices = POINTER(c_int32)() + n = c_int32() + _dll.openmc_cell_get_fill(self._index, fill_type, indices, n) + + if fill_type.value == 1: + if n.value > 1: + return [MaterialView(i) for i in indices[:n.value]] + else: + return MaterialView(indices[0]) + else: + raise NotImplementedError + + @fill.setter + def fill(self, fill): + if isinstance(fill, Iterable): + n = len(fill) + indices = (c_int*n)(*(m._index for m in fill)) + _dll.openmc_cell_set_fill(self._index, 1, 1, indices) + elif isinstance(fill, MaterialView): + materials = [fill] + indices = (c_int*1)(fill._index) + _dll.openmc_cell_set_fill(self._index, 1, 1, indices) + else: + raise NotImplementedError + def set_temperature(self, T, instance=None): """Set the temperature of a cell diff --git a/openmc/capi/error.py b/openmc/capi/error.py index 8e37ea089..7ebf2cb64 100644 --- a/openmc/capi/error.py +++ b/openmc/capi/error.py @@ -62,7 +62,7 @@ def _error_handler(err, func, args): elif err == errcode('e_invalid_type'): raise InvalidTypeError(msg) if err == errcode('e_invalid_id'): - raise IDError(msg) + raise InvalidIDError(msg) elif err == errcode('e_geometry'): raise GeometryError(msg) elif err == errcode('e_data'): diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index dc59560e0..e171e05bd 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -144,11 +144,16 @@ class MaterialFilterView(FilterView): _dll.openmc_material_filter_set_bins(self._index, n, bins) @classmethod - def new(cls, bins=None): + def new(cls, bins=None, filter_id=None): + # Determine filter ID to assign + if filter_id is None: + filter_id = max(filters) + 1 + index = c_int32() _dll.openmc_extend_filters(1, index, None) _dll.openmc_filter_set_type(index, b'material') f = cls(index.value) + f.id = filter_id if bins is not None: f.bins = bins return f diff --git a/openmc/capi/material.py b/openmc/capi/material.py index 81c47f5f9..ba4b65089 100644 --- a/openmc/capi/material.py +++ b/openmc/capi/material.py @@ -13,6 +13,9 @@ from .error import _error_handler __all__ = ['MaterialView', 'materials'] # Material functions +_dll.openmc_extend_materials.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)] +_dll.openmc_extend_materials.restype = c_int +_dll.openmc_extend_materials.errcheck = _error_handler _dll.openmc_get_material_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_material_index.restype = c_int _dll.openmc_get_material_index.errcheck = _error_handler @@ -35,6 +38,9 @@ _dll.openmc_material_set_densities.argtypes = [ c_int32, c_int, POINTER(c_char_p), POINTER(c_double)] _dll.openmc_material_set_densities.restype = c_int _dll.openmc_material_set_densities.errcheck = _error_handler +_dll.openmc_material_set_id.argtypes = [c_int32, c_int32] +_dll.openmc_material_set_id.restype = c_int +_dll.openmc_material_set_id.errcheck = _error_handler class MaterialView(_View): @@ -67,6 +73,10 @@ class MaterialView(_View): _dll.openmc_material_get_id(self._index, mat_id) return mat_id.value + @id.setter + def id(self, mat_id): + _dll.openmc_material_set_id(self._index, mat_id) + @property def nuclides(self): return self._get_densities()[0] @@ -100,7 +110,7 @@ class MaterialView(_View): density_array = as_array(densities, (n.value,)) return nuclide_list, density_array - def add_nuclide(name, density): + def add_nuclide(self, name, density): """Add a nuclide to a material. Parameters @@ -113,6 +123,18 @@ class MaterialView(_View): """ _dll.openmc_material_add_nuclide(self._index, name.encode(), density) + @classmethod + def new(cls, material_id=None): + # Determine ID to assign + if material_id is None: + material_id = max(materials) + 1 + + index = c_int32() + _dll.openmc_extend_materials(1, index, None) + mat = cls(index.value) + mat.id = material_id + return mat + def set_density(self, density): """Set density of a material. diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index db32f98f8..335cbb8e7 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -140,11 +140,17 @@ class TallyView(_View): _dll.openmc_tally_set_scores(self._index, len(scores), scores_) @classmethod - def new(cls): + def new(cls, tally_id=None): + # Determine ID to assign + if tally_id is None: + tally_id = max(tallies) + 1 + index = c_int32() _dll.openmc_extend_tallies(1, index, None) _dll.openmc_tally_set_type(index, b'generic') - return cls(index.value) + tally = cls(index.value) + tally.id = tally_id + return tally class _TallyMapping(Mapping): diff --git a/src/api.F90 b/src/api.F90 index 1c995806b..3b2fdacad 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -34,10 +34,13 @@ module openmc_api private public :: openmc_calculate_volumes public :: openmc_cell_get_id + public :: openmc_cell_get_fill + public :: openmc_cell_set_fill public :: openmc_cell_set_temperature public :: openmc_energy_filter_get_bins public :: openmc_energy_filter_set_bins public :: openmc_extend_filters + public :: openmc_extend_materials public :: openmc_extend_tallies public :: openmc_filter_get_id public :: openmc_filter_get_type @@ -59,6 +62,7 @@ module openmc_api public :: openmc_material_get_densities public :: openmc_material_set_density public :: openmc_material_set_densities + public :: openmc_material_set_id public :: openmc_material_filter_get_bins public :: openmc_material_filter_set_bins public :: openmc_mesh_filter_set_mesh diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 4d77b45e3..b6499f822 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -6,7 +6,7 @@ module geometry_header use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, & MATERIAL_VOID, NONE use dict_header, only: DictCharInt, DictIntInt - use material_header, only: Material, materials, material_dict + use material_header, only: Material, materials, material_dict, n_materials use nuclide_header use sab_header use stl_vector, only: VectorReal @@ -447,6 +447,33 @@ 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 = C_LOC(c % fill) + 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 @@ -463,9 +490,56 @@ 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 fill + 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)) + + c % type = FILL_MATERIAL + do i = 1, n + j = indices(i) + if (j == 0) then + c % material(i) = MATERIAL_VOID + else + if (j >= 1 .and. j <= n_materials) 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 if + end do + case (FILL_UNIVERSE) + c % type = FILL_UNIVERSE + case (FILL_LATTICE) + c % 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_temperature(index, T, instance) result(err) bind(C) ! Set the temperature of a cell - integer(C_INT32_T), value, intent(in) :: index ! cell index in cells + integer(C_INT32_T), value, intent(in) :: index ! index in cells real(C_DOUBLE), value, intent(in) :: T ! temperature integer(C_INT32_T), optional, intent(in) :: instance ! cell instance diff --git a/src/material_header.F90 b/src/material_header.F90 index 1846551a7..eac5ab546 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -14,12 +14,14 @@ module material_header private public :: free_memory_material + public :: openmc_extend_materials public :: openmc_get_material_index public :: openmc_material_add_nuclide public :: openmc_material_get_id public :: openmc_material_get_densities public :: openmc_material_set_density public :: openmc_material_set_densities + public :: openmc_material_set_id !=============================================================================== ! MATERIAL describes a material by its constituent nuclides @@ -28,7 +30,7 @@ module material_header type, public :: Material integer :: id ! unique identifier character(len=104) :: name = "" ! User-defined name - integer :: n_nuclides ! number of nuclides + integer :: n_nuclides = 0 ! number of nuclides integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm @@ -237,6 +239,37 @@ contains ! C API FUNCTIONS !=============================================================================== + function openmc_extend_materials(n, index_start, index_end) result(err) bind(C) + ! Extend the materials array by n elements + integer(C_INT32_T), value, intent(in) :: n + integer(C_INT32_T), optional, intent(out) :: index_start + integer(C_INT32_T), optional, intent(out) :: index_end + integer(C_INT) :: err + + type(Material), allocatable :: temp(:) ! temporary materials array + + if (n_materials == 0) then + ! Allocate materials array + allocate(materials(n)) + else + ! Allocate materials array with increased size + allocate(temp(n_materials + n)) + + ! Move original materials to temporary array + temp(1:n_materials) = materials(:) + + ! Move allocation from temporary array + call move_alloc(FROM=temp, TO=materials) + end if + + ! Return indices in materials array + if (present(index_start)) index_start = n_materials + 1 + if (present(index_end)) index_end = n_materials + n + n_materials = n_materials + n + + err = 0 + end function openmc_extend_materials + function openmc_get_material_index(id, index) result(err) bind(C) ! Returns the index in the materials array of a material with a given ID integer(C_INT32_T), value :: id @@ -269,6 +302,7 @@ contains real(8) :: awr integer, allocatable :: new_nuclide(:) real(8), allocatable :: new_density(:) + logical, allocatable :: new_p0(:) character(:), allocatable :: name_ name_ = to_f_string(name) @@ -277,7 +311,7 @@ contains if (index >= 1 .and. index <= size(materials)) then associate (m => materials(index)) ! Check if nuclide is already in material - do j = 1, size(m % nuclide) + do j = 1, m % n_nuclides k = m % nuclide(j) if (nuclides(k) % name == name_) then awr = nuclides(k) % awr @@ -296,15 +330,20 @@ contains if (err == 0) then ! Extend arrays - n = size(m % nuclide) + n = m % n_nuclides allocate(new_nuclide(n + 1)) - new_nuclide(1:n) = m % nuclide + if (n > 0) new_nuclide(1:n) = m % nuclide call move_alloc(FROM=new_nuclide, TO=m % nuclide) allocate(new_density(n + 1)) - new_density(1:n) = m % atom_density + if (n > 0) new_density(1:n) = m % atom_density call move_alloc(FROM=new_density, TO=m % atom_density) + allocate(new_p0(n + 1)) + if (n > 0) new_p0(1:n) = m % p0 + new_p0(n + 1) = .false. + call move_alloc(FROM=new_p0, TO=m % p0) + ! Append new nuclide/density k = nuclide_dict % get_key(to_lower(name_)) m % nuclide(n + 1) = k @@ -368,6 +407,23 @@ contains end function openmc_material_get_id + function openmc_material_set_id(index, id) result(err) bind(C) + ! Set the ID of a material + integer(C_INT32_T), value, intent(in) :: index + integer(C_INT32_T), value, intent(in) :: id + integer(C_INT) :: err + + if (index >= 1 .and. index <= n_materials) then + materials(index) % id = id + call material_dict % add_key(id, index) + err = 0 + else + err = E_OUT_OF_BOUNDS + call set_errmsg("Index in materials array is out of bounds.") + end if + end function openmc_material_set_id + + function openmc_material_set_density(index, density) result(err) bind(C) ! Set the total density of a material in atom/b-cm integer(C_INT32_T), value, intent(in) :: index @@ -396,14 +452,15 @@ contains integer(C_INT) :: err integer :: i + integer :: stat character(C_CHAR), pointer :: string(:) character(len=:, kind=C_CHAR), allocatable :: name_ if (index >= 1 .and. index <= size(materials)) then associate (m => materials(index)) ! If nuclide/density arrays are not correct size, reallocate - if (n /= size(m % nuclide)) then - deallocate(m % nuclide, m % atom_density, m % p0) + if (n /= m % n_nuclides) then + deallocate(m % nuclide, m % atom_density, m % p0, STAT=stat) allocate(m % nuclide(n), m % atom_density(n), m % p0(n)) end if