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

Merge branch 'develop' into log-energy-grid

Browse source
This commit is contained in:
Paul Romano 2014-12-14 20:56:40 -05:00
commit f08330ecbb
22 changed files with 282 additions and 106 deletions
Download patch file Download diff file Expand all files Collapse all files

17
docs/source/methods/physics.rst
View file

@ -682,17 +682,20 @@ nuclear temperature, which is a function of the incoming energy of the
neutron. The ACE format contains a list of nuclear temperatures versus incoming
energies. The nuclear temperature is interpolated between neighboring incoming
energies using a specified interpolation law. Once the temperature :math:`T` is
determined, we then calculate a candidate outgoing energy based on rule C45 in
the `Monte Carlo Sampler`_:
determined, we then calculate a candidate outgoing energy based on the algorithm
given in LA-UR-14-27694_:
.. math::
:label: evaporation-E
E' = -T \log (\xi_1 \xi_2)
E' = -T \log ((1 - g\xi_1)(1 - g\xi_2))
where :math:`\xi_1, \xi_2` are random numbers sampled on the unit
interval. The outgoing energy is only accepted according to a specified
restriction energy as in equation :eq:`maxwell-restriction`.
where :math:`g = 1 - e^{-w}`, :math:`w = (E - U)/T`, :math:`U` is the
restriction energy, and :math:`\xi_1, \xi_2` are random numbers sampled on the
unit interval. The outgoing energy is only accepted according to the restriction
energy as in equation :eq:`maxwell-restriction`. This algorithm has a much
higher rejection efficiency than the standard technique, i.e. rule C45 in the
`Monte Carlo Sampler`_.
ACE Law 11 - Energy-Dependent Watt Spectrum
+++++++++++++++++++++++++++++++++++++++++++
@ -1591,6 +1594,8 @@ References
.. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721_3rdmcsampler.pdf
.. _LA-UR-14-27694: http://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-14-27694
.. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf
.. _Sutton and Brown: http://www.osti.gov/bridge/product.biblio.jsp?osti_id=307911

7
docs/source/quickinstall.rst
View file

@ -35,9 +35,9 @@ OpenMC from source as is described in :ref:`usersguide_install`.
Installing from Source on Linux or Mac OS X
-------------------------------------------
All OpenMC source code is hosted on GitHub_. If you have git_ and the gfortran_
compiler installed, you can download and install OpenMC be entering the
following commands in a terminal:
All OpenMC source code is hosted on GitHub_. If you have git_, the gfortran_
compiler, and CMake_ installed, you can download and install OpenMC be entering
the following commands in a terminal:
.. code-block:: sh
@ -58,3 +58,4 @@ can be replaced with a local install, e.g.
.. _GitHub: https://github.com/mit-crpg/openmc
.. _git: http://git-scm.com
.. _gfortran: http://gcc.gnu.org/wiki/GFortran
.. _CMake: http://www.cmake.org

19
docs/source/usersguide/input.rst
View file

@ -759,7 +759,10 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
is on the negative side of surface 3 and the positive side of surface 5, the
bounding surfaces would be given as "-3 5".
*Default*: None
.. note:: The surface attribute/element can be omitted to make a cell fill
its entire universe.
*Default*: No surfaces
:rotation:
If the cell is filled with a universe, this element specifies the angles in
@ -1172,7 +1175,7 @@ implemented in openMC:
``<plot>`` Element
------------------
Each plot must contain a combination of the following attributes or
Each plot is specified by a combination of the following attributes or
sub-elements:
:id:
@ -1193,6 +1196,18 @@ sub-elements:
*Default*: ``cell``
:level:
Universe depth to plot at (optional). This parameter controls how many
universe levels deep to pull cell and material ids from when setting plot
colors. If a given location does not have as many levels as specified,
colors will be taken from the lowest level at that location. For example, if
``level`` is set to zero colors will be taken from top-level (universe zero)
cells only. However, if ``level`` is set to 1 colors will be taken from
cells in universes that fill top-level fill-cells, and from top-level cells
that contain materials.
*Default*: Whatever the deepest universe is in the model
:origin:
Specifies the (x,y,z) coordinate of the center of the plot. Should be three
floats separated by spaces.

43
src/geometry.F90
View file

