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

implemented half of Sterling s corrections (inheriting cellfilter, moved partial current scoring)+ bug fix on energy filter

Browse source
This commit is contained in:
Guillaume 2017-06-19 20:07:24 -04:00
parent 81d122657a
commit e584ade684
8 changed files with 95 additions and 113 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/constants.F90
View file

@ -305,7 +305,7 @@ module constants
EVENT_ABSORB = 2
! Tally score type
integer, parameter :: N_SCORE_TYPES = 26
integer, parameter :: N_SCORE_TYPES = 25
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
@ -331,8 +331,7 @@ module constants
SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value
SCORE_DECAY_RATE = -24, & ! delayed neutron precursor decay rate
SCORE_CELL_TO_CELL = -25, & ! cell to cell partial current
SCORE_CELL_TO_CELL_NORMAL_PROJECTION = -26 ! cell to cell partial current projected on surface normal
SCORE_CELL_TO_CELL = -25 ! cell to cell partial current
! Maximum scattering order supported
integer, parameter :: MAX_ANG_ORDER = 10

2
src/endf.F90
View file

@ -68,8 +68,6 @@ contains
string = "fission-q-recoverable"
case (SCORE_CELL_TO_CELL)
string = "cell-to-cell-partial-current"
case (SCORE_CELL_TO_CELL_NORMAL_PROJECTION)
string = "cell-to-cell-partial-current-projection-on-normal"
! Normal ENDF-based reactions
case (TOTAL_XS)

23
src/geometry.F90
View file

@ -611,33 +611,13 @@ contains
! cells on the positive side
call find_cell(p, found, surf%neighbor_pos)
! Save cosine of angle between surface normal and particle direction
p % normal_proj = dot_product(p % coord(1) % uvw, &
& surf % normal(p % coord(1) % xyz)) / &
& sqrt(dot_product(surf % normal(p % coord(1) % xyz)&
&, surf % normal(p % coord(1) % xyz)))
! Score cell to cell partial currents
call score_partial_current(p)
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
! Save cosine of angle between surface normal and particle direction
p % normal_proj = dot_product(p % coord(1) % uvw, &
& surf % normal(p % coord(1) % xyz)) / &
& sqrt(dot_product(surf % normal(p % coord(1) % xyz)&
&, surf % normal(p % coord(1) % xyz)))
call find_cell(p, found, surf%neighbor_neg)
! Score cell to cell partial currents
call score_partial_current(p)
if (found) return
end if
@ -671,9 +651,6 @@ contains
end if
end if
! Score cell to cell partial currents
call score_partial_current(p)
end subroutine cross_surface
!===============================================================================

3
src/input_xml.F90
View file

@ -3923,9 +3923,6 @@ contains
case ('partial_current')
t % type = TALLY_CELL_TO_CELL
t % score_bins(j) = SCORE_CELL_TO_CELL
case ('projected_partial_current')
t % type = TALLY_CELL_TO_CELL
t % score_bins(j) = SCORE_CELL_TO_CELL_NORMAL_PROJECTION
case ('events')
t % score_bins(j) = SCORE_EVENTS
case ('elastic', '(n,elastic)')

2
src/output.F90
View file

@ -734,8 +734,6 @@ contains
score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power"
score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power"
score_names(abs(SCORE_CELL_TO_CELL)) = "Partial current"
score_names(abs(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)) = "Partial &
&current projected on normal of surface"
! Create filename for tally output
filename = trim(path_output) // "tallies.out"

35
src/tally.F90
View file

@ -119,6 +119,7 @@ contains
select case(score_bin)
case (SCORE_FLUX, SCORE_FLUX_YN)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events score to a flux bin. We actually use a collision
@ -1176,12 +1177,6 @@ contains
end if
end if
case(SCORE_CELL_TO_CELL)
score = p % last_wgt * flux
case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)
score = p % last_wgt * flux * p % normal_proj
case default
if (t % estimator == ESTIMATOR_ANALOG) then
! Any other score is assumed to be a MT number. Thus, we just need
@ -2085,12 +2080,6 @@ contains
! Simply count number of scoring events
score = ONE
case(SCORE_CELL_TO_CELL)
score = p % last_wgt * flux
case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)
score = p % last_wgt * flux * p % normal_proj
end select
!#########################################################################
@ -3161,11 +3150,15 @@ contains
integer :: i_tally
integer :: i_filt
integer :: i_bin
integer :: q ! loop index for scoring bins
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
integer :: score_index ! scoring bin index
integer :: j ! loop index for scoring bins
integer :: filter_index ! single index for single bin
integer :: matching_bin ! next valid filter bin
real(8) :: flux ! collision estimate of flux
real(8) :: filter_weight ! combined weight of all filters
real(8) :: score ! analog tally score
logical :: finished ! found all valid bin combinations
type(TallyObject), pointer :: t
@ -3226,8 +3219,22 @@ contains
data(i_bin)
end do
call score_general(p, t, 0, filter_index, &
0, 0d0, flux * filter_weight)
! Determine score
score = flux * filter_weight
! Currently only one score type
SCORE_LOOP: do q = 1, t % n_user_score_bins
! determine what type of score bin
score_bin = t % score_bins(q)
! determine scoring bin index, no offset from nuclide bins
score_index = q
! Expand score if necessary and add to tally results.
call expand_and_score(p, t, score_index, filter_index, score_bin, &
score, q)
end do SCORE_LOOP
! ======================================================================
! Filter logic

128
src/tally_filter.F90
View file

