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

Add Madland fission-Q support to F90

Browse source
This commit is contained in:
Sterling Harper 2016-07-29 13:54:30 -05:00
parent a38f53d26f
commit 9acac0725a
6 changed files with 209 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

27
openmc/data/fission_energy.py
View file

@ -62,6 +62,19 @@ class FissionEnergyRelease(object):
return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons,
self.prompt_photons, self.delayed_photons, self.betas,
self.neutrinos])
@property
def prompt_q(self):
return Sum([self.fragments, self.prompt_neutrons, self.prompt_photons,
lambda E: -E])
@property
def recoverable_q(self):
return Sum([self.recoverable, lambda E: -E])
@property
def total_q(self):
return Sum([self.total, lambda E: -E])
@property
def form(self):
@ -272,6 +285,20 @@ class FissionEnergyRelease(object):
data=self.delayed_photons.coef)
group.create_dataset('betas', data=self.betas.coef)
group.create_dataset('neutrinos', data=self.neutrinos.coef)
q_prompt = (self.fragments + self.prompt_neutrons +
self.prompt_photons + Polynomial((-1.0, 0.0)))
group.create_dataset('q_prompt', data=q_prompt.coef)
q_recoverable = (self.fragments + self.prompt_neutrons +
self.delayed_neutrons + self.prompt_photons +
self.delayed_photons + self.betas +
Polynomial((-1.0, 0.0)))
group.create_dataset('q_recoverable', data=q_recoverable.coef)
q_total = (self.fragments + self.prompt_neutrons +
self.delayed_neutrons + self.prompt_photons +
self.delayed_photons + self.betas + self.neutrinos +
Polynomial((-1.0, 0.0)))
group.create_dataset('q_total', data=q_total.coef)
elif self.form == 'Sher-Beck':
group.attrs['format'] = np.string_('Sher-Beck')
self.fragments.to_hdf5(group, 'fragments')

6
src/constants.F90
View file

@ -289,7 +289,7 @@ module constants
EVENT_ABSORB = 2
! Tally score type
integer, parameter :: N_SCORE_TYPES = 20
integer, parameter :: N_SCORE_TYPES = 22
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
@ -310,7 +310,9 @@ module constants
SCORE_NU_SCATTER_YN = -17, & ! angular flux-weighted nu-scattering moment (0:N)
SCORE_EVENTS = -18, & ! number of events
SCORE_DELAYED_NU_FISSION = -19, & ! delayed neutron production rate
SCORE_INVERSE_VELOCITY = -20 ! flux-weighted inverse velocity
SCORE_INVERSE_VELOCITY = -20, & ! flux-weighted inverse velocity
SCORE_FISS_Q_PROMPT = -21, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -22 ! recoverable fission Q-value
! Maximum scattering order supported
integer, parameter :: MAX_ANG_ORDER = 10

4
src/input_xml.F90
View file

@ -3629,6 +3629,10 @@ contains
t % score_bins(j) = SCORE_KAPPA_FISSION
case ('inverse-velocity')
t % score_bins(j) = SCORE_INVERSE_VELOCITY
case ('fission-q-prompt')
t % score_bins(j) = SCORE_FISS_Q_PROMPT
case ('fission-q-recoverable')
t % score_bins(j) = SCORE_FISS_Q_RECOV
case ('current')
t % score_bins(j) = SCORE_CURRENT
t % type = TALLY_SURFACE_CURRENT

34
src/nuclide_header.F90
View file

@ -88,6 +88,10 @@ module nuclide_header
type(DictIntInt) :: reaction_index ! map MT values to index in reactions
! array; used at tally-time
! Fission energy release
class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas
class(Function1D), allocatable :: fission_q_recov ! neutrons, gammas, betas
contains
procedure :: clear => nuclide_clear
procedure :: print => nuclide_print
@ -192,6 +196,8 @@ module nuclide_header
integer(HID_T) :: rxs_group
integer(HID_T) :: rx_group
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
integer(SIZE_T) :: name_len, name_file_len
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
@ -251,8 +257,8 @@ module nuclide_header
call this % urr_data % from_hdf5(urr_group)
! if the inelastic competition flag indicates that the inelastic cross
! section should be determined from a normal reaction cross section, we need
! to get the index of the reaction
! section should be determined from a normal reaction cross section, we
! need to get the index of the reaction
if (this % urr_data % inelastic_flag > 0) then
do i = 1, size(this % reactions)
if (this % reactions(i) % MT == this % urr_data % inelastic_flag) then
@ -296,6 +302,30 @@ module nuclide_header
call close_group(nu_group)
end if
! Read fission energy release data if present
call h5ltpath_valid_f(group_id, 'fission_energy_release', .true., exists, &
hdf5_err)
if (exists) then
fer_group = open_group(group_id, 'fission_energy_release')
call read_attribute(temp, fer_group, 'format')
if (temp == 'Madland') then
! The data uses the Madland format, i.e. polynomials
! Read the prompt Q-value
allocate(Polynomial :: this % fission_q_prompt)
fer_dset = open_dataset(fer_group, 'q_prompt')
call this % fission_q_prompt % from_hdf5(fer_dset)
call close_dataset(fer_dset)
! Read the recoverable energy Q-value
allocate(Polynomial :: this % fission_q_recov)
fer_dset = open_dataset(fer_group, 'q_recoverable')
call this % fission_q_recov % from_hdf5(fer_dset)
call close_dataset(fer_dset)
end if
call close_group(fer_group)
end if
! Create derived cross section data
call this % create_derived()

2
src/output.F90
View file

@ -791,6 +791,8 @@ contains
score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment"
score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity"
score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power"
score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power"
! Create filename for tally output
filename = trim(path_output) // "tallies.out"

162
src/tally.F90
View file

@ -625,14 +625,14 @@ contains
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scale by kappa-fission
associate (nuc => nuclides(p%event_nuclide))
if (micro_xs(p%event_nuclide)%absorption > ZERO .and. &
nuc%fissionable) then
score = p%absorb_wgt * &
nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(p%event_nuclide)%fission / &
micro_xs(p%event_nuclide)%absorption
! fission scaled by kappa-fission
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
nuc % fissionable) then
score = p % absorb_wgt * &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(p % event_nuclide) % fission / &
micro_xs(p % event_nuclide) % absorption
end if
end associate
else
@ -641,12 +641,12 @@ contains
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
associate (nuc => nuclides(p%event_nuclide))
if (nuc%fissionable) then
score = p%last_wgt * &
nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(p%event_nuclide)%fission / &
micro_xs(p%event_nuclide)%absorption
associate (nuc => nuclides(p % event_nuclide))
if (nuc % fissionable) then
score = p % last_wgt * &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(p % event_nuclide) % fission / &
micro_xs(p % event_nuclide) % absorption
end if
end associate
end if
@ -654,22 +654,23 @@ contains
else
if (i_nuclide > 0) then
associate (nuc => nuclides(i_nuclide))
if (nuc%fissionable) then
score = nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(i_nuclide)%fission * atom_density * flux
if (nuc % fissionable) then
score = nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(i_nuclide) % fission * atom_density * flux
end if
end associate
else
do l = 1, materials(p%material)%n_nuclides
do l = 1, materials(p % material) % n_nuclides
! Determine atom density and index of nuclide
atom_density_ = materials(p%material)%atom_density(l)
i_nuc = materials(p%material)%nuclide(l)
atom_density_ = materials(p % material) % atom_density(l)
i_nuc = materials(p % material) % nuclide(l)
! If nuclide is fissionable, accumulate kappa fission
associate(nuc => nuclides(i_nuc))
if (nuc % fissionable) then
score = score + nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(i_nuc)%fission * atom_density_ * flux
score = score + &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(i_nuc) % fission * atom_density_ * flux
end if
end associate
end do
@ -694,6 +695,123 @@ contains
end if
end if
case (SCORE_FISS_Q_PROMPT)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scaled by Q-value
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
allocated(nuc % fission_q_prompt)) then
score = p % absorb_wgt &
* nuc % fission_q_prompt % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
associate (nuc => nuclides(p % event_nuclide))
if (allocated(nuc % fission_q_prompt)) then
score = p % last_wgt &
* nuc % fission_q_prompt % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption
end if
end associate
end if
else
if (t % estimator == ESTIMATOR_COLLISION) then
E = p % last_E
else
E = p % E
end if
if (i_nuclide > 0) then
if (allocated(nuclides(i_nuclide) % fission_q_prompt)) then
score = micro_xs(i_nuclide) % fission * atom_density * flux &
* nuclides(i_nuclide) % fission_q_prompt % evaluate(E)
else
score = ZERO
end if
else
score = ZERO
do l = 1, materials(p % material) % n_nuclides
atom_density_ = materials(p % material) % atom_density(l)
i_nuc = materials(p % material) % nuclide(l)
if (allocated(nuclides(i_nuc) % fission_q_prompt)) then
score = score + micro_xs(i_nuc) % fission * atom_density_ &
* flux &
* nuclides(i_nuc) % fission_q_prompt % evaluate(E)
end if
end do
end if
end if
case (SCORE_FISS_Q_RECOV)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scaled by Q-value
associate (nuc => nuclides(p % event_nuclide))
if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
allocated(nuc % fission_q_recov)) then
score = p % absorb_wgt &
* nuc % fission_q_recov % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
associate (nuc => nuclides(p % event_nuclide))
if (allocated(nuc % fission_q_recov)) then
score = p % last_wgt &
* nuc % fission_q_recov % evaluate(p % last_E) &
* micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption
end if
end associate
end if
else
if (t % estimator == ESTIMATOR_COLLISION) then
E = p % last_E
else
E = p % E
end if
if (i_nuclide > 0) then
if (allocated(nuclides(i_nuclide) % fission_q_recov)) then
score = micro_xs(i_nuclide) % fission * atom_density * flux &
* nuclides(i_nuclide) % fission_q_recov % evaluate(E)
else
score = ZERO
end if
else
score = ZERO
do l = 1, materials(p % material) % n_nuclides
atom_density_ = materials(p % material) % atom_density(l)
i_nuc = materials(p % material) % nuclide(l)
if (allocated(nuclides(i_nuc) % fission_q_recov)) then
score = score + micro_xs(i_nuc) % fission * atom_density_ &
* flux * nuclides(i_nuc) % fission_q_recov % evaluate(E)
end if
end do
end if
end if
case default
if (t % estimator == ESTIMATOR_ANALOG) then
! Any other score is assumed to be a MT number. Thus, we just need