From 9acac0725a26efc14906d692bc99f0dacabeccea Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 29 Jul 2016 13:54:30 -0500 Subject: [PATCH] Add Madland fission-Q support to F90 --- openmc/data/fission_energy.py | 27 ++++++ src/constants.F90 | 6 +- src/input_xml.F90 | 4 + src/nuclide_header.F90 | 34 ++++++- src/output.F90 | 2 + src/tally.F90 | 162 +++++++++++++++++++++++++++++----- 6 files changed, 209 insertions(+), 26 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 5716a3e85..7415a9b51 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -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') diff --git a/src/constants.F90 b/src/constants.F90 index d2bb6e1cf..5447b5cc2 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3c8291e9c..f661927d7 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8ff83482e..31e069b68 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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() diff --git a/src/output.F90 b/src/output.F90 index 09df23f0d..35524a9a3 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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" diff --git a/src/tally.F90 b/src/tally.F90 index b40a39ce0..1271da7da 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -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