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

Set event_MT during photon transport, allow for photon energy cutoff

Browse source
This commit is contained in:
Paul Romano 2017-06-29 11:44:03 -05:00
parent 974b714565
commit 686362248b
7 changed files with 79 additions and 57 deletions
Download patch file Download diff file Expand all files Collapse all files

29
openmc/settings.py
View file

@ -577,16 +577,15 @@ class Settings(object):
raise ValueError(msg)
for key in cutoff:
if key == 'weight':
cv.check_type('weight cutoff', cutoff['weight'], Real)
cv.check_greater_than('weight cutoff', cutoff['weight'], 0.0)
cv.check_type('weight cutoff', cutoff[key], Real)
cv.check_greater_than('weight cutoff', cutoff[key], 0.0)
elif key == 'weight_avg':
cv.check_type('average survival weight', cutoff['weight_avg'],
Real)
cv.check_type('average survival weight', cutoff[key], Real)
cv.check_greater_than('average survival weight',
cutoff['weight_avg'], 0.0)
elif key == 'energy':
cv.check_type('energy cutoff', cutoff['energy'], Real)
cv.check_greater_than('energy cutoff', cutoff['energy'], 0.0)
cutoff[key], 0.0)
elif key in ['energy', 'energy_photon']:
cv.check_type('energy cutoff', cutoff[key], Real)
cv.check_greater_than('energy cutoff', cutoff[key], 0.0)
else:
msg = 'Unable to set cutoff to "{0}" which is unsupported by '\
'OpenMC'.format(key)
@ -946,17 +945,9 @@ class Settings(object):
def _create_cutoff_subelement(self, root):
if self._cutoff is not None:
element = ET.SubElement(root, "cutoff")
if 'weight' in self._cutoff:
subelement = ET.SubElement(element, "weight")
subelement.text = str(self._cutoff['weight'])
if 'weight_avg' in self._cutoff:
subelement = ET.SubElement(element, "weight_avg")
subelement.text = str(self._cutoff['weight_avg'])
if 'energy' in self._cutoff:
subelement = ET.SubElement(element, "energy")
subelement.text = str(self._cutoff['energy'])
for key, value in self.items():
subelement = ET.SubElement(element, key)
subelement.text = str(value)
def _create_entropy_subelement(self, root):
if self._entropy_mesh is not None:

4
src/constants.F90
View file

@ -226,7 +226,9 @@ module constants
N_3HEA = 193, N_4N2P = 194, N_4N2A = 195, N_4NPA = 196, N_3P = 197, &
N_N3P = 198, N_3N2PA = 199, N_5N2P = 200, N_P0 = 600, N_PC = 649, &
N_D0 = 650, N_DC = 699, N_T0 = 700, N_TC = 749, N_3HE0 = 750, &
N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891
N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891, &
COHERENT = 502, INCOHERENT = 504, PHOTOELECTRIC = 522, &
PAIR_PROD_ELEC = 515, PAIR_PROD = 516, PAIR_PROD_NUC = 517
! ACE table types
integer, parameter :: &

14
src/endf.F90
View file

@ -176,6 +176,20 @@ contains
string = '(n,Xa)'
case (444)
string = '(damage)'
case (COHERENT)
string = 'coherent scatter'
case (INCOHERENT)
string = 'incoherent scatter'
case (PAIR_PROD_ELEC)
string = 'pair production, electron'
case (PAIR_PROD)
string = 'pair production'
case (PAIR_PROD_NUC)
string = 'pair production, nuclear'
case (PHOTOELECTRIC)
string = 'photoelectric'
case (534 : 572)
string = 'photoelectric, ' // trim(SUBSHELLS(MT - 533)) // ' subshell'
case (600 : 648)
string = '(n,p' // trim(to_str(MT-600)) // ')'
case (649)

2
src/global.F90
View file

@ -313,7 +313,7 @@ module global
logical :: survival_biasing = .false.
real(8) :: weight_cutoff = 0.25_8
real(8) :: energy_cutoff = ZERO
real(8) :: energy_cutoff(3) = [ZERO, 1000.0_8, ZERO]
real(8) :: weight_survive = ONE
! ============================================================================

5
src/input_xml.F90
View file

@ -599,7 +599,10 @@ contains
call get_node_value(node_cutoff, "weight_avg", weight_survive)
end if
if (check_for_node(node_cutoff, "energy")) then
call get_node_value(node_cutoff, "energy", energy_cutoff)
call get_node_value(node_cutoff, "energy", energy_cutoff(1))
end if
if (check_for_node(node_cutoff, "energy_photon")) then
call get_node_value(node_cutoff, "energy_photon", energy_cutoff(2))
end if
end if

17
src/photon_header.F90
View file

