From 42ef7ccfac5febf1d62e44bd8074dc15222c0974 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 23 Mar 2016 15:47:51 -0500 Subject: [PATCH 1/6] Add periodic boundary conditions --- docs/source/usersguide/input.rst | 2 +- src/geometry.F90 | 62 ++++++++++++++++++++++++ src/initialize.F90 | 14 ++++++ src/input_xml.F90 | 72 +++++++++++++++++++++++++++- src/relaxng/geometry.rnc | 5 +- src/relaxng/geometry.rng | 12 +++++ src/surface_header.F90 | 1 + tests/test_periodic/geometry.xml | 12 +++++ tests/test_periodic/materials.xml | 13 +++++ tests/test_periodic/results_true.dat | 2 + tests/test_periodic/settings.xml | 13 +++++ tests/test_periodic/tallies.xml | 14 ++++++ tests/test_periodic/test_periodic.py | 11 +++++ tests/testing_harness.py | 3 +- 14 files changed, 231 insertions(+), 5 deletions(-) create mode 100644 tests/test_periodic/geometry.xml create mode 100644 tests/test_periodic/materials.xml create mode 100644 tests/test_periodic/results_true.dat create mode 100644 tests/test_periodic/settings.xml create mode 100644 tests/test_periodic/tallies.xml create mode 100644 tests/test_periodic/test_periodic.py diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 775407d70..ea51723b7 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -898,7 +898,7 @@ Each ```` element can have the following attributes or sub-elements: :boundary: The boundary condition for the surface. This can be "transmission", - "vacuum", or "reflective". + "vacuum", "reflective", or "periodic". *Default*: "transmission" diff --git a/src/geometry.F90 b/src/geometry.F90 index 8a38f982b..9f7781738 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -475,6 +475,68 @@ contains &// trim(to_str(surf%id))) end if return + elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then + ! ======================================================================= + ! PERIODIC BOUNDARY + + ! Do not handle periodic boundary conditions on lower universes + if (p % n_coord /= 1) then + call handle_lost_particle(p, "Cannot period particle " & + // trim(to_str(p % id)) // " off surface in a lower universe.") + return + end if + + ! Score surface currents since reflection causes the direction of the + ! particle to change -- artificially move the particle slightly back in + ! case the surface crossing in coincident with a mesh boundary + + if (active_current_tallies % size() > 0) then + p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw + call score_surface_current(p) + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + end if + + select type (surf) + type is (SurfaceXPlane) + select type (opposite => surfaces(surf % opposite) % obj) + type is (SurfaceXPlane) + p % coord(1) % xyz(1) = opposite % x0 + end select + + type is (SurfaceYPlane) + select type (opposite => surfaces(surf % opposite) % obj) + type is (SurfaceYPlane) + p % coord(1) % xyz(2) = opposite % y0 + end select + + type is (SurfaceZPlane) + select type (opposite => surfaces(surf % opposite) % 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) + + ! Figure out what cell particle is in now + p % n_coord = 1 + call find_cell(p, found) + if (.not. found) then + call handle_lost_particle(p, "Couldn't find particle after hitting & + &periodic boundary on surface " // trim(to_str(surf%id)) // ".") + return + end if + + ! Set previous coordinate going slightly past surface crossing + p % last_xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + + ! Diagnostic message + if (verbosity >= 10 .or. trace) then + call write_message(" Hit periodic boundary on surface " & + // trim(to_str(surf%id))) + end if + return end if ! ========================================================================== diff --git a/src/initialize.F90 b/src/initialize.F90 index 09bedb138..74cab87c4 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -580,6 +580,20 @@ 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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 90c703d27..9e49e3c34 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1109,6 +1109,8 @@ contains integer :: universe_num integer :: n_cells_in_univ integer :: coeffs_reqd + integer :: i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax + real(8) :: xmin, xmax, ymin, ymax, zmin, zmax integer, allocatable :: temp_int_array(:) real(8) :: phi, theta, psi real(8), allocatable :: coeffs(:) @@ -1387,6 +1389,13 @@ contains call fatal_error("No surfaces found in geometry.xml!") end if + xmin = INFINITY + xmax = -INFINITY + ymin = INFINITY + ymax = -INFINITY + zmin = INFINITY + zmax = -INFINITY + ! Allocate cells array allocate(surfaces(n_surfaces)) @@ -1478,10 +1487,28 @@ contains select type(s) type is (SurfaceXPlane) s%x0 = coeffs(1) + + ! Determine outer surfaces + xmin = min(xmin, s % x0) + xmax = max(xmax, s % x0) + if (xmin == s % x0) i_xmin = i + if (xmax == s % x0) i_xmax = i type is (SurfaceYPlane) s%y0 = coeffs(1) + + ! Determine outer surfaces + ymin = min(ymin, s % y0) + ymax = max(ymax, s % y0) + if (ymin == s % y0) i_ymin = i + if (ymax == s % y0) i_ymax = i type is (SurfaceZPlane) s%z0 = coeffs(1) + + ! Determine outer surfaces + zmin = min(zmin, s % z0) + zmax = max(zmax, s % z0) + if (zmin == s % z0) i_zmin = i + if (zmax == s % z0) i_zmax = i type is (SurfacePlane) s%A = coeffs(1) s%B = coeffs(2) @@ -1548,11 +1575,13 @@ contains case ('reflective', 'reflect', 'reflecting') s%bc = BC_REFLECT boundary_exists = .true. + case ('periodic') + s%bc = BC_PERIODIC + boundary_exists = .true. case default call fatal_error("Unknown boundary condition '" // trim(word) // & &"' specified on surface " // trim(to_str(s%id))) end select - ! Add surface to dictionary call surface_dict % add_key(s%id, i) end do @@ -1563,6 +1592,47 @@ contains call fatal_error("No boundary conditions were applied to any surfaces!") end if + ! Determine opposite side for periodic boundaries + do i = 1, size(surfaces) + 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 + else + call fatal_error("Periodic boundary condition applied to & + &interior surface.") + end if + + type is (SurfaceYPlane) + if (i == i_ymin) then + surf % opposite = i_ymax + elseif (i == i_ymax) then + surf % opposite = i_ymin + else + call fatal_error("Periodic boundary condition applied to & + &interior surface.") + end if + + type is (SurfaceZPlane) + if (i == i_zmin) then + surf % opposite = i_zmax + elseif (i == i_zmax) then + surf % opposite = i_zmin + else + call fatal_error("Periodic boundary condition applied to & + &interior surface.") + end if + + class default + call fatal_error("Periodic boundary condition applied to & + &non-planar surface.") + end select + end if + end do + ! ========================================================================== ! READ LATTICES FROM GEOMETRY.XML diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index 8d25789f5..6cb6f7c15 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -21,8 +21,9 @@ element geometry { (element type { xsd:string { maxLength = "15" } } | attribute type { xsd:string { maxLength = "15" } }) & (element coeffs { list { xsd:double+ } } | attribute coeffs { list { xsd:double+ } }) & - (element boundary { ( "transmit" | "reflective" | "vacuum" ) } | - attribute boundary { ( "transmit" | "reflective" | "vacuum" ) })? + (element boundary { ( "transmit" | "reflective" | "vacuum" | "periodic" ) } | + attribute boundary { ( "transmit" | "reflective" | "vacuum" | "periodic" ) })? & + (element opposite { xsd:int } | attribute opposite { xsd:int })? }* & element lattice { diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index d40401b28..3ff0f67c6 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -173,6 +173,7 @@ transmit reflective vacuum + periodic @@ -180,10 +181,21 @@ transmit reflective vacuum + periodic + + + + + + + + + + diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 468655217..0b5d3c86b 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -15,6 +15,7 @@ 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 contains procedure :: sense diff --git a/tests/test_periodic/geometry.xml b/tests/test_periodic/geometry.xml new file mode 100644 index 000000000..6ecfec197 --- /dev/null +++ b/tests/test_periodic/geometry.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/tests/test_periodic/materials.xml b/tests/test_periodic/materials.xml new file mode 100644 index 000000000..a7bf4faf4 --- /dev/null +++ b/tests/test_periodic/materials.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/tests/test_periodic/results_true.dat b/tests/test_periodic/results_true.dat new file mode 100644 index 000000000..f65dbafd1 --- /dev/null +++ b/tests/test_periodic/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.542742E+00 4.410461E-02 diff --git a/tests/test_periodic/settings.xml b/tests/test_periodic/settings.xml new file mode 100644 index 000000000..af09407ae --- /dev/null +++ b/tests/test_periodic/settings.xml @@ -0,0 +1,13 @@ + + + + 1000 + 4 + 0 + + + + -5. -5. -5. 5. 5. 5. + + + diff --git a/tests/test_periodic/tallies.xml b/tests/test_periodic/tallies.xml new file mode 100644 index 000000000..595d7c0dd --- /dev/null +++ b/tests/test_periodic/tallies.xml @@ -0,0 +1,14 @@ + + + + regular + -200. -1e50 + 200. 1e50 + 50 1 + + + collision + + fission + + diff --git a/tests/test_periodic/test_periodic.py b/tests/test_periodic/test_periodic.py new file mode 100644 index 000000000..b584632f0 --- /dev/null +++ b/tests/test_periodic/test_periodic.py @@ -0,0 +1,11 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness + + +if __name__ == '__main__': + harness = TestHarness('statepoint.4.h5') + harness.main() diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 78e5553e8..e65976885 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -133,9 +133,10 @@ class TestHarness(object): def _cleanup(self): """Delete statepoints, tally, and test files.""" - output = glob.glob(os.path.join(os.getcwd(), 'statepoint.*.*')) + output = glob.glob(os.path.join(os.getcwd(), 'statepoint.*.h5')) output.append(os.path.join(os.getcwd(), 'tallies.out')) output.append(os.path.join(os.getcwd(), 'results_test.dat')) + output.append(os.path.join(os.getcwd(), 'summary.h5')) for f in output: if os.path.exists(f): os.remove(f) From 5e910498b39b12ab9fe7ae5ef4a14a80a28940d7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 11 May 2016 21:39:56 -0500 Subject: [PATCH 2/6] Check for non-matching periodic boundary conditions --- src/geometry.F90 | 4 ++-- src/input_xml.F90 | 8 ++++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 9f7781738..62c5036a9 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -432,7 +432,7 @@ contains ! Score surface currents since reflection causes the direction of the ! particle to change -- artificially move the particle slightly back in - ! case the surface crossing in coincident with a mesh boundary + ! case the surface crossing is coincident with a mesh boundary if (active_current_tallies % size() > 0) then p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw @@ -488,7 +488,7 @@ contains ! Score surface currents since reflection causes the direction of the ! particle to change -- artificially move the particle slightly back in - ! case the surface crossing in coincident with a mesh boundary + ! case the surface crossing is coincident with a mesh boundary if (active_current_tallies % size() > 0) then p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9e49e3c34..242d3459a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1630,6 +1630,14 @@ contains call fatal_error("Periodic boundary condition applied to & &non-planar surface.") end select + + ! Make sure opposite surface is also periodic + associate (surf => surfaces(i) % obj) + if (surfaces(surf % opposite) % 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 + end associate end if end do From 48ac499d8a96ccc2ceadc72675e0b8a0b7df3d61 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 24 May 2016 16:24:26 -0500 Subject: [PATCH 3/6] Respond to @smharper suggestions on #647 --- docs/source/usersguide/input.rst | 8 +++-- openmc/surface.py | 53 +++++++++++++++++--------------- src/geometry.F90 | 12 +++++--- 3 files changed, 43 insertions(+), 30 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index ea51723b7..e9ce42fd9 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -897,8 +897,12 @@ Each ```` element can have the following attributes or sub-elements: *Default*: None :boundary: - The boundary condition for the surface. This can be "transmission", - "vacuum", "reflective", or "periodic". + The boundary condition for the surface. This can be "transmission", + "vacuum", "reflective", or "periodic". Periodic boundary conditions can + only be applied to x-, y-, and z-planes. Only axis-aligned periodicity is + supported, i.e., x-planes an only be paired with x-planes. Specify which + planes are periodic and the code will automatically identify which planes + are paired together. *Default*: "transmission" diff --git a/openmc/surface.py b/openmc/surface.py index 37e7c2ffd..239b868aa 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -37,7 +37,9 @@ class Surface(object): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles - freely pass through the surface. + freely pass through the surface. Note that periodic boundary conditions + can only be applied to x-, y-, and z-planes, and only axis-aligned + periodicity is supported. name : str, optional Name of the surface. If not specified, the name will be the empty string. @@ -192,7 +194,7 @@ class Plane(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -217,7 +219,7 @@ class Plane(Surface): The 'C' parameter for the plane d : float The 'D' parameter for the plane - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -290,7 +292,8 @@ class XPlane(Plane): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles - freely pass through the surface. + freely pass through the surface. Only axis-aligned periodicity is + supported, i.e., x-planes can only be paired with x-planes. x0 : float, optional Location of the plane. Defaults to 0. name : str, optional @@ -374,7 +377,8 @@ class YPlane(Plane): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles - freely pass through the surface. + freely pass through the surface. Only axis-aligned periodicity is + supported, i.e., x-planes can only be paired with x-planes. y0 : float, optional Location of the plane name : str, optional @@ -459,7 +463,8 @@ class ZPlane(Plane): boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles - freely pass through the surface. + freely pass through the surface. Only axis-aligned periodicity is + supported, i.e., x-planes can only be paired with x-planes. z0 : float, optional Location of the plane. Defaults to 0. name : str, optional @@ -541,7 +546,7 @@ class Cylinder(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -555,7 +560,7 @@ class Cylinder(Surface): ---------- r : float Radius of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -597,7 +602,7 @@ class XCylinder(Cylinder): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -617,7 +622,7 @@ class XCylinder(Cylinder): y-coordinate of the center of the cylinder z0 : float z-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -700,7 +705,7 @@ class YCylinder(Cylinder): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -720,7 +725,7 @@ class YCylinder(Cylinder): x-coordinate of the center of the cylinder z0 : float z-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -803,7 +808,7 @@ class ZCylinder(Cylinder): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -823,7 +828,7 @@ class ZCylinder(Cylinder): x-coordinate of the center of the cylinder y0 : float y-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -905,7 +910,7 @@ class Sphere(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -930,7 +935,7 @@ class Sphere(Surface): z-coordinate of the center of the sphere R : float Radius of the sphere - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1033,7 +1038,7 @@ class Cone(Surface): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1058,7 +1063,7 @@ class Cone(Surface): z-coordinate of the apex R2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1130,7 +1135,7 @@ class XCone(Cone): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1155,7 +1160,7 @@ class XCone(Cone): z-coordinate of the apex R2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1186,7 +1191,7 @@ class YCone(Cone): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1211,7 +1216,7 @@ class YCone(Cone): z-coordinate of the apex R2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -1242,7 +1247,7 @@ class ZCone(Cone): surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -1267,7 +1272,7 @@ class ZCone(Cone): z-coordinate of the apex R2 : float Parameter related to the aperature - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + boundary_type : {'transmission, 'vacuum', 'reflective'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict diff --git a/src/geometry.F90 b/src/geometry.F90 index 62c5036a9..2059ea42f 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -378,6 +378,7 @@ contains real(8) :: v ! y-component of direction real(8) :: w ! z-component of direction real(8) :: norm ! "norm" of surface normal + real(8) :: xyz(3) ! Saved global coordinate integer :: i_surface ! index in surfaces logical :: found ! particle found in universe? class(Surface), pointer :: surf @@ -435,9 +436,10 @@ contains ! case the surface crossing is coincident with a mesh boundary if (active_current_tallies % size() > 0) then + xyz = p % coord(1) % xyz p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw call score_surface_current(p) - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + p % coord(1) % xyz = xyz end if ! Reflect particle off surface @@ -481,8 +483,9 @@ contains ! Do not handle periodic boundary conditions on lower universes if (p % n_coord /= 1) then - call handle_lost_particle(p, "Cannot period particle " & - // trim(to_str(p % id)) // " off surface in a lower universe.") + call handle_lost_particle(p, "Cannot transfer particle " & + // trim(to_str(p % id)) // " across surface in a lower universe.& + & Boundary conditions must be applied to universe 0.") return end if @@ -491,9 +494,10 @@ contains ! case the surface crossing is coincident with a mesh boundary if (active_current_tallies % size() > 0) then + xyz = p % coord(1) % xyz p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw call score_surface_current(p) - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + p % coord(1) % xyz = xyz end if select type (surf) From 1af96416eaa502d714076dfd85f68a59b01b52de Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 25 May 2016 07:22:43 -0500 Subject: [PATCH 4/6] Extend input so that user can specify periodic surface pairs. Also fix a bug. --- openmc/surface.py | 32 ++++++++++++++++ src/geometry.F90 | 8 ++-- src/initialize.F90 | 14 ------- src/input_xml.F90 | 56 ++++++++++++++++++---------- src/surface_header.F90 | 6 +-- tests/test_periodic/geometry.xml | 12 ------ tests/test_periodic/inputs_true.dat | 1 + tests/test_periodic/materials.xml | 13 ------- tests/test_periodic/results_true.dat | 2 +- tests/test_periodic/settings.xml | 13 ------- tests/test_periodic/test_periodic.py | 53 +++++++++++++++++++++++++- 11 files changed, 129 insertions(+), 81 deletions(-) delete mode 100644 tests/test_periodic/geometry.xml create mode 100644 tests/test_periodic/inputs_true.dat delete mode 100644 tests/test_periodic/materials.xml delete mode 100644 tests/test_periodic/settings.xml diff --git a/openmc/surface.py b/openmc/surface.py index 239b868aa..ca2d5d6b4 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -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 diff --git a/src/geometry.F90 b/src/geometry.F90 index 2059ea42f..21c8baa6b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 74cab87c4..09bedb138 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 242d3459a..8a9b2299f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 0b5d3c86b..68e5144b7 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -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 diff --git a/tests/test_periodic/geometry.xml b/tests/test_periodic/geometry.xml deleted file mode 100644 index 6ecfec197..000000000 --- a/tests/test_periodic/geometry.xml +++ /dev/null @@ -1,12 +0,0 @@ - - - - - - - - - - - - diff --git a/tests/test_periodic/inputs_true.dat b/tests/test_periodic/inputs_true.dat new file mode 100644 index 000000000..d50d0b859 --- /dev/null +++ b/tests/test_periodic/inputs_true.dat @@ -0,0 +1 @@ +af589996f2930337afe34ba9894098ff5efe3b29b6e927117220b718bf29b630ffdbc931754d465a8e8100125a8aa997dbe10aab322b43f69d59710573996a6d \ No newline at end of file diff --git a/tests/test_periodic/materials.xml b/tests/test_periodic/materials.xml deleted file mode 100644 index a7bf4faf4..000000000 --- a/tests/test_periodic/materials.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - - - - - - - - - - diff --git a/tests/test_periodic/results_true.dat b/tests/test_periodic/results_true.dat index f65dbafd1..b0bdb22c7 100644 --- a/tests/test_periodic/results_true.dat +++ b/tests/test_periodic/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.542742E+00 4.410461E-02 +1.040109E+00 6.527490E-02 diff --git a/tests/test_periodic/settings.xml b/tests/test_periodic/settings.xml deleted file mode 100644 index af09407ae..000000000 --- a/tests/test_periodic/settings.xml +++ /dev/null @@ -1,13 +0,0 @@ - - - - 1000 - 4 - 0 - - - - -5. -5. -5. 5. 5. 5. - - - diff --git a/tests/test_periodic/test_periodic.py b/tests/test_periodic/test_periodic.py index b584632f0..558514575 100644 --- a/tests/test_periodic/test_periodic.py +++ b/tests/test_periodic/test_periodic.py @@ -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() From 63e355f5a02730d6f2d1c4dbf5373d4e9dea7395 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 25 May 2016 10:09:02 -0500 Subject: [PATCH 5/6] Fix typo in documentation pointed out by @smharper --- docs/source/usersguide/input.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index e9ce42fd9..da9896fbb 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -900,7 +900,7 @@ Each ```` element can have the following attributes or sub-elements: The boundary condition for the surface. This can be "transmission", "vacuum", "reflective", or "periodic". Periodic boundary conditions can only be applied to x-, y-, and z-planes. Only axis-aligned periodicity is - supported, i.e., x-planes an only be paired with x-planes. Specify which + supported, i.e., x-planes can only be paired with x-planes. Specify which planes are periodic and the code will automatically identify which planes are paired together. From 55740283934463d2e58255b1a5956c5afc757ffd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 25 May 2016 10:33:47 -0500 Subject: [PATCH 6/6] Update documentation and fix Relax NG schema --- docs/source/usersguide/input.rst | 4 ++++ src/relaxng/geometry.rnc | 2 +- src/relaxng/geometry.rng | 4 ++-- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index da9896fbb..d1a01b65b 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -906,6 +906,10 @@ Each ```` element can have the following attributes or sub-elements: *Default*: "transmission" + :periodic_surface_id: + If a periodic boundary condition is applied, this attribute identifies the + ``id`` of the corresponding periodic sufrace. + The following quadratic surfaces can be modeled: :x-plane: diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index 6cb6f7c15..35d5ef8b2 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -23,7 +23,7 @@ element geometry { (element coeffs { list { xsd:double+ } } | attribute coeffs { list { xsd:double+ } }) & (element boundary { ( "transmit" | "reflective" | "vacuum" | "periodic" ) } | attribute boundary { ( "transmit" | "reflective" | "vacuum" | "periodic" ) })? & - (element opposite { xsd:int } | attribute opposite { xsd:int })? + (element periodic_surface_id { xsd:int } | attribute periodic_surface_id { xsd:int })? }* & element lattice { diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index 3ff0f67c6..b53d0e8db 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -188,10 +188,10 @@ - + - +