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

Partial support for cells defined with union and complement operators.

Browse source 
The Python API still has no concept of operators. The Summary class has a hack
to at least be able to read a summary file for now, basically skipping any
operators that are present.
This commit is contained in:
Paul Romano 2015-09-24 13:55:11 +07:00
parent da0d13e18c
commit acf2e51272
9 changed files with 286 additions and 132 deletions
Download patch file Download diff file Expand all files Collapse all files

17
openmc/summary.py
View file

@ -213,7 +213,7 @@ class Summary(object):
fill = self._f['geometry/cells'][key]['lattice'].value
if 'surfaces' in self._f['geometry/cells'][key].keys():
surfaces = self._f['geometry/cells'][key]['surfaces'][...]
surfaces = self._f['geometry/cells'][key]['surfaces'].value.decode()
else:
surfaces = []
@ -241,12 +241,15 @@ class Summary(object):
self._cell_fills[index] = (fill_type, fill)
# Iterate over all Surfaces and add them to the Cell
for surface_halfspace in surfaces:
halfspace = np.sign(surface_halfspace)
surface_id = abs(surface_halfspace)
surface = self.get_surface_by_id(surface_id)
cell.add_surface(surface, halfspace)
for token in surfaces.split():
try:
surface_halfspace = int(token)
halfspace = np.sign(surface_halfspace)
surface_id = abs(surface_halfspace)
surface = self.get_surface_by_id(surface_id)
cell.add_surface(surface, halfspace)
except ValueError:
pass
# Add the Cell to the global dictionary of all Cells
self.cells[index] = cell

109
src/geometry.F90
View file

