From 42ef7ccfac5febf1d62e44bd8074dc15222c0974 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 23 Mar 2016 15:47:51 -0500 Subject: [PATCH] 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 775407d70b..ea51723b7b 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 8a38f982b5..9f77817387 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 09bedb1388..74cab87c4b 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 90c703d27d..9e49e3c349 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 8d25789f5a..6cb6f7c158 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 d40401b281..3ff0f67c69 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 4686552176..0b5d3c86be 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 0000000000..6ecfec1972 --- /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 0000000000..a7bf4faf4e --- /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 0000000000..f65dbafd1b --- /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 0000000000..af09407ae5 --- /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 0000000000..595d7c0ddd --- /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 0000000000..b584632f08 --- /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 78e5553e8c..e659768856 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)