@ -1,6 +1,6 @@
module photon_header
use hdf5, only: HID_T, HSIZE_T
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use constants, only: ZERO, HALF, SUBSHELLS
use dict_header, only: DictIntInt
@ -10,7 +10,7 @@ module photon_header
real(8), allocatable :: compton_profile_pz(:)
type ElectronSubshell
character(3) :: label
integer :: index_subshell ! index in SUBSHELLS
integer :: threshold
real(8) :: n_electrons
real(8) :: binding_energy
@ -24,7 +24,8 @@ module photon_header
end type ElectronSubshell
type PhotonInteraction
integer :: Z ! atomic number
character(3) :: name ! atomic symbol, e.g. 'Zr'
integer :: Z ! atomic number
! Microscopic cross sections
real(8), allocatable :: energy(:)
@ -82,15 +83,22 @@ contains
integer(HID_T) :: rgroup, tgroup
integer(HID_T) :: dset_id
integer(HSIZE_T) :: dims(1), dims2(2)
integer(SIZE_T) :: name_len
integer :: n_energy
integer :: n_shell
integer :: n_profile
integer :: n_transition
! integer, allocatable :: designators(:)
character(3), allocatable :: designators(:)
real(8) :: c
real(8), allocatable :: matrix(:,:)
! Get name of nuclide from group
name_len = len(this % name)
this % name = get_name(group_id, name_len)
! Get rid of leading '/'
this % name = trim(this % name(2:))
! Get atomic number
call read_attribute(this % Z, group_id, 'Z')
@ -160,6 +168,7 @@ contains
do j = 1, size(SUBSHELLS)
if (designators(i) == SUBSHELLS(j)) then
call this % shell_dict % add_key(j, i)
this % shells(i) % index_subshell = j
exit
end if
end do

65
src/physics.F90
View file

@ -45,42 +45,42 @@ contains
! Sample reaction for the material the particle is in
if (p % type == NEUTRON) then
call sample_reaction(p)
elseif (p % type == PHOTON) then
call sample_neutron_reaction(p)
else
call sample_photon_reaction(p)
end if
! Kill particle if energy falls below cutoff
if (p % E < energy_cutoff(p % type)) then
p % alive = .false.
p % wgt = ZERO
p % last_wgt = ZERO
end if
! Display information about collision
if (verbosity >= 10 .or. trace) then
call write_message(" " // trim(reaction_name(p % event_MT)) &
&// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) &
&// ". Energy = " // trim(to_str(p % E)) // " eV.")
if (p % type == NEUTRON) then
call write_message(" " // trim(reaction_name(p % event_MT)) &
&// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) &
&// ". Energy = " // trim(to_str(p % E)) // " eV.")
else
call write_message(" " // trim(reaction_name(p % event_MT)) &
&// " with " // trim(adjustl(elements(p % event_nuclide) % name)) &
&// ". Energy = " // trim(to_str(p % E)) // " eV.")
end if
end if
! check for very low energy
if (p % E < 1.0e-100_8) then
p % alive = .false.
if (master) call warning("Killing neutron with extremely low energy")
end if
! Advance URR seed stream 'N' times after energy changes
if (p % E /= p % last_E) then
call prn_set_stream(STREAM_URR_PTABLE)
call advance_prn_seed(size(nuclides, kind=8))
call prn_set_stream(STREAM_TRACKING)
endif
end subroutine collision
!===============================================================================
! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for
! each nuclide within a material and then samples a reaction for that nuclide
! and calls the appropriate routine to process the physics. Note that there is
! special logic when suvival biasing is turned on since fission and
! disappearance are treated implicitly.
! SAMPLE_NEUTRON_REACTION samples a nuclide based on the macroscopic cross
! sections for each nuclide within a material and then samples a reaction for
! that nuclide and calls the appropriate routine to process the physics. Note
! that there is special logic when suvival biasing is turned on since fission
! and disappearance are treated implicitly.
!===============================================================================
subroutine sample_reaction(p)
subroutine sample_neutron_reaction(p)
type(Particle), intent(inout) :: p
@ -128,20 +128,19 @@ contains
call scatter(p, i_nuclide, i_nuc_mat)
! Play russian roulette if survival biasing is turned on
if (survival_biasing) then
call russian_roulette(p)
if (.not. p % alive) return
end if
! Kill neutron under certain energy
if (p % E < energy_cutoff) then
p % alive = .false.
p % wgt = ZERO
p % last_wgt = ZERO
! Advance URR seed stream 'N' times after energy changes
if (p % E /= p % last_E) then
call prn_set_stream(STREAM_URR_PTABLE)
call advance_prn_seed(size(nuclides, kind=8))
call prn_set_stream(STREAM_TRACKING)
end if
end subroutine sample_reaction
end subroutine sample_neutron_reaction
!===============================================================================
! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross
@ -184,6 +183,7 @@ contains
if (prob > cutoff) then
call rayleigh_scatter(elm, alpha, mu)
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
p % event_MT = COHERENT
return
end if
@ -193,6 +193,7 @@ contains
call compton_scatter(elm, alpha, alpha_out, mu, .true.)
p % E = alpha_out*MASS_ELECTRON
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
p % event_MT = INCOHERENT
return
end if
@ -219,6 +220,7 @@ contains
! E_electron = p % E - elm % shells(i_shell) % binding_energy
call atomic_relaxation(p, elm, i_shell)
p % event_MT = 533 + elm % shells(i_shell) % index_subshell
p % alive = .false.
return
end if
@ -239,6 +241,7 @@ contains
! Set energy
p % E = MASS_ELECTRON
p % event_MT = PAIR_PROD
! Create photon in opposite direction
call p % create_secondary(-p % coord(1) % uvw, MASS_ELECTRON, &