mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Merge pull request #948 from paulromano/source-norm
Fix source normalization for fixed-source simulations
This commit is contained in:
Adam Nelson
GitHub
commit
6a16da00d2
32 changed files with 602 additions and 596 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -366,7 +366,7 @@ contains
|
|||
|
||||
subroutine cmfd_tally_reset()
|
||||
|
||||
use output, only: write_message
|
||||
use error, only: write_message
|
||||
|
||||
integer :: i ! loop counter
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
319
src/geometry.F90
319
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.
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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) = &
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
278
src/tracking.F90
278
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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,8 +0,0 @@
|
|||
<?xml version="1.0"?>
|
||||
<geometry>
|
||||
|
||||
<!-- Sphere with radius 10 -->
|
||||
<surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
|
||||
<cell id="1" material="1" region="-1" />
|
||||
|
||||
</geometry>
|
||||
31
tests/test_fixed_source/inputs_true.dat
Normal file
31
tests/test_fixed_source/inputs_true.dat
Normal file
|
|
@ -0,0 +1,31 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<geometry>
|
||||
<cell id="1" material="1" region="-1" universe="1" />
|
||||
<surface boundary="vacuum" coeffs="0.0 0.0 0.0 10.0" id="1" type="sphere" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
<material id="1">
|
||||
<density units="g/cc" value="7.5" />
|
||||
<nuclide ao="1.0" name="O16" />
|
||||
<nuclide ao="0.0001" name="U238" />
|
||||
</material>
|
||||
</materials>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<settings>
|
||||
<run_mode>fixed source</run_mode>
|
||||
<particles>100</particles>
|
||||
<batches>10</batches>
|
||||
<source strength="10.0">
|
||||
<space type="point">
|
||||
<parameters>0.0 0.0 0.0</parameters>
|
||||
</space>
|
||||
</source>
|
||||
<temperature_default>294</temperature_default>
|
||||
</settings>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<tallies>
|
||||
<tally id="1">
|
||||
<scores>flux</scores>
|
||||
</tally>
|
||||
</tallies>
|
||||
|
|
@ -1,10 +0,0 @@
|
|||
<?xml version="1.0"?>
|
||||
<materials>
|
||||
|
||||
<material id="1">
|
||||
<density value="7.5" units="g/cc" />
|
||||
<nuclide name="O16" ao="1.0" />
|
||||
<nuclide name="U238" ao="0.0001" />
|
||||
</material>
|
||||
|
||||
</materials>
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,14 +0,0 @@
|
|||
<?xml version="1.0"?>
|
||||
<settings>
|
||||
|
||||
<run_mode>fixed source</run_mode>
|
||||
<batches>10</batches>
|
||||
<particles>100</particles>
|
||||
|
||||
<temperature_default>294</temperature_default>
|
||||
|
||||
<source>
|
||||
<space type="point" parameters="0 0 0" />
|
||||
</source>
|
||||
|
||||
</settings>
|
||||
|
|
@ -1,8 +0,0 @@
|
|||
<?xml version="1.0"?>
|
||||
<tallies>
|
||||
|
||||
<tally id="1">
|
||||
<scores>flux</scores>
|
||||
</tally>
|
||||
|
||||
</tallies>
|
||||
|
|
@ -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()
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue