From 996a03d5ab90d5a13af442afe7e9fab5421a3aa7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 25 Oct 2014 16:37:30 -0400 Subject: [PATCH 01/13] Add new evaporation spectrum sampling method outlined in LA-UR-14-27694. Still need documentation. --- src/physics.F90 | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 8fa0575cd..d136d7306 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -797,7 +797,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 +859,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 +972,7 @@ contains case default call fatal_error("Not a recognized resonance scattering treatment!") end select - + end subroutine sample_target_velocity !=============================================================================== @@ -1835,14 +1835,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 +1853,8 @@ contains end if end do + E_out = x*T + case (11) ! ======================================================================= ! ENERGY-DEPENDENT WATT SPECTRUM From 332970ea2e03d2512e75575c9939b50b8d576f31 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Oct 2014 21:53:39 -0400 Subject: [PATCH 02/13] Update documentation for evaporation spectrum. --- docs/source/methods/physics.rst | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index 54dc91345..db6cfd89e 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -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 From 26fb936f2a086ec0a6157a9eedb2d6adfbb15936 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Nov 2014 22:32:18 -0500 Subject: [PATCH 03/13] Fix problem with -t, --track command-line flag. --- src/initialize.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index e0a5bf8b0..a418b3ec3 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -434,7 +434,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 From a704b412ee11436a40796dff49253474b9f2b2f1 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 12 Nov 2014 19:24:58 -0500 Subject: [PATCH 04/13] Added support for plot level depth, removed id from filename --- src/input_xml.F90 | 18 ++++++++++++++---- src/output.F90 | 6 ++++++ src/plot.F90 | 28 ++++++++++++++++++++++------ src/plot_header.F90 | 4 ++++ 4 files changed, 46 insertions(+), 10 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e992a0eae..340e430c7 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2735,11 +2735,9 @@ contains 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 +2817,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")) & diff --git a/src/output.F90 b/src/output.F90 index 4d0ef7edd..ce5059feb 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1428,6 +1428,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" diff --git a/src/plot.F90 b/src/plot.F90 index 1297b746c..ddc4641b2 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -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 diff --git a/src/plot_header.F90 b/src/plot_header.F90 index 68eb89a11..7ce0aa450 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -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 From 94d3652c36b1ffe7c16f7add3ceffeb08efdf74f Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 12 Nov 2014 19:25:22 -0500 Subject: [PATCH 05/13] Updated documentation --- docs/source/usersguide/input.rst | 14 +++++++++++++- src/relaxng/plots.rnc | 1 + 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index f7e70ef20..73d2e7f86 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1170,7 +1170,7 @@ implemented in openMC: ```` 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: @@ -1191,6 +1191,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 lcation. 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. diff --git a/src/relaxng/plots.rnc b/src/relaxng/plots.rnc index 27b2ae7f7..5ef14f229 100644 --- a/src/relaxng/plots.rnc +++ b/src/relaxng/plots.rnc @@ -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+ } } | From 4bda190ee8bd5196ff05e595234d108137e31617 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 12 Nov 2014 19:29:05 -0500 Subject: [PATCH 06/13] Fixed plot filename issue --- src/input_xml.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 340e430c7..2adda87a7 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2730,7 +2730,7 @@ 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) From 4dbca051f74a2ed62902d44a99eaf1b89420b6a6 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 13 Nov 2014 02:43:16 -0500 Subject: [PATCH 07/13] Allow cells with no surfaces. --- src/geometry.F90 | 3 +- src/input_xml.F90 | 19 +++--- src/relaxng/geometry.rnc | 6 +- tests/test_infinite_cell/geometry.xml | 17 ++++++ tests/test_infinite_cell/materials.xml | 14 +++++ tests/test_infinite_cell/results.py | 25 ++++++++ tests/test_infinite_cell/results_true.dat | 2 + tests/test_infinite_cell/settings.xml | 16 +++++ .../test_infinite_cell/test_infinite_cell.py | 59 +++++++++++++++++++ 9 files changed, 148 insertions(+), 13 deletions(-) create mode 100644 tests/test_infinite_cell/geometry.xml create mode 100644 tests/test_infinite_cell/materials.xml create mode 100644 tests/test_infinite_cell/results.py create mode 100644 tests/test_infinite_cell/results_true.dat create mode 100644 tests/test_infinite_cell/settings.xml create mode 100644 tests/test_infinite_cell/test_infinite_cell.py diff --git a/src/geometry.F90 b/src/geometry.F90 index 043ded05b..2d9f3355a 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e992a0eae..928e78911 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index b7975d05f..e7f81c485 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -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:intr* } })? & (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+ } }) & diff --git a/tests/test_infinite_cell/geometry.xml b/tests/test_infinite_cell/geometry.xml new file mode 100644 index 000000000..77ef6110a --- /dev/null +++ b/tests/test_infinite_cell/geometry.xml @@ -0,0 +1,17 @@ + + + + + + + + 11 12 + 12 11 + + + + + + + diff --git a/tests/test_infinite_cell/materials.xml b/tests/test_infinite_cell/materials.xml new file mode 100644 index 000000000..2e5b48381 --- /dev/null +++ b/tests/test_infinite_cell/materials.xml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + diff --git a/tests/test_infinite_cell/results.py b/tests/test_infinite_cell/results.py new file mode 100644 index 000000000..be13ee66f --- /dev/null +++ b/tests/test_infinite_cell/results.py @@ -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) diff --git a/tests/test_infinite_cell/results_true.dat b/tests/test_infinite_cell/results_true.dat new file mode 100644 index 000000000..45eaa4f91 --- /dev/null +++ b/tests/test_infinite_cell/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +9.998895E-02 2.846817E-04 diff --git a/tests/test_infinite_cell/settings.xml b/tests/test_infinite_cell/settings.xml new file mode 100644 index 000000000..a6fd5da19 --- /dev/null +++ b/tests/test_infinite_cell/settings.xml @@ -0,0 +1,16 @@ + + + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/test_infinite_cell/test_infinite_cell.py b/tests/test_infinite_cell/test_infinite_cell.py new file mode 100644 index 000000000..6fdbf8745 --- /dev/null +++ b/tests/test_infinite_cell/test_infinite_cell.py @@ -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() From dad06b49fc35cb216aa1629fa6b5dde605608a95 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 13 Nov 2014 03:55:16 -0500 Subject: [PATCH 08/13] Updated docs for cells with no surfaces --- docs/source/usersguide/input.rst | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index f7e70ef20..fe495b117 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -757,7 +757,10 @@ Each ```` 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 From 991cdb9cccf7c07ac4e1d22a169222d42c43171c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 14 Nov 2014 16:54:01 -0500 Subject: [PATCH 09/13] Change !$omp critical directives to !$omp atomic where possible. --- src/geometry.F90 | 40 +++++++++++++++++++--------------------- src/physics.F90 | 13 +++++-------- src/tally.F90 | 45 +++++++++++++++------------------------------ src/tracking.F90 | 9 +++------ 4 files changed, 42 insertions(+), 65 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 043ded05b..3e124983a 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -11,7 +11,7 @@ module geometry use tally, only: score_surface_current implicit none - + contains !=============================================================================== @@ -264,13 +264,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 +287,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 +406,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 +643,7 @@ contains return end if end if - + end subroutine cross_surface !=============================================================================== @@ -906,7 +905,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 +954,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 +1003,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 +1050,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 +1097,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 +1116,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 +1146,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 +1165,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 +1195,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 +1214,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 +1272,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 +1564,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 diff --git a/src/physics.F90 b/src/physics.F90 index 8fa0575cd..7988970ef 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 !=============================================================================== diff --git a/src/tally.F90 b/src/tally.F90 index 87f5a2cab..7c7747ff5 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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 diff --git a/src/tracking.F90 b/src/tracking.F90 index d108325fa..8eca3871c 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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 From bab1d485f2bd2df4fdb952046a1d8273fc73c756 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 15 Nov 2014 07:00:02 -0500 Subject: [PATCH 10/13] Fixed typo for #341 --- src/relaxng/geometry.rnc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index e7f81c485..62827f04c 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -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:intr* } })? & + (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+ } })? }* From 3706b38bf20713cd03ce0c5f1c3406dd8286968d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 18 Nov 2014 22:15:35 -0500 Subject: [PATCH 11/13] Make global_tallies allocatable instead of automatic. --- src/global.F90 | 23 ++++++++++++----------- src/tally_initialize.F90 | 3 +++ 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index f4c50b833..874b8d0d1 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -62,7 +62,7 @@ module global ! Cross section arrays type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables - type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings + type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings ! Cross section caches type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide @@ -119,7 +119,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(:) @@ -300,7 +300,7 @@ module global logical :: write_initial_source = .false. ! ============================================================================ - ! CMFD VARIABLES + ! CMFD VARIABLES ! Main object type(cmfd_type) :: cmfd @@ -310,11 +310,11 @@ module global ! CMFD communicator integer :: cmfd_comm - + ! Timing objects type(Timer) :: time_cmfd ! timer for whole cmfd calculation type(Timer) :: time_cmfdbuild ! timer for matrix build - type(Timer) :: time_cmfdsolve ! timer for solver + type(Timer) :: time_cmfdsolve ! timer for solver ! Flag for active core map logical :: cmfd_coremap = .false. @@ -390,7 +390,7 @@ module global ! RESONANCE SCATTERING VARIABLES logical :: treat_res_scat = .false. ! is resonance scattering treated? - integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers + integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, & @@ -399,14 +399,14 @@ module global contains !=============================================================================== -! FREE_MEMORY deallocates and clears all global allocatable arrays in the +! FREE_MEMORY deallocates and clears all global allocatable arrays in the ! program !=============================================================================== subroutine free_memory() - + integer :: i ! Loop Index - + ! Deallocate cells, surfaces, materials if (allocated(cells)) deallocate(cells) if (allocated(universes)) deallocate(universes) @@ -449,6 +449,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 @@ -488,7 +489,7 @@ contains ! Deallocate track_identifiers if (allocated(track_identifiers)) deallocate(track_identifiers) - + ! Deallocate dictionaries call cell_dict % clear() call universe_dict % clear() @@ -525,7 +526,7 @@ contains if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width) deallocate(ufs_mesh) end if - + end subroutine free_memory end module global diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 23ab5b563..b7dc5ed98 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -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() From de2e293c0aa17366bbdb2916156b525be171e3df Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 19 Nov 2014 17:50:58 -0500 Subject: [PATCH 12/13] Fix typo from last merge. --- 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 72fda1905..7132a67cd 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1198,7 +1198,7 @@ sub-elements: 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 lcation. For example, if + 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 From 233cdb3cd00b8a8642d4a740382808fb01500a20 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 11 Dec 2014 11:03:16 -0500 Subject: [PATCH 13/13] Mention CMake in quick install instructions. --- docs/source/quickinstall.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/source/quickinstall.rst b/docs/source/quickinstall.rst index 9dd6be6d1..3fb9036d3 100644 --- a/docs/source/quickinstall.rst +++ b/docs/source/quickinstall.rst @@ -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