diff --git a/CMakeLists.txt b/CMakeLists.txt index 8ca62b3e58..c5150745b6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -377,7 +377,6 @@ set(LIBOPENMC_FORTRAN_SRC src/output.F90 src/particle_header.F90 src/particle_restart.F90 - src/particle_restart_write.F90 src/physics_common.F90 src/physics.F90 src/physics_mg.F90 diff --git a/src/api.F90 b/src/api.F90 index 4538bbf3ee..326b755b6a 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -112,6 +112,7 @@ contains legendre_to_tabular = .true. legendre_to_tabular_points = 33 n_batch_interval = 1 + n_lost_particles = 0 n_particles = 0 n_source_points = 0 n_state_points = 0 diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 5eb2158e91..521be5e3f2 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -644,7 +644,7 @@ contains subroutine compute_dhat() use constants, only: CMFD_NOACCEL, ZERO - use output, only: write_message + use error, only: write_message use string, only: to_str integer :: nx ! maximum number of cells in x direction diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 140166e0a4..9cdb751a40 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -366,7 +366,7 @@ contains subroutine cmfd_tally_reset() - use output, only: write_message + use error, only: write_message integer :: i ! loop counter diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index a1827febd4..dbfabb2540 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -41,8 +41,7 @@ contains subroutine read_cmfd_xml() use constants, only: ZERO, ONE - use error, only: fatal_error, warning - use output, only: write_message + use error, only: fatal_error, warning, write_message use string, only: to_lower use xml_interface use, intrinsic :: ISO_FORTRAN_ENV diff --git a/src/error.F90 b/src/error.F90 index 92902a04cb..27004d729a 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -5,12 +5,14 @@ module error use constants use message_passing + use settings, only: verbosity implicit none private public :: fatal_error public :: warning + public :: write_message ! Error codes integer(C_INT), public, bind(C) :: E_UNASSIGNED = -1 @@ -192,4 +194,59 @@ contains end subroutine fatal_error +!=============================================================================== +! WRITE_MESSAGE displays an informational message to the log file and the +! standard output stream. +!=============================================================================== + + subroutine write_message(message, level) + character(*), intent(in) :: message ! message to write + integer, intent(in), optional :: level ! verbosity level + + integer :: i_start ! starting position + integer :: i_end ! ending position + integer :: line_wrap ! length of line + integer :: length ! length of message + integer :: last_space ! index of last space (relative to start) + + ! Set length of line + line_wrap = 80 + + ! Only allow master to print to screen + if (.not. master .and. present(level)) return + + if (.not. present(level) .or. level <= verbosity) then + ! Determine length of message + length = len_trim(message) + + i_start = 0 + do + if (length - i_start < line_wrap + 1) then + ! Remainder of message will fit on line + write(OUTPUT_UNIT, fmt='(1X,A)') message(i_start+1:length) + exit + + else + ! Determine last space in current line + last_space = index(message(i_start+1:i_start+line_wrap), & + ' ', BACK=.true.) + if (last_space == 0) then + i_end = min(length + 1, i_start+line_wrap) - 1 + write(OUTPUT_UNIT, fmt='(1X,A)') message(i_start+1:i_end) + else + i_end = i_start + last_space + write(OUTPUT_UNIT, fmt='(1X,A)') message(i_start+1:i_end-1) + end if + + ! Write up to last space + + ! Advance starting position + i_start = i_end + if (i_start > length) exit + end if + end do + end if + + end subroutine write_message + end module error diff --git a/src/geometry.F90 b/src/geometry.F90 index 1cc27e09e7..ced62e6813 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -1,18 +1,14 @@ module geometry use constants - use error, only: fatal_error, warning + use error, only: fatal_error, warning, write_message use geometry_header - use output, only: write_message use particle_header, only: LocalCoord, Particle - use particle_restart_write, only: write_particle_restart use simulation_header use settings use surface_header use stl_vector, only: VectorInt use string, only: to_str - use tally, only: score_surface_current - use tally_header implicit none @@ -357,7 +353,7 @@ contains else ! Particle is outside the lattice. if (lat % outer == NO_OUTER_UNIVERSE) then - call handle_lost_particle(p, "Particle " // trim(to_str(p %id)) & + call p % mark_as_lost("Particle " // trim(to_str(p %id)) & // " is outside lattice " // trim(to_str(lat % id)) & // " but the lattice has no defined outer universe.") return @@ -380,276 +376,6 @@ contains end subroutine find_cell -!=============================================================================== -! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of -! the geometry, is reflected, or crosses into a new lattice or cell -!=============================================================================== - - subroutine cross_surface(p) - type(Particle), intent(inout) :: p - - real(8) :: u ! x-component of direction - real(8) :: v ! y-component of direction - real(8) :: w ! z-component of direction - real(8) :: norm ! "norm" of surface normal - real(8) :: d ! distance between point and plane - real(8) :: xyz(3) ! Saved global coordinate - integer :: i_surface ! index in surfaces - logical :: rotational ! if rotational periodic BC applied - logical :: found ! particle found in universe? - class(Surface), pointer :: surf - - i_surface = abs(p % surface) - surf => surfaces(i_surface)%obj - if (verbosity >= 10 .or. trace) then - call write_message(" Crossing surface " // trim(to_str(surf % id))) - end if - - if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE LEAKS OUT OF PROBLEM - - ! Kill particle - p % alive = .false. - - ! Score any surface current tallies -- note that the particle is moved - ! forward slightly so that if the mesh boundary is on the surface, it is - ! still processed - - if (active_current_tallies % size() > 0) then - ! TODO: Find a better solution to score surface currents than - ! physically moving the particle forward slightly - - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) - end if - - ! Score to global leakage tally - global_tally_leakage = global_tally_leakage + p % wgt - - ! Display message - if (verbosity >= 10 .or. trace) then - call write_message(" Leaked out of surface " & - &// trim(to_str(surf % id))) - end if - return - - elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE REFLECTS FROM SURFACE - - ! Do not handle reflective boundary conditions on lower universes - if (p % n_coord /= 1) then - call handle_lost_particle(p, "Cannot reflect particle " & - &// trim(to_str(p % id)) // " off surface in a lower universe.") - return - end if - - ! Score surface currents since reflection causes the direction of the - ! particle to change -- artificially move the particle slightly back in - ! case the surface crossing is coincident with a mesh boundary - - if (active_current_tallies % size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) - p % coord(1) % xyz = xyz - end if - - ! Reflect particle off surface - call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw) - - ! Make sure new particle direction is normalized - u = p%coord(1)%uvw(1) - v = p%coord(1)%uvw(2) - w = p%coord(1)%uvw(3) - norm = sqrt(u*u + v*v + w*w) - p%coord(1)%uvw(:) = [u, v, w] / norm - - ! Reassign particle's cell and surface - p % coord(1) % cell = p % last_cell(p % last_n_coord) - p % surface = -p % surface - - ! If a reflective surface is coincident with a lattice or universe - ! boundary, it is necessary to redetermine the particle's coordinates in - ! the lower universes. - - p % n_coord = 1 - call find_cell(p, found) - if (.not. found) then - call handle_lost_particle(p, "Couldn't find particle after reflecting& - & from surface " // trim(to_str(surf%id)) // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Reflected from surface " & - &// trim(to_str(surf%id))) - end if - return - elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then - ! ======================================================================= - ! PERIODIC BOUNDARY - - ! Do not handle periodic boundary conditions on lower universes - if (p % n_coord /= 1) then - call handle_lost_particle(p, "Cannot transfer particle " & - // trim(to_str(p % id)) // " across surface in a lower universe.& - & Boundary conditions must be applied to universe 0.") - return - end if - - ! Score surface currents since reflection causes the direction of the - ! particle to change -- artificially move the particle slightly back in - ! case the surface crossing is coincident with a mesh boundary - - if (active_current_tallies % size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_surface_current(p) - p % coord(1) % xyz = xyz - end if - - rotational = .false. - select type (surf) - type is (SurfaceXPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceXPlane) - p % coord(1) % xyz(1) = opposite % x0 - type is (SurfaceYPlane) - rotational = .true. - - ! Rotate direction - u = p % coord(1) % uvw(1) - v = p % coord(1) % uvw(2) - p % coord(1) % uvw(1) = v - p % coord(1) % uvw(2) = -u - - ! Rotate position - p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0 - p % coord(1) % xyz(2) = opposite % y0 - end select - - type is (SurfaceYPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceYPlane) - p % coord(1) % xyz(2) = opposite % y0 - type is (SurfaceXPlane) - rotational = .true. - - ! Rotate direction - u = p % coord(1) % uvw(1) - v = p % coord(1) % uvw(2) - p % coord(1) % uvw(1) = -v - p % coord(1) % uvw(2) = u - - ! Rotate position - p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0 - p % coord(1) % xyz(1) = opposite % x0 - end select - - type is (SurfaceZPlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfaceZPlane) - p % coord(1) % xyz(3) = opposite % z0 - end select - - type is (SurfacePlane) - select type (opposite => surfaces(surf % i_periodic) % obj) - type is (SurfacePlane) - ! Get surface normal for opposite plane - xyz(:) = opposite % normal(p % coord(1) % xyz) - - ! Determine distance to plane - norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3) - d = opposite % evaluate(p % coord(1) % xyz) / norm - - ! Move particle along normal vector based on distance - p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz - end select - end select - - ! Reassign particle's surface - if (rotational) then - p % surface = surf % i_periodic - else - p % surface = sign(surf % i_periodic, p % surface) - end if - - ! Figure out what cell particle is in now - p % n_coord = 1 - call find_cell(p, found) - if (.not. found) then - call handle_lost_particle(p, "Couldn't find particle after hitting & - &periodic boundary on surface " // trim(to_str(surf%id)) // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Hit periodic boundary on surface " & - // trim(to_str(surf%id))) - end if - return - end if - - ! ========================================================================== - ! SEARCH NEIGHBOR LISTS FOR NEXT CELL - - if (p % surface > 0 .and. allocated(surf%neighbor_pos)) then - ! If coming from negative side of surface, search all the neighboring - ! cells on the positive side - - call find_cell(p, found, surf%neighbor_pos) - if (found) return - - elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then - ! If coming from positive side of surface, search all the neighboring - ! cells on the negative side - - call find_cell(p, found, surf%neighbor_neg) - if (found) return - - end if - - ! ========================================================================== - ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS - - ! Remove lower coordinate levels and assignment of surface - p % surface = NONE - p % n_coord = 1 - call find_cell(p, found) - - if (run_mode /= MODE_PLOTTING .and. (.not. found)) then - ! If a cell is still not found, there are two possible causes: 1) there is - ! a void in the model, and 2) the particle hit a surface at a tangent. If - ! the particle is really traveling tangent to a surface, if we move it - ! forward a tiny bit it should fix the problem. - - p % n_coord = 1 - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call find_cell(p, found) - - ! Couldn't find next cell anywhere! This probably means there is an actual - ! undefined region in the geometry. - - if (.not. found) then - call handle_lost_particle(p, "After particle " // trim(to_str(p % id)) & - // " crossed surface " // trim(to_str(surf%id)) & - // " it could not be located in any cell and it did not leak.") - return - end if - end if - - end subroutine cross_surface - !=============================================================================== ! CROSS_LATTICE moves a particle into a new lattice element !=============================================================================== @@ -690,7 +416,7 @@ contains call find_cell(p, found) if (.not. found) then if (p % alive) then ! Particle may have been killed in find_cell - call handle_lost_particle(p, "Could not locate particle " & + call p % mark_as_lost("Could not locate particle " & // trim(to_str(p % id)) // " after crossing a lattice boundary.") return end if @@ -713,9 +439,8 @@ contains ! Search for particle call find_cell(p, found) if (.not. found) then - call handle_lost_particle(p, "Could not locate particle " & - // trim(to_str(p % id)) & - // " after crossing a lattice boundary.") + call p % mark_as_lost("Could not locate particle " // & + trim(to_str(p % id)) // " after crossing a lattice boundary.") return end if end if @@ -992,7 +717,7 @@ contains end select LAT_TYPE if (d_lat < ZERO) then - call handle_lost_particle(p, "Particle " // trim(to_str(p % id)) & + call p % mark_as_lost("Particle " // trim(to_str(p % id)) & //" had a negative distance to a lattice boundary. d = " & //trim(to_str(d_lat))) end if @@ -1093,38 +818,6 @@ contains end subroutine neighbor_lists -!=============================================================================== -! HANDLE_LOST_PARTICLE -!=============================================================================== - - subroutine handle_lost_particle(p, message) - - type(Particle), intent(inout) :: p - character(*) :: message - - integer(8) :: tot_n_particles - - ! Print warning and write lost particle file - call warning(message) - call write_particle_restart(p) - - ! Increment number of lost particles - p % alive = .false. -!$omp atomic - n_lost_particles = n_lost_particles + 1 - - ! Count the total number of simulated particles (on this processor) - tot_n_particles = n_batches * gen_per_batch * work - - ! Abort the simulation if the maximum number of lost particles has been - ! reached - if (n_lost_particles >= MAX_LOST_PARTICLES .and. & - n_lost_particles >= REL_MAX_LOST_PARTICLES * tot_n_particles) then - call fatal_error("Maximum number of lost particles has been reached.") - end if - - end subroutine handle_lost_particle - !=============================================================================== ! CALC_OFFSETS calculates and stores the offsets in all fill cells. This ! routine is called once upon initialization. diff --git a/src/initialize.F90 b/src/initialize.F90 index 21d568c72a..392e8d6bab 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -10,7 +10,7 @@ module initialize use bank_header, only: Bank use constants use set_header, only: SetInt - use error, only: fatal_error, warning + use error, only: fatal_error, warning, write_message use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& root_universe use hdf5_interface, only: file_open, read_attribute, file_close, & @@ -19,7 +19,7 @@ module initialize use material_header, only: Material use message_passing use mgxs_data, only: read_mgxs, create_macro_xs - use output, only: print_version, write_message, print_usage + use output, only: print_version, print_usage use random_lcg, only: openmc_set_seed use settings #ifdef _OPENMP diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 91cd57ae29..e95762a0b3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -10,7 +10,7 @@ module input_xml use distribution_multivariate use distribution_univariate use endf, only: reaction_name - use error, only: fatal_error, warning + use error, only: fatal_error, warning, write_message use geometry, only: calc_offsets, maximum_levels, count_instance, & neighbor_lists use geometry_header @@ -23,7 +23,7 @@ module input_xml use mgxs_header use multipole, only: multipole_read use nuclide_header - use output, only: write_message, title, header, print_plot + use output, only: title, header, print_plot use plot_header use random_lcg, only: prn, openmc_set_seed use surface_header diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index f2aaa00de8..1b14e6e9f0 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -3,13 +3,12 @@ module mgxs_data use constants use algorithm, only: find use dict_header, only: DictCharInt - use error, only: fatal_error + use error, only: fatal_error, write_message use geometry_header, only: get_temperatures, cells use hdf5_interface use material_header, only: Material, materials, n_materials use mgxs_header use nuclide_header, only: n_nuclides - use output, only: write_message use set_header, only: SetChar use settings use stl_vector, only: VectorReal diff --git a/src/output.F90 b/src/output.F90 index c200775840..d09beda1b6 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -202,61 +202,6 @@ contains end subroutine print_usage -!=============================================================================== -! WRITE_MESSAGE displays an informational message to the log file and the -! standard output stream. -!=============================================================================== - - subroutine write_message(message, level) - character(*), intent(in) :: message ! message to write - integer, intent(in), optional :: level ! verbosity level - - integer :: i_start ! starting position - integer :: i_end ! ending position - integer :: line_wrap ! length of line - integer :: length ! length of message - integer :: last_space ! index of last space (relative to start) - - ! Set length of line - line_wrap = 80 - - ! Only allow master to print to screen - if (.not. master .and. present(level)) return - - if (.not. present(level) .or. level <= verbosity) then - ! Determine length of message - length = len_trim(message) - - i_start = 0 - do - if (length - i_start < line_wrap + 1) then - ! Remainder of message will fit on line - write(ou, fmt='(1X,A)') message(i_start+1:length) - exit - - else - ! Determine last space in current line - last_space = index(message(i_start+1:i_start+line_wrap), & - ' ', BACK=.true.) - if (last_space == 0) then - i_end = min(length + 1, i_start+line_wrap) - 1 - write(ou, fmt='(1X,A)') message(i_start+1:i_end) - else - i_end = i_start + last_space - write(ou, fmt='(1X,A)') message(i_start+1:i_end-1) - end if - - ! Write up to last space - - ! Advance starting position - i_start = i_end - if (i_start > length) exit - end if - end do - end if - - end subroutine write_message - !=============================================================================== ! PRINT_PARTICLE displays the attributes of a particle !=============================================================================== diff --git a/src/particle_header.F90 b/src/particle_header.F90 index bda7f03b43..552f9b042d 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -1,20 +1,27 @@ module particle_header - use bank_header, only: Bank - use constants, only: NEUTRON, ONE, NONE, ZERO, MAX_SECONDARY, & - MAX_DELAYED_GROUPS, ERROR_REAL - use error, only: fatal_error + use hdf5, only: HID_T + + use bank_header, only: Bank, source_bank + use constants + use error, only: fatal_error, warning use geometry_header, only: root_universe + use hdf5_interface + use settings + use simulation_header + use string, only: to_str implicit none + private + !=============================================================================== ! LOCALCOORD describes the location of a particle local to a single ! universe. When the geometry consists of nested universes, a particle will have ! a list of coordinates in each level !=============================================================================== - type LocalCoord + type, public :: LocalCoord ! Indices in various arrays for this level integer :: cell = NONE @@ -39,7 +46,7 @@ module particle_header ! geometry !=============================================================================== - type Particle + type, public :: Particle ! Basic data integer(8) :: id ! Unique ID integer :: type ! Particle type (n, p, e, etc) @@ -107,20 +114,84 @@ module particle_header type(Bank) :: secondary_bank(MAX_SECONDARY) contains - procedure :: initialize => initialize_particle - procedure :: clear => clear_particle - procedure :: initialize_from_source + procedure :: clear procedure :: create_secondary + procedure :: initialize + procedure :: initialize_from_source + procedure :: mark_as_lost + procedure :: write_restart end type Particle contains !=============================================================================== -! INITIALIZE_PARTICLE sets default attributes for a particle from the source -! bank +! RESET_COORD clears data from a single coordinate level !=============================================================================== - subroutine initialize_particle(this) + elemental subroutine reset_coord(this) + class(LocalCoord), intent(inout) :: this + + this % cell = NONE + this % universe = NONE + this % lattice = NONE + this % lattice_x = NONE + this % lattice_y = NONE + this % lattice_z = NONE + this % rotated = .false. + + end subroutine reset_coord + +!=============================================================================== +! CLEAR_PARTICLE resets all coordinate levels for the particle +!=============================================================================== + + subroutine clear(this) + class(Particle) :: this + + integer :: i + + ! remove any coordinate levels + do i = 1, MAX_COORD + call this % coord(i) % reset() + end do + end subroutine clear + +!=============================================================================== +! CREATE_SECONDARY stores the current phase space attributes of the particle in +! the secondary bank and increments the number of sites in the secondary bank. +!=============================================================================== + + subroutine create_secondary(this, uvw, type, run_CE) + class(Particle), intent(inout) :: this + real(8), intent(in) :: uvw(3) + integer, intent(in) :: type + logical, intent(in) :: run_CE + + integer(8) :: n + + ! Check to make sure that the hard-limit on secondary particles is not + ! exceeded. + if (this % n_secondary == MAX_SECONDARY) then + call fatal_error("Too many secondary particles created.") + end if + + n = this % n_secondary + 1 + this % secondary_bank(n) % wgt = this % wgt + this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz + this % secondary_bank(n) % uvw(:) = uvw + this % n_secondary = n + this % secondary_bank(this % n_secondary) % E = this % E + if (.not. run_CE) then + this % secondary_bank(this % n_secondary) % E = real(this % g, 8) + end if + + end subroutine create_secondary + +!=============================================================================== +! INITIALIZE sets default attributes for a particle from the source bank +!=============================================================================== + + subroutine initialize(this) class(Particle) :: this @@ -154,40 +225,7 @@ contains this % n_coord = 1 this % last_n_coord = 1 - end subroutine initialize_particle - -!=============================================================================== -! CLEAR_PARTICLE resets all coordinate levels for the particle -!=============================================================================== - - subroutine clear_particle(this) - - class(Particle) :: this - integer :: i - - ! remove any coordinate levels - do i = 1, MAX_COORD - call this % coord(i) % reset() - end do - - end subroutine clear_particle - -!=============================================================================== -! RESET_COORD clears data from a single coordinate level -!=============================================================================== - - elemental subroutine reset_coord(this) - class(LocalCoord), intent(inout) :: this - - this % cell = NONE - this % universe = NONE - this % lattice = NONE - this % lattice_x = NONE - this % lattice_y = NONE - this % lattice_z = NONE - this % rotated = .false. - - end subroutine reset_coord + end subroutine initialize !=============================================================================== ! INITIALIZE_FROM_SOURCE initializes a particle from data stored in a source @@ -225,34 +263,90 @@ contains end subroutine initialize_from_source !=============================================================================== -! CREATE_SECONDARY stores the current phase space attributes of the particle in -! the secondary bank and increments the number of sites in the secondary bank. +! MARK_AS_LOST !=============================================================================== - subroutine create_secondary(this, uvw, type, run_CE) + subroutine mark_as_lost(this, message) class(Particle), intent(inout) :: this - real(8), intent(in) :: uvw(3) - integer, intent(in) :: type - logical, intent(in) :: run_CE + character(*) :: message - integer(8) :: n + integer(8) :: tot_n_particles - ! Check to make sure that the hard-limit on secondary particles is not - ! exceeded. - if (this % n_secondary == MAX_SECONDARY) then - call fatal_error("Too many secondary particles created.") + ! Print warning and write lost particle file + call warning(message) + call this % write_restart() + + ! Increment number of lost particles + this % alive = .false. +!$omp atomic + n_lost_particles = n_lost_particles + 1 + + ! Count the total number of simulated particles (on this processor) + tot_n_particles = current_batch * gen_per_batch * work + + ! Abort the simulation if the maximum number of lost particles has been + ! reached + if (n_lost_particles >= MAX_LOST_PARTICLES .and. & + n_lost_particles >= REL_MAX_LOST_PARTICLES * tot_n_particles) then + call fatal_error("Maximum number of lost particles has been reached.") end if - n = this % n_secondary + 1 - this % secondary_bank(n) % wgt = this % wgt - this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz - this % secondary_bank(n) % uvw(:) = uvw - this % n_secondary = n - this % secondary_bank(this % n_secondary) % E = this % E - if (.not. run_CE) then - this % secondary_bank(this % n_secondary) % E = real(this % g, 8) - end if + end subroutine mark_as_lost - end subroutine create_secondary +!=============================================================================== +! WRITE_RESTART creates a particle restart file +!=============================================================================== + + subroutine write_restart(this) + class(Particle), intent(in) :: this + + integer(HID_T) :: file_id + character(MAX_FILE_LEN) :: filename + + ! Dont write another restart file if in particle restart mode + if (run_mode == MODE_PARTICLE) return + + ! Set up file name + filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) & + // '_' // trim(to_str(this % id)) // '.h5' + +!$omp critical (WriteParticleRestart) + ! Create file + file_id = file_create(filename) + + associate (src => source_bank(current_work)) + ! Write filetype and version info + call write_attribute(file_id, 'filetype', 'particle restart') + call write_attribute(file_id, 'version', VERSION_PARTICLE_RESTART) + call write_attribute(file_id, "openmc_version", VERSION) +#ifdef GIT_SHA1 + call write_attribute(file_id, "git_sha1", GIT_SHA1) +#endif + + ! Write data to file + call write_dataset(file_id, 'current_batch', current_batch) + call write_dataset(file_id, 'generations_per_batch', gen_per_batch) + call write_dataset(file_id, 'current_generation', current_gen) + call write_dataset(file_id, 'n_particles', n_particles) + select case(run_mode) + case (MODE_FIXEDSOURCE) + call write_dataset(file_id, 'run_mode', 'fixed source') + case (MODE_EIGENVALUE) + call write_dataset(file_id, 'run_mode', 'eigenvalue') + case (MODE_PARTICLE) + call write_dataset(file_id, 'run_mode', 'particle restart') + end select + call write_dataset(file_id, 'id', this % id) + call write_dataset(file_id, 'weight', src % wgt) + call write_dataset(file_id, 'energy', src % E) + call write_dataset(file_id, 'xyz', src % xyz) + call write_dataset(file_id, 'uvw', src % uvw) + end associate + + ! Close file + call file_close(file_id) +!$omp end critical (WriteParticleRestart) + + end subroutine write_restart end module particle_header diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 53332c794f..2d15f2ddc4 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -4,10 +4,11 @@ module particle_restart use bank_header, only: Bank use constants + use error, only: write_message use hdf5_interface, only: file_open, file_close, read_dataset use mgxs_header, only: energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides - use output, only: write_message, print_particle + use output, only: print_particle use particle_header, only: Particle use random_lcg, only: set_particle_seed use settings diff --git a/src/particle_restart_write.F90 b/src/particle_restart_write.F90 deleted file mode 100644 index 5addf6056f..0000000000 --- a/src/particle_restart_write.F90 +++ /dev/null @@ -1,74 +0,0 @@ -module particle_restart_write - - use bank_header, only: Bank, source_bank - use hdf5_interface - use particle_header, only: Particle - use settings - use simulation_header - use string, only: to_str - - use hdf5 - - implicit none - private - public :: write_particle_restart - -contains - -!=============================================================================== -! WRITE_PARTICLE_RESTART is the main routine that writes out the particle file -!=============================================================================== - - subroutine write_particle_restart(p) - type(Particle), intent(in) :: p - - integer(HID_T) :: file_id - character(MAX_FILE_LEN) :: filename - - ! Dont write another restart file if in particle restart mode - if (run_mode == MODE_PARTICLE) return - - ! Set up file name - filename = trim(path_output) // 'particle_' // trim(to_str(current_batch)) & - // '_' // trim(to_str(p%id)) // '.h5' - -!$omp critical (WriteParticleRestart) - ! Create file - file_id = file_create(filename) - - associate (src => source_bank(current_work)) - ! Write filetype and version info - call write_attribute(file_id, 'filetype', 'particle restart') - call write_attribute(file_id, 'version', VERSION_PARTICLE_RESTART) - call write_attribute(file_id, "openmc_version", VERSION) -#ifdef GIT_SHA1 - call write_attribute(file_id, "git_sha1", GIT_SHA1) -#endif - - ! Write data to file - call write_dataset(file_id, 'current_batch', current_batch) - call write_dataset(file_id, 'generations_per_batch', gen_per_batch) - call write_dataset(file_id, 'current_generation', current_gen) - call write_dataset(file_id, 'n_particles', n_particles) - select case(run_mode) - case (MODE_FIXEDSOURCE) - call write_dataset(file_id, 'run_mode', 'fixed source') - case (MODE_EIGENVALUE) - call write_dataset(file_id, 'run_mode', 'eigenvalue') - case (MODE_PARTICLE) - call write_dataset(file_id, 'run_mode', 'particle restart') - end select - call write_dataset(file_id, 'id', p%id) - call write_dataset(file_id, 'weight', src%wgt) - call write_dataset(file_id, 'energy', src%E) - call write_dataset(file_id, 'xyz', src%xyz) - call write_dataset(file_id, 'uvw', src%uvw) - end associate - - ! Close file - call file_close(file_id) -!$omp end critical (WriteParticleRestart) - - end subroutine write_particle_restart - -end module particle_restart_write diff --git a/src/physics.F90 b/src/physics.F90 index fd6a0087c9..d13313e05a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -4,15 +4,13 @@ module physics use constants use cross_section, only: elastic_xs_0K use endf, only: reaction_name - use error, only: fatal_error, warning + use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math use mesh_header, only: meshes use message_passing use nuclide_header - use output, only: write_message use particle_header, only: Particle - use particle_restart_write, only: write_particle_restart use physics_common use random_lcg, only: prn, advance_prn_seed, prn_set_stream use reaction_header, only: Reaction @@ -172,7 +170,7 @@ contains ! Check to make sure that a nuclide was sampled if (i_nuc_mat > mat % n_nuclides) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Did not sample any nuclide during collision.") end if @@ -384,7 +382,7 @@ contains ! Check to make sure inelastic scattering reaction sampled if (i > size(nuc % reactions)) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Did not sample any reaction for nuclide " & &// trim(nuc % name)) end if @@ -1098,7 +1096,7 @@ contains ! Determine indices on ufs mesh for current location call m % get_bin(p % coord(1) % xyz, mesh_bin) if (mesh_bin == NO_BIN_FOUND) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Source site outside UFS mesh!") end if @@ -1251,7 +1249,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) + ! call p % write_restart() call fatal_error("Resampled energy distribution maximum number of " & // "times for nuclide " // nuc % name) end if @@ -1275,7 +1273,7 @@ contains ! check for large number of resamples n_sample = n_sample + 1 if (n_sample == MAX_SAMPLE) then - ! call write_particle_restart(p) + ! call p % write_restart() call fatal_error("Resampled energy distribution maximum number of " & // "times for nuclide " // nuc % name) end if diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 42094f6d06..eb82085941 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -4,16 +4,14 @@ module physics_mg use bank_header use constants - use error, only: fatal_error, warning + use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math, only: rotate_angle use mesh_header, only: meshes use mgxs_header use message_passing use nuclide_header, only: material_xs - use output, only: write_message use particle_header, only: Particle - use particle_restart_write, only: write_particle_restart use physics_common use random_lcg, only: prn use scattdata_header @@ -196,7 +194,7 @@ contains call m % get_bin(p % coord(1) % xyz, mesh_bin) if (mesh_bin == NO_BIN_FOUND) then - call write_particle_restart(p) + call p % write_restart() call fatal_error("Source site outside UFS mesh!") end if diff --git a/src/plot.F90 b/src/plot.F90 index fe2ddd13b7..5eb894ce8d 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -5,11 +5,11 @@ module plot use hdf5 use constants - use error, only: fatal_error + use error, only: fatal_error, write_message use geometry, only: find_cell, check_cell_overlap use geometry_header, only: Cell, root_universe, cells use hdf5_interface - use output, only: write_message, time_stamp + use output, only: time_stamp use material_header, only: materials use particle_header, only: LocalCoord, Particle use plot_header diff --git a/src/simulation.F90 b/src/simulation.F90 index 5eacbd40b8..b04e967f27 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -16,12 +16,12 @@ module simulation #ifdef _OPENMP use eigenvalue, only: join_bank_from_threads #endif - use error, only: fatal_error + use error, only: fatal_error, write_message use geometry_header, only: n_cells use message_passing use mgxs_header, only: energy_bins, energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides - use output, only: write_message, header, print_columns, & + use output, only: header, print_columns, & print_batch_keff, print_generation, print_runtime, & print_results, print_overlap_check, write_tallies use particle_header, only: Particle diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 5041f0793c..825b88129f 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -13,7 +13,7 @@ module simulation_header ! GEOMETRY-RELATED VARIABLES ! Number of lost particles - integer :: n_lost_particles + integer :: n_lost_particles = 0 real(8) :: log_spacing ! spacing on logarithmic grid diff --git a/src/state_point.F90 b/src/state_point.F90 index 6d647b72e6..f7a8928987 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -19,13 +19,13 @@ module state_point use constants use eigenvalue, only: openmc_get_keff use endf, only: reaction_name - use error, only: fatal_error, warning + use error, only: fatal_error, warning, write_message use hdf5_interface use mesh_header, only: RegularMesh, meshes, n_meshes use message_passing use mgxs_header, only: nuclides_MG use nuclide_header, only: nuclides - use output, only: write_message, time_stamp + use output, only: time_stamp use random_lcg, only: seed use settings use simulation_header diff --git a/src/summary.F90 b/src/summary.F90 index 248b48e865..0f45ef1bbd 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -4,6 +4,7 @@ module summary use constants use endf, only: reaction_name + use error, only: write_message use geometry_header use hdf5_interface use material_header, only: Material, n_materials @@ -11,7 +12,7 @@ module summary use message_passing use mgxs_header, only: nuclides_MG use nuclide_header - use output, only: time_stamp, write_message + use output, only: time_stamp use settings, only: run_CE use surface_header use string, only: to_str diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 1778474e04..21b4c8a586 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -9,7 +9,8 @@ module tally_header use dict_header, only: DictIntInt use message_passing, only: n_procs use nuclide_header, only: nuclide_dict - use settings, only: reduce_tallies + use settings, only: reduce_tallies, run_mode + use source_header, only: external_source use stl_vector, only: VectorInt use string, only: to_lower, to_f_string, str_to_int, to_str use tally_filter_header, only: TallyFilterContainer, filters, n_filters @@ -161,6 +162,7 @@ contains integer :: i, j real(C_DOUBLE) :: val + real(C_DOUBLE) :: total_source ! Increment number of realizations if (reduce_tallies) then @@ -169,10 +171,17 @@ contains this % n_realizations = this % n_realizations + n_procs end if + ! Calculate total source strength for normalization + if (run_mode == MODE_FIXEDSOURCE) then + total_source = sum(external_source(:) % strength) + else + total_source = ONE + end if + ! Accumulate each result do j = 1, size(this % results, 3) do i = 1, size(this % results, 2) - val = this % results(RESULT_VALUE, i, j)/total_weight + val = this % results(RESULT_VALUE, i, j)/total_weight * total_source this % results(RESULT_VALUE, i, j) = ZERO this % results(RESULT_SUM, i, j) = & diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 772e739bad..9e9b06d205 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -8,8 +8,8 @@ module trigger use constants use eigenvalue, only: openmc_get_keff + use error, only: warning, write_message use string, only: to_str - use output, only: warning, write_message use mesh_header, only: RegularMesh, meshes use message_passing, only: master use settings diff --git a/src/tracking.F90 b/src/tracking.F90 index 726a263adf..65769a7f18 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -2,11 +2,10 @@ module tracking use constants use cross_section, only: calculate_xs - use error, only: fatal_error, warning + use error, only: fatal_error, warning, write_message use geometry_header, only: cells - use geometry, only: find_cell, distance_to_boundary, cross_surface, & - cross_lattice, check_cell_overlap - use output, only: write_message + use geometry, only: find_cell, distance_to_boundary, cross_lattice, & + check_cell_overlap use message_passing use mgxs_header use nuclide_header @@ -17,6 +16,7 @@ module tracking use settings use simulation_header use string, only: to_str + use surface_header use tally_header use tally, only: score_analog_tally, score_tracklength_tally, & score_collision_tally, score_surface_current, & @@ -278,4 +278,274 @@ contains end subroutine transport +!=============================================================================== +! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of +! the geometry, is reflected, or crosses into a new lattice or cell +!=============================================================================== + + subroutine cross_surface(p) + type(Particle), intent(inout) :: p + + real(8) :: u ! x-component of direction + real(8) :: v ! y-component of direction + real(8) :: w ! z-component of direction + real(8) :: norm ! "norm" of surface normal + real(8) :: d ! distance between point and plane + real(8) :: xyz(3) ! Saved global coordinate + integer :: i_surface ! index in surfaces + logical :: rotational ! if rotational periodic BC applied + logical :: found ! particle found in universe? + class(Surface), pointer :: surf + + i_surface = abs(p % surface) + surf => surfaces(i_surface)%obj + if (verbosity >= 10 .or. trace) then + call write_message(" Crossing surface " // trim(to_str(surf % id))) + end if + + if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then + ! ======================================================================= + ! PARTICLE LEAKS OUT OF PROBLEM + + ! Kill particle + p % alive = .false. + + ! Score any surface current tallies -- note that the particle is moved + ! forward slightly so that if the mesh boundary is on the surface, it is + ! still processed + + if (active_current_tallies % size() > 0) then + ! TODO: Find a better solution to score surface currents than + ! physically moving the particle forward slightly + + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + call score_surface_current(p) + end if + + ! Score to global leakage tally + global_tally_leakage = global_tally_leakage + p % wgt + + ! Display message + if (verbosity >= 10 .or. trace) then + call write_message(" Leaked out of surface " & + &// trim(to_str(surf % id))) + end if + return + + elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then + ! ======================================================================= + ! PARTICLE REFLECTS FROM SURFACE + + ! Do not handle reflective boundary conditions on lower universes + if (p % n_coord /= 1) then + call p % mark_as_lost("Cannot reflect particle " & + // trim(to_str(p % id)) // " off surface in a lower universe.") + return + end if + + ! Score surface currents since reflection causes the direction of the + ! particle to change -- artificially move the particle slightly back in + ! case the surface crossing is coincident with a mesh boundary + + if (active_current_tallies % size() > 0) then + xyz = p % coord(1) % xyz + p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw + call score_surface_current(p) + p % coord(1) % xyz = xyz + end if + + ! Reflect particle off surface + call surf%reflect(p%coord(1)%xyz, p%coord(1)%uvw) + + ! Make sure new particle direction is normalized + u = p%coord(1)%uvw(1) + v = p%coord(1)%uvw(2) + w = p%coord(1)%uvw(3) + norm = sqrt(u*u + v*v + w*w) + p%coord(1)%uvw(:) = [u, v, w] / norm + + ! Reassign particle's cell and surface + p % coord(1) % cell = p % last_cell(p % last_n_coord) + p % surface = -p % surface + + ! If a reflective surface is coincident with a lattice or universe + ! boundary, it is necessary to redetermine the particle's coordinates in + ! the lower universes. + + p % n_coord = 1 + call find_cell(p, found) + if (.not. found) then + call p % mark_as_lost("Couldn't find particle after reflecting& + & from surface " // trim(to_str(surf % id)) // ".") + return + end if + + ! Set previous coordinate going slightly past surface crossing + p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + + ! Diagnostic message + if (verbosity >= 10 .or. trace) then + call write_message(" Reflected from surface " & + &// trim(to_str(surf%id))) + end if + return + elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then + ! ======================================================================= + ! PERIODIC BOUNDARY + + ! Do not handle periodic boundary conditions on lower universes + if (p % n_coord /= 1) then + call p % mark_as_lost("Cannot transfer particle " & + // trim(to_str(p % id)) // " across surface in a lower universe.& + & Boundary conditions must be applied to universe 0.") + return + end if + + ! Score surface currents since reflection causes the direction of the + ! particle to change -- artificially move the particle slightly back in + ! case the surface crossing is coincident with a mesh boundary + + if (active_current_tallies % size() > 0) then + xyz = p % coord(1) % xyz + p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw + call score_surface_current(p) + p % coord(1) % xyz = xyz + end if + + rotational = .false. + select type (surf) + type is (SurfaceXPlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfaceXPlane) + p % coord(1) % xyz(1) = opposite % x0 + type is (SurfaceYPlane) + rotational = .true. + + ! Rotate direction + u = p % coord(1) % uvw(1) + v = p % coord(1) % uvw(2) + p % coord(1) % uvw(1) = v + p % coord(1) % uvw(2) = -u + + ! Rotate position + p % coord(1) % xyz(1) = surf % x0 + p % coord(1) % xyz(2) - opposite % y0 + p % coord(1) % xyz(2) = opposite % y0 + end select + + type is (SurfaceYPlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfaceYPlane) + p % coord(1) % xyz(2) = opposite % y0 + type is (SurfaceXPlane) + rotational = .true. + + ! Rotate direction + u = p % coord(1) % uvw(1) + v = p % coord(1) % uvw(2) + p % coord(1) % uvw(1) = -v + p % coord(1) % uvw(2) = u + + ! Rotate position + p % coord(1) % xyz(2) = surf % y0 + p % coord(1) % xyz(1) - opposite % x0 + p % coord(1) % xyz(1) = opposite % x0 + end select + + type is (SurfaceZPlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfaceZPlane) + p % coord(1) % xyz(3) = opposite % z0 + end select + + type is (SurfacePlane) + select type (opposite => surfaces(surf % i_periodic) % obj) + type is (SurfacePlane) + ! Get surface normal for opposite plane + xyz(:) = opposite % normal(p % coord(1) % xyz) + + ! Determine distance to plane + norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3) + d = opposite % evaluate(p % coord(1) % xyz) / norm + + ! Move particle along normal vector based on distance + p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz + end select + end select + + ! Reassign particle's surface + if (rotational) then + p % surface = surf % i_periodic + else + p % surface = sign(surf % i_periodic, p % surface) + end if + + ! Figure out what cell particle is in now + p % n_coord = 1 + call find_cell(p, found) + if (.not. found) then + call p % mark_as_lost("Couldn't find particle after hitting & + &periodic boundary on surface " // trim(to_str(surf % id)) // ".") + return + end if + + ! Set previous coordinate going slightly past surface crossing + p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + + ! Diagnostic message + if (verbosity >= 10 .or. trace) then + call write_message(" Hit periodic boundary on surface " & + // trim(to_str(surf%id))) + end if + return + end if + + ! ========================================================================== + ! SEARCH NEIGHBOR LISTS FOR NEXT CELL + + if (p % surface > 0 .and. allocated(surf%neighbor_pos)) then + ! If coming from negative side of surface, search all the neighboring + ! cells on the positive side + + call find_cell(p, found, surf%neighbor_pos) + if (found) return + + elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then + ! If coming from positive side of surface, search all the neighboring + ! cells on the negative side + + call find_cell(p, found, surf%neighbor_neg) + if (found) return + + end if + + ! ========================================================================== + ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS + + ! Remove lower coordinate levels and assignment of surface + p % surface = NONE + p % n_coord = 1 + call find_cell(p, found) + + if (run_mode /= MODE_PLOTTING .and. (.not. found)) then + ! If a cell is still not found, there are two possible causes: 1) there is + ! a void in the model, and 2) the particle hit a surface at a tangent. If + ! the particle is really traveling tangent to a surface, if we move it + ! forward a tiny bit it should fix the problem. + + p % n_coord = 1 + p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw + call find_cell(p, found) + + ! Couldn't find next cell anywhere! This probably means there is an actual + ! undefined region in the geometry. + + if (.not. found) then + call p % mark_as_lost("After particle " // trim(to_str(p % id)) & + // " crossed surface " // trim(to_str(surf % id)) & + // " it could not be located in any cell and it did not leak.") + return + end if + end if + + end subroutine cross_surface + end module tracking diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 48f93f0e6a..5985a236b3 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -8,11 +8,12 @@ module volume_calc #endif use constants + use error, only: write_message use geometry, only: find_cell use geometry_header, only: universes, cells use hdf5_interface, only: file_create, file_close, write_attribute, & create_group, close_group, write_dataset - use output, only: write_message, header, time_stamp + use output, only: header, time_stamp use material_header, only: materials use message_passing use nuclide_header, only: nuclides diff --git a/tests/test_fixed_source/geometry.xml b/tests/test_fixed_source/geometry.xml deleted file mode 100644 index bc56030e18..0000000000 --- a/tests/test_fixed_source/geometry.xml +++ /dev/null @@ -1,8 +0,0 @@ - - - - - - - - diff --git a/tests/test_fixed_source/inputs_true.dat b/tests/test_fixed_source/inputs_true.dat new file mode 100644 index 0000000000..2e0d57b0c1 --- /dev/null +++ b/tests/test_fixed_source/inputs_true.dat @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + fixed source + 100 + 10 + + + 0.0 0.0 0.0 + + + 294 + + + + + flux + + diff --git a/tests/test_fixed_source/materials.xml b/tests/test_fixed_source/materials.xml deleted file mode 100644 index 6e4249da32..0000000000 --- a/tests/test_fixed_source/materials.xml +++ /dev/null @@ -1,10 +0,0 @@ - - - - - - - - - - diff --git a/tests/test_fixed_source/results_true.dat b/tests/test_fixed_source/results_true.dat index 6246dfaf62..c911b4fd25 100644 --- a/tests/test_fixed_source/results_true.dat +++ b/tests/test_fixed_source/results_true.dat @@ -1,6 +1,6 @@ tally 1: -4.448476E+02 -1.984373E+04 +4.448476E+03 +1.984373E+06 leakage: 9.790000E+00 9.588300E+00 diff --git a/tests/test_fixed_source/settings.xml b/tests/test_fixed_source/settings.xml deleted file mode 100644 index e7e7f2f5b9..0000000000 --- a/tests/test_fixed_source/settings.xml +++ /dev/null @@ -1,14 +0,0 @@ - - - - fixed source - 10 - 100 - - 294 - - - - - - diff --git a/tests/test_fixed_source/tallies.xml b/tests/test_fixed_source/tallies.xml deleted file mode 100644 index 87e08d6c73..0000000000 --- a/tests/test_fixed_source/tallies.xml +++ /dev/null @@ -1,8 +0,0 @@ - - - - - flux - - - diff --git a/tests/test_fixed_source/test_fixed_source.py b/tests/test_fixed_source/test_fixed_source.py index 8d61a46d0a..db58e2f706 100644 --- a/tests/test_fixed_source/test_fixed_source.py +++ b/tests/test_fixed_source/test_fixed_source.py @@ -5,17 +5,18 @@ import os import sys import numpy as np sys.path.insert(0, os.pardir) -from testing_harness import TestHarness -from openmc import StatePoint +from testing_harness import PyAPITestHarness +import openmc +import openmc.stats -class FixedSourceTestHarness(TestHarness): +class FixedSourceTestHarness(PyAPITestHarness): def _get_results(self): """Digest info in the statepoint and return as a string.""" # Read the statepoint file. statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] outstr = '' - with StatePoint(statepoint) as sp: + with openmc.StatePoint(statepoint) as sp: # Write out tally data. for i, tally_ind in enumerate(sp.tallies): tally = sp.tallies[tally_ind] @@ -36,5 +37,28 @@ class FixedSourceTestHarness(TestHarness): if __name__ == '__main__': - harness = FixedSourceTestHarness('statepoint.10.h5') + mat = openmc.Material() + mat.add_nuclide('O16', 1.0) + mat.add_nuclide('U238', 0.0001) + mat.set_density('g/cc', 7.5) + + surf = openmc.Sphere(R=10.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-surf) + + model = openmc.model.Model() + model.geometry.root_universe = openmc.Universe(cells=[cell]) + model.materials.append(mat) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 10 + model.settings.particles = 100 + model.settings.temperature = {'default': 294} + model.settings.source = openmc.Source(space=openmc.stats.Point(), + strength=10.0) + + tally = openmc.Tally() + tally.scores = ['flux'] + model.tallies.append(tally) + + harness = FixedSourceTestHarness('statepoint.10.h5', model) harness.main()