@ -87,14 +87,12 @@ module tally_filter
!===============================================================================
! CELLTOFILTER specifies which geometric cells particles enter.
!===============================================================================
type, extends(TallyFilter) :: CellToFilter
integer, allocatable :: cells(:)
type(DictIntInt) :: map
type, extends(CellFilter) :: CellToFilter
! integer, allocatable :: cells(:)
! type(DictIntInt) :: map
contains
procedure :: get_next_bin => get_next_bin_cell_to
procedure :: to_statepoint => to_statepoint_cell_to
procedure :: text_label => text_label_cell_to
procedure :: initialize => initialize_cell_to
procedure :: to_statepoint_cell => to_statepoint_cell_to
procedure :: text_label_cell => text_label_cell_to
end type CellToFilter
!===============================================================================
@ -779,7 +777,7 @@ contains
integer :: i
integer, allocatable :: cell_ids(:)
call write_dataset(filter_group, "type", "cell")
call write_dataset(filter_group, "type", "cell from")
call write_dataset(filter_group, "n_bins", this % n_bins)
allocate(cell_ids(size(this % cells)))
@ -822,40 +820,40 @@ contains
!===============================================================================
! CellToFilter methods
!===============================================================================
subroutine get_next_bin_cell_to(this, p, estimator, current_bin, &
next_bin, weight)
class(CellToFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
integer, value, intent(in) :: current_bin
integer, intent(out) :: next_bin
real(8), intent(out) :: weight
integer :: i, start
! Find the coordinate level of the last bin we found.
if (current_bin == NO_BIN_FOUND) then
start = 1
else
do i = 1, p % n_coord
if (p % coord(i) % cell == this % cells(current_bin)) then
start = i + 1
exit
end if
end do
end if
! Starting one coordinate level deeper, find the next bin.
next_bin = NO_BIN_FOUND
weight = ERROR_REAL
do i = start, p % n_coord
if (this % map % has_key(p % coord(i) % cell)) then
next_bin = this % map % get_key(p % coord(i) % cell)
weight = ONE
exit
end if
end do
end subroutine get_next_bin_cell_to
! subroutine get_next_bin_cell_to(this, p, estimator, current_bin, &
! next_bin, weight)
! class(CellToFilter), intent(in) :: this
! type(Particle), intent(in) :: p
! integer, intent(in) :: estimator
! integer, value, intent(in) :: current_bin
! integer, intent(out) :: next_bin
! real(8), intent(out) :: weight
!
! integer :: i, start
!
! ! Find the coordinate level of the last bin we found.
! if (current_bin == NO_BIN_FOUND) then
! start = 1
! else
! do i = 1, p % n_coord
! if (p % coord(i) % cell == this % cells(current_bin)) then
! start = i + 1
! exit
! end if
! end do
! end if
!
! ! Starting one coordinate level deeper, find the next bin.
! next_bin = NO_BIN_FOUND
! weight = ERROR_REAL
! do i = start, p % n_coord
! if (this % map % has_key(p % coord(i) % cell)) then
! next_bin = this % map % get_key(p % coord(i) % cell)
! weight = ONE
! exit
! end if
! end do
! end subroutine get_next_bin_cell_to
subroutine to_statepoint_cell_to(this, filter_group)
class(CellToFilter), intent(in) :: this
@ -864,7 +862,7 @@ contains
integer :: i
integer, allocatable :: cell_ids(:)
call write_dataset(filter_group, "type", "cell")
call write_dataset(filter_group, "type", "cell to")
call write_dataset(filter_group, "n_bins", this % n_bins)
allocate(cell_ids(size(this % cells)))
@ -874,28 +872,28 @@ contains
call write_dataset(filter_group, "bins", cell_ids)
end subroutine to_statepoint_cell_to
subroutine initialize_cell_to(this)
class(CellToFilter), intent(inout) :: this
integer :: i, id
! Convert ids to indices.
do i = 1, this % n_bins
id = this % cells(i)
if (cell_dict % has_key(id)) then
this % cells(i) = cell_dict % get_key(id)
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
&// " specified on tally filter.")
end if
end do
! Generate mapping from cell indices to filter bins.
do i = 1, this % n_bins
call this % map % add_key(this % cells(i), i)
end do
end subroutine initialize_cell_to
! subroutine initialize_cell_to(this)
! class(CellToFilter), intent(inout) :: this
!
! integer :: i, id
!
! ! Convert ids to indices.
! do i = 1, this % n_bins
! id = this % cells(i)
!
! if (cell_dict % has_key(id)) then
! this % cells(i) = cell_dict % get_key(id)
! else
! call fatal_error("Could not find cell " // trim(to_str(id)) &
! &// " specified on tally filter.")
! end if
! end do
!
! ! Generate mapping from cell indices to filter bins.
! do i = 1, this % n_bins
! call this % map % add_key(this % cells(i), i)
! end do
! end subroutine initialize_cell_to
function text_label_cell_to(this, bin) result(label)
class(CellToFilter), intent(in) :: this

10
src/tracking.F90
View file

@ -15,7 +15,7 @@ module tracking
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current, &
score_track_derivative, &
score_track_derivative, score_partial_current, &
score_collision_derivative, zero_flux_derivs
use track_output, only: initialize_particle_track, write_particle_track, &
add_particle_track, finalize_particle_track
@ -166,8 +166,16 @@ contains
! Saving last cell for tallying purposes
p % last_cell = last_cell
! Update last_ information. This is needed to use the same filters as
! the ones implemented for regular tallies
p % last_uvw = p % coord(p % n_coord) % uvw
p % last_E = p % E
call cross_surface(p, last_cell)
p % event = EVENT_SURFACE
! Score cell to cell partial currents
call score_partial_current(p)
end if
else
! ====================================================================