@ -11,7 +11,7 @@ module geometry
use tally, only: score_surface_current
implicit none
contains
!===============================================================================
@ -59,7 +59,8 @@ contains
end if
end do SURFACE_LOOP
! If we've reached here, then the sense matched on every surface
! 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
@ -264,13 +265,13 @@ contains
lattice_edge = .true.
end if
end if
if (lattice_edge) then
! In this case the neutron is leaving the lattice, so we move it
! out, remove all lower coordinate levels and then search from
! universe 0.
p % coord => p % coord0
call deallocate_coord(p % coord % next)
@ -287,7 +288,7 @@ contains
p % last_material = p % material
p % material = c % material
! We'll still make a new coordinate for the particle, as
! We'll still make a new coordinate for the particle, as
! distance_to_boundary will still need to track through lattice
! widths even though there's nothing in them but this material
@ -406,10 +407,9 @@ contains
! Score to global leakage tally
if (tallies_on) then
!$omp critical
!$omp atomic
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + p % wgt
!$omp end critical
end if
! Display message
@ -644,7 +644,7 @@ contains
return
end if
end if
end subroutine cross_surface
!===============================================================================
@ -906,7 +906,7 @@ contains
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
d = INFINITY
elseif (on_surface) then
! particle is on the cylinder, thus one distance is
@ -955,7 +955,7 @@ contains
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
d = INFINITY
elseif (on_surface) then
! particle is on the cylinder, thus one distance is
@ -1004,7 +1004,7 @@ contains
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
d = INFINITY
elseif (on_surface) then
! particle is on the cylinder, thus one distance is
@ -1051,7 +1051,7 @@ contains
if (quad < ZERO) then
! no intersection with sphere
d = INFINITY
d = INFINITY
elseif (on_surface) then
! particle is on the sphere, thus one distance is
@ -1098,7 +1098,7 @@ contains
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is positive/negative
@ -1117,7 +1117,7 @@ contains
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
@ -1147,7 +1147,7 @@ contains
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is positive/negative
@ -1166,7 +1166,7 @@ contains
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
@ -1196,7 +1196,7 @@ contains
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
d = INFINITY
elseif (on_surface) then
! particle is on the cone, thus one distance is positive/negative
@ -1215,7 +1215,7 @@ contains
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
@ -1273,7 +1273,7 @@ contains
! logic here checks whether the relative difference is within floating
! point precision.
if (d < dist) then
if (d < dist) then
if (abs(d - dist)/dist >= FP_REL_PRECISION) then
dist = d
if (u > 0) then
@ -1565,9 +1565,8 @@ contains
! Increment number of lost particles
p % alive = .false.
!$omp critical
!$omp atomic
n_lost_particles = n_lost_particles + 1
!$omp end critical
! Abort the simulation if the maximum number of lost particles has been
! reached

3
src/global.F90
View file

@ -118,7 +118,7 @@ module global
! 2) track-length estimate of k-eff
! 3) leakage fraction
type(TallyResult), target :: global_tallies(N_GLOBAL_TALLIES)
type(TallyResult), allocatable, target :: global_tallies(:)
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)
@ -447,6 +447,7 @@ contains
if (allocated(entropy_p)) deallocate(entropy_p)
! Deallocate tally-related arrays
if (allocated(global_tallies)) deallocate(global_tallies)
if (allocated(meshes)) deallocate(meshes)
if (allocated(tallies)) then
! First call the clear routines

1
src/initialize.F90
View file