@ -17,56 +17,77 @@ module geometry
contains
!===============================================================================
! SIMPLE_CELL_CONTAINS determines whether a given the current coordinates of the
! particle are inside a cell defined as the intersection of a series of surfaces
! CELL_CONTAINS determines if a cell contains the particle at a given
! location. The bounds of the cell are detemined by a logical expression
! involving surface half-spaces. At initialization, the expression was converted
! to RPN notation. In cell_contains, we evaluate the RPN expression using a
! stack, similar to how a RPN calculator would work.
!===============================================================================
function simple_cell_contains(c, p) result(in_cell)
pure function cell_contains(c, p) result(in_cell)
type(Cell), intent(in) :: c
type(Particle), intent(inout) :: p
type(Particle), intent(in) :: p
logical :: in_cell
integer :: i ! index of surfaces in cell
integer :: i_surface ! index in surfaces array (with sign)
logical :: specified_sense ! specified sense of surface in list
integer :: i
integer :: token
logical :: b1, b2
integer :: i_stack
logical :: actual_sense ! sense of particle wrt surface
class(Surface), pointer :: s
logical :: stack(size(c%rpn))
SURFACE_LOOP: do i = 1, c % n_surfaces
! Lookup surface
i_surface = c % surfaces(i)
! Check if the particle is currently on the specified surface
if (i_surface == p % surface) then
! Particle is heading into the cell
cycle
elseif (i_surface == -p % surface) then
! Particle is heading out of the cell
in_cell = .false.
return
i_stack = 0
do i = 1, size(c%rpn)
token = c%rpn(i)
if (token < OP_UNION) then
! If the token is not an operator, evaluate the sense of particle with
! respect to the surface and see if the token matches the sense. If the
! particle's surface attribute is set and matches the token, that
! overrides the determination based on sense().
i_stack = i_stack + 1
if (token == p%surface) then
stack(i_stack) = .true.
elseif (-token == p%surface) then
stack(i_stack) = .false.
else
actual_sense = surfaces(abs(token))%obj%sense(&
p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
stack(i_stack) = (actual_sense .eqv. (token > 0))
end if
else
! If the token is a binary operator (intersection/union), apply it to
! the last two items on the stack. If the token is a unary operator
! (complement), apply it to the last item on the stack.
b1 = stack(i_stack)
select case (token)
case (OP_UNION)
b2 = stack(i_stack - 1)
stack(i_stack - 1) = b1 .or. b2
i_stack = i_stack - 1
case (OP_INTERSECTION)
b2 = stack(i_stack - 1)
stack(i_stack - 1) = b1 .and. b2
i_stack = i_stack - 1
case (OP_COMPLEMENT)
stack(i_stack) = .not. b1
end select
end if
end do
! Determine the specified sense of the surface in the cell and the actual
! sense of the particle with respect to the surface
s => surfaces(abs(i_surface))%obj
actual_sense = s%sense(p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
specified_sense = (c % surfaces(i) > 0)
! Compare sense of point to specified sense
if (actual_sense .neqv. specified_sense) then
in_cell = .false.
return
end if
end do SURFACE_LOOP
! If we've reached here, then the sense matched on every surface or there
! are no surfaces.
in_cell = .true.
end function simple_cell_contains
if (i_stack == 1) then
! The one remaining logical on the stack indicates whether the particle is
! in the cell.
in_cell = stack(i_stack)
else
! This case occurs if there is no surface specification since i_stack will
! still be zero.
in_cell = .true.
end if
end function cell_contains
!===============================================================================
! CHECK_CELL_OVERLAP checks for overlapping cells at the current particle's
! position using simple_cell_contains and the LocalCoord's built up by find_cell
! position using cell_contains and the LocalCoord's built up by find_cell
!===============================================================================
subroutine check_cell_overlap(p)
@ -93,7 +114,7 @@ contains
index_cell = univ % cells(i)
c => cells(index_cell)
if (simple_cell_contains(c, p)) then
if (cell_contains(c, p)) then
! the particle should only be contained in one cell per level
if (index_cell /= p % coord(j) % cell) then
call fatal_error("Overlapping cells detected: " &
@ -162,7 +183,7 @@ contains
c => cells(index_cell)
! Move on to the next cell if the particle is not inside this cell
if (.not. simple_cell_contains(c, p)) cycle
if (.not. cell_contains(c, p)) cycle
! Set cell on this level
p % coord(j) % cell = index_cell
@ -556,7 +577,7 @@ contains
! =======================================================================
! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL
SURFACE_LOOP: do i = 1, cl % n_surfaces
SURFACE_LOOP: do i = 1, size(cl%surfaces)
! check for operators
index_surf = cl%surfaces(i)
coincident = (index_surf == p % surface)
@ -823,10 +844,11 @@ contains
c => cells(i)
! loop over each surface specification
do j = 1, c % n_surfaces
do j = 1, size(c%surfaces)
i_surface = c % surfaces(j)
positive = (i_surface > 0)
i_surface = abs(i_surface)
if (i_surface >= OP_UNION) cycle
if (positive) then
count_positive(i_surface) = count_positive(i_surface) + 1
else
@ -853,10 +875,11 @@ contains
c => cells(i)
! loop over each surface specification
do j = 1, c % n_surfaces
do j = 1, size(c%surfaces)
i_surface = c % surfaces(j)
positive = (i_surface > 0)
i_surface = abs(i_surface)
if (i_surface >= OP_UNION) cycle
if (positive) then
count_positive(i_surface) = count_positive(i_surface) + 1

11
src/geometry_header.F90
View file

@ -128,14 +128,13 @@ module geometry_header
! the geom
integer :: material ! Material within cell (0 for
! universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: offset (:) ! Distribcell offset for tally
! counter
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell
! -- note that parentheses, union,
! etc operators will be listed here
! too
integer, allocatable :: surfaces(:) ! List of surfaces bounding cell --
! note that parentheses, union, etc
! operators will be listed here too
integer, allocatable :: rpn(:) ! Reverse Polish notation for surface
! expression
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)

10
src/initialize.F90
View file

@ -571,7 +571,7 @@ contains
! ADJUST SURFACE LIST FOR EACH CELL
c => cells(i)
do j = 1, c%n_surfaces
do j = 1, size(c%surfaces)
id = c%surfaces(j)
if (id < OP_UNION) then
if (surface_dict%has_key(abs(id))) then
@ -584,6 +584,14 @@ contains
end if
end do
do j = 1, size(c%rpn)
id = c%rpn(j)
if (id < OP_UNION) then
i_array = surface_dict%get_key(abs(id))
c%rpn(j) = sign(i_array, id)
end if
end do
! =======================================================================
! ADJUST UNIVERSE INDEX FOR EACH CELL

116
src/input_xml.F90
View file

@ -13,8 +13,9 @@ module input_xml
use plot_header
use random_lcg, only: prn
use surface_header
use stl_vector, only: VectorInt
use string, only: to_lower, to_str, str_to_int, str_to_real, &
starts_with, ends_with
starts_with, ends_with, tokenize
use tally_header, only: TallyObject, TallyFilter
use tally_initialize, only: add_tallies
use xml_interface
@ -993,6 +994,7 @@ contains
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
character(MAX_LINE_LEN) :: surface_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
class(Lattice), pointer :: lat
@ -1004,6 +1006,8 @@ contains
type(NodeList), pointer :: node_surf_list => null()
type(NodeList), pointer :: node_rlat_list => null()
type(NodeList), pointer :: node_hlat_list => null()
type(VectorInt) :: tokens
type(VectorInt) :: rpn
! Display output message
call write_message("Reading geometry XML file...", 5)
@ -1116,16 +1120,26 @@ contains
! Allocate array for surfaces and copy
if (check_for_node(node_cell, "surfaces")) then
n = get_arraysize_integer(node_cell, "surfaces")
else
n = 0
end if
c % n_surfaces = n
call get_node_value(node_cell, "surfaces", surface_spec)
if (len_trim(surface_spec) > 0) then
! Create surfaces array from string
call tokenize(surface_spec, tokens)
if (n > 0) then
allocate(c % surfaces(n))
call get_node_array(node_cell, "surfaces", c % surfaces)
! Use shunting-yard algorithm to determine RPN for surface algorithm
call generate_rpn(tokens, rpn)
! Copy surface spec and RPN form to cell arrays
allocate(c % surfaces(tokens%size()))
allocate(c % rpn(rpn%size()))
c % surfaces(:) = tokens%data(1:tokens%size())
c % rpn(:) = rpn%data(1:rpn%size())
call tokens%clear()
call rpn%clear()
end if
end if
if (.not. allocated(c%surfaces)) allocate(c%surfaces(0))
if (.not. allocated(c%rpn)) allocate(c%rpn(0))
! Rotation matrix
if (check_for_node(node_cell, "rotation")) then
@ -4766,4 +4780,88 @@ contains
end subroutine expand_natural_element
!===============================================================================
! GENERATE_RPN implements the shunting-yard algorithm to generate a Reverse
! polish notation (RPN) expression for the surface specification of a cell given
! the infix notation.
!===============================================================================
subroutine generate_rpn(tokens, output)
type(VectorInt), intent(in) :: tokens ! infix notation
type(VectorInt), intent(inout) :: output ! RPN notation
integer :: i
integer :: token
integer :: op
type(VectorInt) :: stack
do i = 1, tokens%size()
token = tokens%data(i)
if (token < OP_UNION) then
! If token is not an operator, add it to output
call output%push_back(token)
elseif (token < OP_RIGHT_PAREN) then
! Regular operators union, intersection, complement
do while (stack%size() > 0)
op = stack%data(stack%size())
if (op < OP_RIGHT_PAREN .and. &
((token == OP_COMPLEMENT .and. token < op) .or. &
(token /= OP_COMPLEMENT .and. token <= op))) then
! While there is an operator, op, on top of the stack, if the token
! is left-associative and its precedence is less than or equal to
! that of op or if the token is right-associative and its precedence
! is less than that of op, move op to the output queue and push the
! token on to the stack. Note that only complement is
! right-associative.
call output%push_back(op)
call stack%pop_back()
else
exit
end if
end do
call stack%push_back(token)
elseif (token == OP_LEFT_PAREN) then
! If the token is a left parenthesis, push it onto the stack
call stack%push_back(token)
else
! If the token is a right parenthesis, move operators from the stack to
! the output queue until reaching the left parenthesis.
do
! If we run out of operators without finding a left parenthesis, it
! means there are mismatched parentheses.
if (stack%size() == 0) then
call fatal_error('Mimatched parentheses in surface specification')
end if
op = stack%data(stack%size())
if (op == OP_LEFT_PAREN) exit
call output%push_back(op)
call stack%pop_back()
end do
! Pop the left parenthesis.
call stack%pop_back()
end if
end do
! While there are operators on the stack, move them to the output queue
do while (stack%size() > 0)
op = stack%data(stack%size())
! If the operator is a parenthesis, it is mismatched
if (op >= OP_RIGHT_PAREN) then
call fatal_error('Mimatched parentheses in surface specification')
end if
call output%push_back(op)
call stack%pop_back()
end do
end subroutine generate_rpn
end module input_xml

2
src/relaxng/geometry.rnc
View file

@ -9,7 +9,7 @@ element geometry {
(element material { ( xsd:int | "void" ) } |
attribute material { ( xsd:int | "void" ) })
) &
(element surfaces { list { xsd:int* } } | attribute surfaces { list { xsd:int* } })? &
(element surfaces { xsd:string } | attribute surfaces { xsd:string })? &
(element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? &
(element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })?
}*

12
src/relaxng/geometry.rng
View file

@ -63,18 +63,10 @@
<optional>
<choice>
<element name="surfaces">
<list>
<zeroOrMore>
<data type="int"/>
</zeroOrMore>
</list>
<data type="string"/>
</element>
<attribute name="surfaces">
<list>
<zeroOrMore>
<data type="int"/>
</zeroOrMore>
</list>
<data type="string"/>
</attribute>
</choice>
</optional>

110
src/string.F90
View file

@ -1,8 +1,10 @@
module string
use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL
use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL, &
OP_LEFT_PAREN, OP_RIGHT_PAREN, OP_COMPLEMENT, OP_INTERSECTION, OP_UNION
use error, only: fatal_error, warning
use global, only: master
use stl_vector, only: VectorInt
implicit none
@ -63,63 +65,81 @@ contains
end subroutine split_string
!===============================================================================
! SPLIT_STRING_WL takes a string that includes logical expressions for a list of
! bounding surfaces in a cell and splits it into separate words. The characters
! (, ), :, and # count as separate words since they represent operators.
!
! Arguments:
! string = input line
! words = array of words
! n = number of words
! TOKENIZE takes a string that includes logical expressions for a list of
! bounding surfaces in a cell and splits it into separate tokens. The characters
! (, ), ^, and ~ count as separate tokens since they represent operators.
!===============================================================================
subroutine split_string_wl(string, words, n)
subroutine tokenize(string, tokens)
character(*), intent(in) :: string
type(VectorInt), intent(inout) :: tokens
character(*), intent(in) :: string
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
integer :: i ! current index
integer :: i_start ! starting index of word
character(len=len_trim(string)) :: string_
character(1) :: chr ! current character
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
! Remove leading blanks
string_ = adjustl(string)
i_start = 0
i_end = 0
n = 0
do i = 1, len_trim(string)
chr = string(i:i)
i = 1
do while (i <= len_trim(string_))
! Check for special characters
if (index('():#', chr) > 0) then
if (index('()^~ ', string_(i:i)) > 0) then
! If the special character appears immediately after a non-operator,
! create a token with the surface half-space
if (i_start > 0) then
i_end = i - 1
n = n + 1
words(n) = string(i_start:i_end)
call tokens%push_back(int(str_to_int(&
string_(i_start:i - 1)), 4))
end if
n = n + 1
words(n) = chr
select case (string_(i:i))
case ('(')
call tokens%push_back(OP_LEFT_PAREN)
case (')')
call tokens%push_back(OP_RIGHT_PAREN)
case ('^')
call tokens%push_back(OP_UNION)
case ('~')
call tokens%push_back(OP_COMPLEMENT)
case (' ')
! Find next non-space character
do while (string_(i+1:i+1) == ' ')
i = i + 1
end do
! If previous token is not an operator and next token is not a left
! parenthese or union operator, that implies that the whitespace is to
! be interpreted as an intersection operator
if (i_start > 0) then
if (index(')^', string_(i+1:i+1)) == 0) then
call tokens%push_back(OP_INTERSECTION)
end if
end if
end select
i_start = 0
i_end = 0
cycle
else
! Check for invalid characters
if (index('-0123456789', string_(i:i)) == 0) then
call fatal_error("Invalid character '" // string_(i:i) // "' in &
&surface specification.")
end if
! If we haven't yet reached the start of a word, start a new word
if (i_start == 0) i_start = i
end if
if ((i_start == 0) .and. (chr /= ' ')) then
i_start = i
end if
if (i_start > 0) then
if (chr == ' ') i_end = i - 1
if (i == len_trim(string)) i_end = i
if (i_end > 0) then
n = n + 1
words(n) = string(i_start:i_end)
! reset indices
i_start = 0
i_end = 0
end if
end if
i = i + 1
end do
end subroutine split_string_wl
! If we've reached the end and we're still in a word, create a token from it
! and add it to the list
if (i_start > 0) then
call tokens%push_back(int(str_to_int(&
string_(i_start:len_trim(string_))), 4))
end if
end subroutine tokenize
!===============================================================================
! CONCATENATE takes an array of words and concatenates them together in one

31
src/summary.F90
View file

@ -107,13 +107,13 @@ contains
integer :: i, j, k, m
integer, allocatable :: lattice_universes(:,:,:)
integer, allocatable :: surface_ids(:)
integer(HID_T) :: geom_group
integer(HID_T) :: cells_group, cell_group
integer(HID_T) :: surfaces_group, surface_group
integer(HID_T) :: universes_group, univ_group
integer(HID_T) :: lattices_group, lattice_group
real(8), allocatable :: coeffs(:)
character(MAX_LINE_LEN) :: surface_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
type(Universe), pointer :: u
@ -176,15 +176,26 @@ contains
end select
! Write list of bounding surfaces
if (c%n_surfaces > 0) then
allocate(surface_ids(c%n_surfaces))
do j = 1, c%n_surfaces
k = c%surfaces(j)
surface_ids(j) = sign(surfaces(abs(k))%obj%id, k)
end do
call write_dataset(cell_group, "surfaces", surface_ids)
deallocate(surface_ids)
end if
surface_spec = ""
do j = 1, size(c%surfaces)
k = c%surfaces(j)
if (k < OP_UNION) then
surface_spec = trim(surface_spec) // " " // to_str(&
sign(surfaces(abs(k))%obj%id, k))
else
select case(k)
case (OP_LEFT_PAREN)
surface_spec = trim(surface_spec) // " ("
case (OP_RIGHT_PAREN)
surface_spec = trim(surface_spec) // " ("
case (OP_COMPLEMENT)
surface_spec = trim(surface_spec) // " ~"
case (OP_UNION)
surface_spec = trim(surface_spec) // " ^"
end select
end if
end do
call write_dataset(cell_group, "surfaces", adjustl(surface_spec))
call close_group(cell_group)
end do CELL_LOOP