crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Extend input so that user can specify periodic surface pairs. Also fix a bug.

Browse source
This commit is contained in:
Paul Romano 2016-05-25 07:22:43 -05:00
parent 48ac499d8a
commit 1af96416ea
11 changed files with 129 additions and 81 deletions
Download patch file Download diff file Expand all files Collapse all files

32
openmc/surface.py
View file

@ -222,6 +222,9 @@ class Plane(Surface):
boundary_type : {'transmission, 'vacuum', 'reflective'}
Boundary condition that defines the behavior for particles hitting the
surface.
periodic_surface : openmc.Surface
If a periodic boundary condition is used, the surface with which this
one is periodic with
coefficients : dict
Dictionary of surface coefficients
id : int
@ -239,6 +242,7 @@ class Plane(Surface):
self._type = 'plane'
self._coeff_keys = ['A', 'B', 'C', 'D']
self._periodic_surface = None
self.a = A
self.b = B
self.c = C
@ -260,6 +264,10 @@ class Plane(Surface):
def d(self):
return self.coefficients['D']
@property
def periodic_surface(self):
return self._periodic_surface
@a.setter
def a(self, A):
check_type('A coefficient', A, Real)
@ -280,6 +288,21 @@ class Plane(Surface):
check_type('D coefficient', D, Real)
self._coefficients['D'] = D
@periodic_surface.setter
def periodic_surface(self, periodic_surface):
check_type('periodic surface', periodic_surface, Plane)
self._periodic_surface = periodic_surface
periodic_surface._periodic_surface = self
def create_xml_subelement(self):
element = super(Plane, self).create_xml_subelement()
# Add periodic surface pair information
if self.boundary_type == 'periodic':
if self.periodic_surface is not None:
element.set("periodic_surface_id", str(self.periodic_surface.id))
return element
class XPlane(Plane):
"""A plane perpendicular to the x axis of the form :math:`x - x_0 = 0`
@ -306,6 +329,9 @@ class XPlane(Plane):
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surface.
periodic_surface : openmc.Surface
If a periodic boundary condition is used, the surface with which this
one is periodic with
coefficients : dict
Dictionary of surface coefficients
id : int
@ -391,6 +417,9 @@ class YPlane(Plane):
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surface.
periodic_surface : openmc.Surface
If a periodic boundary condition is used, the surface with which this
one is periodic with
coefficients : dict
Dictionary of surface coefficients
id : int
@ -477,6 +506,9 @@ class ZPlane(Plane):
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surface.
periodic_surface : openmc.Surface
If a periodic boundary condition is used, the surface with which this
one is periodic with
coefficients : dict
Dictionary of surface coefficients
id : int

8
src/geometry.F90
View file

@ -502,26 +502,26 @@ contains
select type (surf)
type is (SurfaceXPlane)
select type (opposite => surfaces(surf % opposite) % obj)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceXPlane)
p % coord(1) % xyz(1) = opposite % x0
end select
type is (SurfaceYPlane)
select type (opposite => surfaces(surf % opposite) % obj)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceYPlane)
p % coord(1) % xyz(2) = opposite % y0
end select
type is (SurfaceZPlane)
select type (opposite => surfaces(surf % opposite) % obj)
select type (opposite => surfaces(surf % i_periodic) % obj)
type is (SurfaceZPlane)
p % coord(1) % xyz(3) = opposite % z0
end select
end select
! Reassign particle's surface
p % surface = sign(surfaces(surf % opposite) % obj % id, p % surface)
p % surface = sign(surf % i_periodic, p % surface)
! Figure out what cell particle is in now
p % n_coord = 1

14
src/initialize.F90
View file

@ -580,20 +580,6 @@ contains
class(Lattice), pointer :: lat => null()
type(TallyObject), pointer :: t => null()
! Adjust opposite surfaces for periodic boundaries
do i = 1, size(surfaces)
associate (surf => surfaces(i) % obj)
if (surf % bc == BC_PERIODIC) then
if (surface_dict % has_key(surf % opposite)) then
surf % opposite = surface_dict % get_key(surf % opposite)
else
call fatal_error("Could not find opposite surface " // &
trim(to_str(surf % opposite)) // ".")
end if
end if
end associate
end do
do i = 1, n_cells
! =======================================================================
! ADJUST REGION SPECIFICATION FOR EACH CELL

56
src/input_xml.F90
View file

@ -1578,6 +1578,12 @@ contains
case ('periodic')
s%bc = BC_PERIODIC
boundary_exists = .true.
! Check for specification of periodic surface
if (check_for_node(node_surf, "periodic_surface_id")) then
call get_node_value(node_surf, "periodic_surface_id", &
s % i_periodic)
end if
case default
call fatal_error("Unknown boundary condition '" // trim(word) // &
&"' specified on surface " // trim(to_str(s%id)))
@ -1597,33 +1603,45 @@ contains
if (surfaces(i) % obj % bc == BC_PERIODIC) then
select type (surf => surfaces(i) % obj)
type is (SurfaceXPlane)
if (i == i_xmin) then
surf % opposite = i_xmax
elseif (i == i_xmax) then
surf % opposite = i_xmin
if (surf % i_periodic == NONE) then
if (i == i_xmin) then
surf % i_periodic = i_xmax
elseif (i == i_xmax) then
surf % i_periodic = i_xmin
else
call fatal_error("Periodic boundary condition applied to &
&interior surface.")
end if
else
call fatal_error("Periodic boundary condition applied to &
&interior surface.")
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
end if
type is (SurfaceYPlane)
if (i == i_ymin) then
surf % opposite = i_ymax
elseif (i == i_ymax) then
surf % opposite = i_ymin
if (surf % i_periodic == NONE) then
if (i == i_ymin) then
surf % i_periodic = i_ymax
elseif (i == i_ymax) then
surf % i_periodic = i_ymin
else
call fatal_error("Periodic boundary condition applied to &
&interior surface.")
end if
else
call fatal_error("Periodic boundary condition applied to &
&interior surface.")
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
end if
type is (SurfaceZPlane)
if (i == i_zmin) then
surf % opposite = i_zmax
elseif (i == i_zmax) then
surf % opposite = i_zmin
if (surf % i_periodic == NONE) then
if (i == i_zmin) then
surf % i_periodic = i_zmax
elseif (i == i_zmax) then
surf % i_periodic = i_zmin
else
call fatal_error("Periodic boundary condition applied to &
&interior surface.")
end if
else
call fatal_error("Periodic boundary condition applied to &
&interior surface.")
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
end if
class default
@ -1633,7 +1651,7 @@ contains
! Make sure opposite surface is also periodic
associate (surf => surfaces(i) % obj)
if (surfaces(surf % opposite) % obj % bc /= BC_PERIODIC) then
if (surfaces(surf % i_periodic) % obj % bc /= BC_PERIODIC) then
call fatal_error("Could not find matching surface for periodic &
&boundary on surface " // trim(to_str(surf % id)) // ".")
end if

6
src/surface_header.F90
View file

@ -1,6 +1,6 @@
module surface_header
use constants, only: ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT
use constants, only: NONE, ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT
implicit none
@ -15,8 +15,8 @@ module surface_header
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
integer :: opposite ! Opposite surface for periodic boundary
character(len=104) :: name = "" ! User-defined name
integer :: i_periodic = NONE ! Index of corresponding periodic surface
character(len=104) :: name = "" ! User-defined name
contains
procedure :: sense
procedure :: reflect

12
tests/test_periodic/geometry.xml
View file

@ -1,12 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<surface id="1" type="x-plane" coeffs="-5." boundary="periodic" />
<surface id="2" type="x-plane" coeffs="5." boundary="periodic" />
<surface id="3" type="y-plane" coeffs="-5." boundary="periodic" />
<surface id="4" type="y-plane" coeffs="5." boundary="periodic" />
<surface id="5" type="z-plane" coeffs="-5." boundary="reflective" />
<surface id="6" type="z-plane" coeffs="5." boundary="reflective" />
<surface id="7" type="z-cylinder" coeffs="-2.5 2.5 2.0" />
<cell id="1" material="1" region="1 -2 3 -4 5 -6 7" />
<cell id="2" material="2" region="5 -6 -7" />
</geometry>

1
tests/test_periodic/inputs_true.dat Normal file
View file

@ -0,0 +1 @@
af589996f2930337afe34ba9894098ff5efe3b29b6e927117220b718bf29b630ffdbc931754d465a8e8100125a8aa997dbe10aab322b43f69d59710573996a6d

13
tests/test_periodic/materials.xml
View file

@ -1,13 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="1">
<density units="g/cc" value="1.0" />
<nuclide name="H-1" xs="71c" ao="2.0" />
<nuclide name="O-16" xs="71c" ao="2.0" />
<sab name="HH2O" xs="71t" />
</material>
<material id="2">
<density units="g/cc" value="4.5" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
</materials>

2
tests/test_periodic/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.542742E+00 4.410461E-02
1.040109E+00 6.527490E-02

13
tests/test_periodic/settings.xml
View file

@ -1,13 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<settings>
<eigenvalue>
<particles>1000</particles>
<batches>4</batches>
<inactive>0</inactive>
</eigenvalue>
<source strength="1.0">
<space type="box">
<parameters>-5. -5. -5. 5. 5. 5.</parameters>
</space>
</source>
</settings>

53
tests/test_periodic/test_periodic.py
View file

@ -3,9 +3,58 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
from testing_harness import PyAPITestHarness
import openmc
class PeriodicTest(PyAPITestHarness):
def _build_inputs(self):
# Define materials
water = openmc.Material(1)
water.add_nuclide('H-1', 2.0)
water.add_nuclide('O-16', 1.0)
water.add_s_alpha_beta('HH2O', '71t')
water.set_density('g/cc', 1.0)
fuel = openmc.Material(2)
fuel.add_nuclide('U-235', 1.0)
fuel.set_density('g/cc', 4.5)
materials = openmc.Materials((water, fuel))
materials.default_xs = '71c'
materials.export_to_xml()
# Define geometry
x_min = openmc.XPlane(1, x0=-5., boundary_type='periodic')
x_max = openmc.XPlane(2, x0=5., boundary_type='periodic')
x_max.periodic_surface = x_min
y_min = openmc.YPlane(3, y0=-5., boundary_type='periodic')
y_max = openmc.YPlane(4, y0=5., boundary_type='periodic')
z_min = openmc.ZPlane(5, z0=-5., boundary_type='reflective')
z_max = openmc.ZPlane(6, z0=5., boundary_type='reflective')
z_cyl = openmc.ZCylinder(7, x0=-2.5, y0=2.5, R=2.0)
outside_cyl = openmc.Cell(1, fill=water, region=(
+x_min & -x_max & +y_min & -y_max & +z_min & -z_max & +z_cyl))
inside_cyl = openmc.Cell(2, fill=fuel, region=+z_min & -z_max & -z_cyl)
root_universe = openmc.Universe(0, cells=(outside_cyl, inside_cyl))
geometry = openmc.Geometry()
geometry.root_universe = root_universe
geometry.export_to_xml()
# Define settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 4
settings.inactive = 0
settings.source = openmc.Source(space=openmc.stats.Box(
*outside_cyl.region.bounding_box))
settings.export_to_xml()
if __name__ == '__main__':
harness = TestHarness('statepoint.4.h5')
harness = PeriodicTest('statepoint.4.h5')
harness.main()