@ -432,7 +432,6 @@ contains
i = i + 1
case ('-t', '-track', '--track')
write_all_tracks = .true.
i = i + 1
case default
call fatal_error("Unknown command line option: " // argv(i))
end select

39
src/input_xml.F90
View file

@ -1011,17 +1011,18 @@ contains
call fatal_error("Cannot specify material and fill simultaneously")
end if
! Check to make sure that surfaces were specified
if (.not. check_for_node(node_cell, "surfaces")) then
call fatal_error("No surfaces specified for cell " &
&// trim(to_str(c % id)))
end if
! Allocate array for surfaces and copy
n = get_arraysize_integer(node_cell, "surfaces")
if (check_for_node(node_cell, "surfaces")) then
n = get_arraysize_integer(node_cell, "surfaces")
else
n = 0
end if
c % n_surfaces = n
allocate(c % surfaces(n))
call get_node_array(node_cell, "surfaces", c % surfaces)
if (n > 0) then
allocate(c % surfaces(n))
call get_node_array(node_cell, "surfaces", c % surfaces)
end if
! Rotation matrix
if (check_for_node(node_cell, "rotation")) then
@ -2730,16 +2731,14 @@ contains
end select
! Set output file path
filename = "plot"
filename = trim(to_str(pl % id)) // "_plot"
if (check_for_node(node_plot, "filename")) &
call get_node_value(node_plot, "filename", filename)
select case (pl % type)
case (PLOT_TYPE_SLICE)
pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
"_" // trim(filename) // ".ppm"
pl % path_plot = trim(path_input) // trim(filename) // ".ppm"
case (PLOT_TYPE_VOXEL)
pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // &
"_" // trim(filename) // ".voxel"
pl % path_plot = trim(path_input) // trim(filename) // ".voxel"
end select
! Copy plot pixel size
@ -2819,6 +2818,18 @@ contains
end if
end if
! Copy plot cell universe level
if (check_for_node(node_plot, "level")) then
call get_node_value(node_plot, "level", pl % level)
if (pl % level < 0) then
call fatal_error("Bad universe level in plot " &
&// trim(to_str(pl % id)))
end if
else
pl % level = PLOT_LEVEL_LOWEST
end if
! Copy plot color type and initialize all colors randomly
temp_str = "cell"
if (check_for_node(node_plot, "color")) &

6
src/output.F90
View file

@ -1422,6 +1422,12 @@ contains
! Plot id
write(ou,100) "Plot ID:", trim(to_str(pl % id))
! Plot filename
write(ou,100) "Plot file:", trim(pl % path_plot)
! Plot level
write(ou,100) "Universe depth:", trim(to_str(pl % level))
! Plot type
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Plot Type:", "Slice"

24
src/physics.F90
View file

@ -256,22 +256,19 @@ contains
p % last_wgt = p % wgt
! Score implicit absorption estimate of keff
!$omp critical
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % absorb_wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
!$omp end critical
else
! See if disappearance reaction happens
if (micro_xs(i_nuclide) % absorption > &
prn() * micro_xs(i_nuclide) % total) then
! Score absorption estimate of keff
!$omp critical
!$omp atomic
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + p % wgt * &
micro_xs(i_nuclide) % nu_fission / micro_xs(i_nuclide) % absorption
!$omp end critical
p % alive = .false.
p % event = EVENT_ABSORB
@ -797,7 +794,7 @@ contains
sampling_scheme = 'cxs'
end if
! otherwise, use free gas model
! otherwise, use free gas model
else
if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then
v_target = ZERO
@ -859,7 +856,7 @@ contains
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
& maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
@ -972,7 +969,7 @@ contains
case default
call fatal_error("Not a recognized resonance scattering treatment!")
end select
end subroutine sample_target_velocity
!===============================================================================
@ -1835,14 +1832,15 @@ contains
lc = 2 + 2*NR + 2*NE
U = edist % data(lc + 1)
y = (E_in - U)/T
v = 1 - exp(-y)
! sample outgoing energy based on evaporation spectrum probability
! density function
n_sample = 0
do
r1 = prn()
r2 = prn()
E_out = -T * log(r1*r2)
if (E_out <= E_in - U) exit
x = -log((1 - v*prn())*(1 - v*prn()))
if (x <= y) exit
! check for large number of rejections
n_sample = n_sample + 1
@ -1852,6 +1850,8 @@ contains
end if
end do
E_out = x*T
case (11)
! =======================================================================
! ENERGY-DEPENDENT WATT SPECTRUM

28
src/plot.F90
View file

@ -7,7 +7,7 @@ module plot
use global
use mesh, only: get_mesh_indices
use output, only: write_message
use particle_header, only: deallocate_coord, Particle
use particle_header, only: deallocate_coord, Particle, LocalCoord
use plot_header
use ppmlib, only: Image, init_image, allocate_image, &
deallocate_image, set_pixel
@ -32,7 +32,7 @@ contains
! Display output message
call write_message("Processing plot " // trim(to_str(pl % id)) &
&// "...", 5)
&// ": " // trim(pl % path_plot) // " ...", 5)
if (pl % type == PLOT_TYPE_SLICE) then
! create 2d image
@ -58,7 +58,9 @@ contains
integer, intent(out) :: id
logical :: found_cell
type(Cell), pointer :: c => null()
integer :: level
type(Cell), pointer :: c => null()
type(LocalCoord), pointer :: coord => null()
call deallocate_coord(p % coord0 % next)
p % coord => p % coord0
@ -66,6 +68,16 @@ contains
call find_cell(p, found_cell)
if (check_overlaps) call check_cell_overlap(p)
! Loop through universes and stop on any specified level
level = 0
coord => p % coord0
do
if (level == pl % level) exit
if (.not. associated(coord % next)) exit
coord => coord % next
level = level + 1
end do
if (.not. found_cell) then
! If no cell, revert to default color
rgb = pl % not_found % rgb
@ -73,19 +85,23 @@ contains
else
if (pl % color_by == PLOT_COLOR_MATS) then
! Assign color based on material
c => cells(p % coord % cell)
c => cells(coord % cell)
if (c % material == MATERIAL_VOID) then
! By default, color void cells white
rgb = 255
id = -1
else if (c % type == CELL_FILL) then
! If we stopped on a middle universe level, treat as if not found
rgb = pl % not_found % rgb
id = -1
else
rgb = pl % colors(c % material) % rgb
id = materials(c % material) % id
end if
else if (pl % color_by == PLOT_COLOR_CELLS) then
! Assign color based on cell
rgb = pl % colors(p % coord % cell) % rgb
id = cells(p % coord % cell) % id
rgb = pl % colors(coord % cell) % rgb
id = cells(coord % cell) % id
else
rgb = 0
id = -1

4
src/plot_header.F90
View file

@ -27,6 +27,7 @@ module plot_header
integer :: basis ! direction of plot slice
integer :: pixels(3) ! pixel width/height of plot slice
integer :: meshlines_width ! pixel width of meshlines
integer :: level ! universe depth to plot the cells of
type(StructuredMesh), pointer :: meshlines_mesh => null() ! mesh to plot
type(ObjectColor) :: meshlines_color ! Color for meshlines
type(ObjectColor) :: not_found ! color for positions where no cell found
@ -36,6 +37,9 @@ module plot_header
! Plot type
integer, parameter :: PLOT_TYPE_SLICE = 1
integer, parameter :: PLOT_TYPE_VOXEL = 2
! Plot level
integer, parameter :: PLOT_LEVEL_LOWEST = -1
! Plot basis plane
integer, parameter :: PLOT_BASIS_XY = 1

6
src/relaxng/geometry.rnc
View file

@ -7,7 +7,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 { list { xsd:int* } } | attribute surfaces { list { xsd:int* } })? &
(element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? &
(element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })?
}*
@ -23,8 +23,8 @@ element geometry {
& element lattice {
(element id { xsd:int } | attribute id { xsd:int }) &
(element type { ( "rectangular" | "hexagonal" ) } |
attribute type { ( "rectangular" | "hexagonal" ) })? &
(element type { ( "rect" | "rectangle" | "rectangular" | "hexagonal" ) } |
attribute type { ( "rect" | "rectangle" | "rectangular" | "hexagonal" ) })? &
(element dimension { list { xsd:positiveInteger+ } } |
attribute dimension { list { xsd:positiveInteger+ } }) &
(element lower_left { list { xsd:double+ } } | attribute lower_left { list { xsd:double+ } }) &

1
src/relaxng/plots.rnc
View file

@ -7,6 +7,7 @@ element plots {
attribute type { "slice" | "voxel" })? &
(element color { ( "cell" | "mat" | "material" ) } |
attribute color { ( "cell" | "mat" | "material" ) })? &
(element level { xsd:int } | attribute level { xsd:int })? &
(element origin { list { xsd:double+ } } |
attribute origin { list { xsd:double+ } })? &
(element width { list { xsd:double+ } } |

45
src/tally.F90
View file

@ -281,10 +281,9 @@ contains
! get the score and tally it
score = last_wgt * calc_pn(n, mu)
!$omp critical
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
@ -347,10 +346,9 @@ contains
! get the score and tally it
score = wgt * calc_pn(n, mu)
!$omp critical
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do
j = j + t % moment_order(j)
cycle SCORE_LOOP
@ -542,10 +540,9 @@ contains
end select
! Add score to tally
!$omp critical
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -617,10 +614,9 @@ contains
i_filter = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
! Add score to tally
!$omp critical
!$omp atomic
t % results(i_score, i_filter) % value = &
t % results(i_score, i_filter) % value + score
!$omp end critical
end do
! reset outgoing energy bin and score index
@ -1015,10 +1011,9 @@ contains
end if
! Add score to tally
!$omp critical
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -1214,10 +1209,9 @@ contains
end select
! Add score to tally
!$omp critical
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -1366,10 +1360,9 @@ contains
end select
! Add score to tally
!$omp critical
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do MATERIAL_SCORE_LOOP
@ -1787,10 +1780,9 @@ contains
end if
! Add score to tally
!$omp critical
!$omp atomic
t % results(score_index, filter_index) % value = &
t % results(score_index, filter_index) % value + score
!$omp end critical
end do SCORE_LOOP
@ -2021,10 +2013,9 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
else
@ -2035,10 +2026,9 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
end if
@ -2053,10 +2043,9 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
else
@ -2067,10 +2056,9 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
end if
@ -2085,10 +2073,9 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
else
@ -2099,10 +2086,9 @@ contains
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1
!$omp critical
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
end do
end if
@ -2224,10 +2210,9 @@ contains
end if
! Add to surface current tally
!$omp critical
!$omp atomic
t % results(1, filter_index) % value = &
t % results(1, filter_index) % value + p % wgt
!$omp end critical
end if
! Calculate new coordinates

3
src/tally_initialize.F90
View file

@ -19,6 +19,9 @@ contains
subroutine configure_tallies()
! Allocate global tallies
allocate(global_tallies(N_GLOBAL_TALLIES))
call setup_tally_arrays()
call setup_tally_maps()

9
src/tracking.F90
View file

@ -62,9 +62,8 @@ contains
n_event = 0
! Add paricle's starting weight to count for normalizing tallies later
!$omp critical
!$omp atomic
total_weight = total_weight + p % wgt
!$omp end critical
! Force calculation of cross-sections by setting last energy to zero
micro_xs % last_E = ZERO
@ -112,11 +111,10 @@ contains
call score_tracklength_tally(p, distance)
! Score track-length estimate of k-eff
!$omp critical
!$omp atomic
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * &
material_xs % nu_fission
!$omp end critical
if (d_collision > d_boundary) then
! ====================================================================
@ -140,11 +138,10 @@ contains
! PARTICLE HAS COLLISION
! Score collision estimate of keff
!$omp critical
!$omp atomic
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + p % wgt * &
material_xs % nu_fission / material_xs % total
!$omp end critical
! score surface current tallies -- this has to be done before the collision
! since the direction of the particle will change and we need to use the

17
tests/test_infinite_cell/geometry.xml Normal file
View file

@ -0,0 +1,17 @@
<?xml version="1.0"?>
<geometry>
<cell id="11" universe="11" material="1"/>
<cell id="12" universe="12" material="2" surfaces=""/>
<lattice id="21" type="rect" dimension="2 2" lower_left="-2.0 -2.0"
width="2.0 2.0" outside="2">
<universes>
11 12
12 11
</universes>
</lattice>
<surface id="101" type="z-cylinder" coeffs="0.0 0.0 5.0" boundary="vacuum"/>
<cell id="101" universe="0" fill="21" surfaces="-101"/>
</geometry>

14
tests/test_infinite_cell/materials.xml Normal file
View file

@ -0,0 +1,14 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="71c" ao="1.0" />
</material>
<material id="2">
<density value="4.5" units="g/cc" />
<nuclide name="U-238" xs="71c" ao="1.0" />
</material>
</materials>

25
tests/test_infinite_cell/results.py Normal file
View file

@ -0,0 +1,25 @@
#!/usr/bin/env python
import sys
# import statepoint
sys.path.insert(0, '../../src/utils')
import statepoint
# read in statepoint file
if len(sys.argv) > 1:
sp = statepoint.StatePoint(sys.argv[1])
else:
sp = statepoint.StatePoint('statepoint.10.binary')
sp.read_results()
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

2
tests/test_infinite_cell/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
9.998895E-02 2.846817E-04

16
tests/test_infinite_cell/settings.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

59
tests/test_infinite_cell/test_infinite_cell.py Normal file
View file

@ -0,0 +1,59 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
parser = OptionParser()
parser.add_option('--mpi_exec', dest='mpi_exec', default='')
parser.add_option('--mpi_np', dest='mpi_np', default='3')
parser.add_option('--exe', dest='exe')
(opts, args) = parser.parse_args()
cwd = os.getcwd()
def test_run():
if opts.mpi_exec != '':
proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0, 'OpenMC did not exit successfully.'
def test_created_statepoint():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
'Statepoint file is not a binary or hdf5 file.'
def test_results():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree.'
def teardown():
output = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
output.append(os.path.join(cwd, 'results_test.dat'))
for f in output:
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
# test for openmc executable
if opts.exe is None:
raise Exception('Must specify OpenMC executable from command line with --exe.')
# run tests
try:
test_run()
test_created_statepoint()
test_results()
finally:
teardown()