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 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 diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index c25a9f146..2278e11ff 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -759,7 +759,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 @@ -1172,7 +1175,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: @@ -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. diff --git a/src/geometry.F90 b/src/geometry.F90 index 043ded05b..c79f0e611 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 5db2632f4..624511912 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 10b178e0a..2cf5b7eb9 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 39b9fd887..6402463d5 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 @@ -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")) & diff --git a/src/output.F90 b/src/output.F90 index 835a7d1b0..8deee0e5c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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" diff --git a/src/physics.F90 b/src/physics.F90 index 8fa0575cd..3439a8390 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 !=============================================================================== @@ -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 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 diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index b7975d05f..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: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+ } }) & 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+ } } | 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/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() 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 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()