Move score_fission_eout to C++

This commit is contained in:
Sterling Harper 2019-02-05 23:16:41 -05:00
parent 64121662af
commit ee4f4e84aa
2 changed files with 165 additions and 197 deletions

View file

@ -55,6 +55,14 @@ module tally
integer(C_INT), value :: score_index
end subroutine
subroutine score_fission_eout(p, i_tally, i_score, score_bin) bind(C)
import Particle, C_INT
type(Particle) :: p
integer(C_INT), value :: i_tally
integer(C_INT), value :: i_score
integer(C_INT), value :: score_bin
end subroutine
subroutine score_analog_tally_ce(p) bind(C)
import Particle
type(Particle), intent(in) :: p
@ -1229,7 +1237,7 @@ contains
! Add derivative information on score for differential tallies.
if (t % deriv() /= C_NONE) then
call apply_derivative_to_score(p, t, i_nuclide, atom_density, &
call apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, &
score_bin, score)
end if
@ -2037,205 +2045,18 @@ contains
end associate
end subroutine score_general_mg
!===============================================================================
! SCORE_FISSION_EOUT handles a special case where we need to store neutron
! production rate with an outgoing energy filter (think of a fission matrix). In
! this case, we may need to score to multiple bins if there were multiple
! neutrons produced with different energies.
!===============================================================================
subroutine score_fission_eout(p, i_tally, i_score, score_bin)
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
integer, intent(in) :: i_score ! index for score
integer, intent(in) :: score_bin
integer :: i ! index of outgoing energy filter
integer :: j ! index of delayedgroup filter
integer :: d ! delayed group
integer :: g ! another delayed group
integer :: d_bin ! delayed group bin index
integer :: k ! loop index for bank sites
integer :: l ! loop index for tally filters
integer :: f ! index in filters array
integer :: b ! index of filter bin
integer :: i_match ! matching bin index on energyout filter
integer :: i_bin ! index of matching filter bin
integer :: bin_energyout ! original outgoing energy bin
integer :: i_filter ! index for matching filter bin combination
real(8) :: filter_weight ! combined weight of all filters
real(8) :: score ! actual score
real(8) :: E_out ! energy of fission bank site
integer :: g_out ! energy group of fission bank site
associate (t => tallies(i_tally) % obj)
! save original outgoing energy bin and score index
i = t % filter(t % energyout_filter()) + 1
i_bin = filter_matches(i) % i_bin()
bin_energyout = filter_matches(i) % bins_data(i_bin)
! declare the energyout filter type
select type(eo_filt => filters(i) % obj)
type is (EnergyoutFilter)
! Since the creation of fission sites is weighted such that it is
! expected to create n_particles sites, we need to multiply the
! score by keff to get the true nu-fission rate. Otherwise, the sum
! of all nu-fission rates would be ~1.0.
! loop over number of particles banked
do k = 1, p % n_bank
! get the delayed group
g = fission_bank_delayed_group(n_bank - p % n_bank + k)
! determine score based on bank site weight and keff
score = keff * fission_bank_wgt(n_bank - p % n_bank + k)
! Add derivative information for differential tallies. Note that the
! i_nuclide and atom_density arguments do not matter since this is an
! analog estimator.
if (t % deriv() /= C_NONE) then
call apply_derivative_to_score(p, t, 0, ZERO, SCORE_NU_FISSION, score)
end if
if (.not. run_CE .and. eo_filt % matches_transport_groups) then
! determine outgoing energy group from fission bank
g_out = int(fission_bank_E(n_bank - p % n_bank + k))
! modify the value so that g_out = 1 corresponds to the highest
! energy bin
g_out = eo_filt % n_bins - g_out + 1
! change outgoing energy bin
call filter_matches(i) % bins_set_data(i_bin, g_out)
else
! determine outgoing energy from fission bank
if (run_CE) then
E_out = fission_bank_E(n_bank - p % n_bank + k)
else
E_out = energy_bin_avg(int(fission_bank_E(n_bank - p % n_bank + k)))
end if
! If this outgoing energy falls within the energyout filter's range,
! set the appropriate filter_matches bin.
i_match = eo_filt % search(E_out)
if (i_match == -1) cycle
call filter_matches(i) % bins_set_data(i_bin, i_match)
end if
! Case for tallying prompt neutrons
if (score_bin == SCORE_NU_FISSION .or. &
(score_bin == SCORE_PROMPT_NU_FISSION .and. g == 0)) then
! determine scoring index and weight for this filter combination
i_filter = 1
do l = 1, t % n_filters()
i_filter = i_filter + (filter_matches(t % filter(l) + 1) &
% bins_data(filter_matches(t % filter(l) + 1) % i_bin()) - 1) * &
t % stride(l)
end do
! Add score to tally
!$omp atomic
t % results(RESULT_VALUE, i_score, i_filter) = &
t % results(RESULT_VALUE, i_score, i_filter) + score
! Case for tallying delayed emissions
else if (score_bin == SCORE_DELAYED_NU_FISSION .and. g /= 0) then
! Get the index of delayed group filter
j = t % delayedgroup_filter()
! if the delayed group filter is present, tally to corresponding
! delayed group bin if it exists
if (j > 0) then
! declare the delayed group filter type
select type(dg_filt => filters(t % filter(j) + 1) % obj)
type is (DelayedGroupFilter)
! loop over delayed group bins until the corresponding bin is
! found
do d_bin = 1, dg_filt % n_bins
d = dg_filt % groups(d_bin)
! check whether the delayed group of the particle is equal to
! the delayed group of this bin
if (d == g) then
! Reset scoring index and filter weight
i_filter = 1
filter_weight = ONE
! determine scoring index and weight for this filter
! combination
do l = 1, t % n_filters()
f = t % filter(l) + 1
b = filter_matches(f) % i_bin()
i_filter = i_filter + (filter_matches(f) &
% bins_data(b) - 1) * t % stride(l)
filter_weight = filter_weight * filter_matches(f) &
% weights_data(b)
end do
call score_fission_delayed_dg(i_tally, d_bin, &
score * filter_weight, i_score)
end if
end do
end select
! if the delayed group filter is not present, add score to tally
else
! Reset scoring index and filter weight
i_filter = 1
filter_weight = ONE
! determine scoring index and weight for this filter combination
do l = 1, t % n_filters()
f = t % filter(l) + 1
b = filter_matches(f) % i_bin()
i_filter = i_filter + (filter_matches(f) % bins_data(b) - 1) &
* t % stride(l)
filter_weight = filter_weight * filter_matches(f) &
% weights_data(b)
end do
! Add score to tally
!$omp atomic
t % results(RESULT_VALUE, i_score, i_filter) = &
t % results(RESULT_VALUE, i_score, i_filter) + score * filter_weight
end if
end if
end do
end select
! reset outgoing energy bin and score index
call filter_matches(i) % bins_set_data(i_bin, bin_energyout)
end associate
end subroutine score_fission_eout
!===============================================================================
! APPLY_DERIVATIVE_TO_SCORE multiply the given score by its relative derivative
!===============================================================================
subroutine apply_derivative_to_score(p, t, i_nuclide, atom_density, &
score_bin, score)
type(Particle), intent(in) :: p
type(TallyObject), intent(in) :: t
integer, intent(in) :: i_nuclide
real(8), intent(in) :: atom_density ! atom/b-cm
integer, intent(in) :: score_bin
real(8), intent(inout) :: score
subroutine apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, &
score_bin, score) bind(C)
type(Particle), intent(in) :: p
integer(C_INT), value, intent(in) :: i_tally
integer(C_INT), value, intent(in) :: i_nuclide
real(C_DOUBLE), value, intent(in) :: atom_density ! atom/b-cm
integer(C_INT), value, intent(in) :: score_bin
real(C_DOUBLE), intent(inout) :: score
type(TallyDerivative), pointer :: deriv
integer :: l
@ -2244,6 +2065,8 @@ contains
real(8) :: flux_deriv
real(8) :: dsig_s, dsig_a, dsig_f, cum_dsig
associate (t => tallies(i_tally) % obj)
if (score == ZERO) return
! If our score was previously c then the new score is
@ -2778,7 +2601,7 @@ contains
&analog and collision estimators.")
end select
end select
!end associate
end associate
end subroutine apply_derivative_to_score
!===============================================================================

View file

@ -44,6 +44,10 @@ extern "C" void
score_general_mg(Particle* p, int i_tally, int start_index, int filter_index,
int i_nuclide, double atom_density, double flux);
extern "C" void
apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
double atom_density, int score_bin, double* score);
extern "C" int
energy_filter_search(const EnergyFilter* filt, double val);
@ -699,6 +703,147 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
dg_match.bins_[i_bin-1] = original_bin;
}
//! Helper function for nu-fission tallies with energyout filters.
//
//! In this case, we may need to score to multiple bins if there were multiple
//! neutrons produced with different energies.
extern "C" void
score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin)
{
//TODO: off-by-one
const Tally& tally {*model::tallies[i_tally-1]};
auto results = tally_results(i_tally);
auto i_eout_filt = tally.filters()[tally.energyout_filter_-1];
auto i_bin = simulation::filter_matches[i_eout_filt].i_bin_;
auto bin_energyout = simulation::filter_matches[i_eout_filt].bins_[i_bin-1];
const EnergyoutFilter& eo_filt
{*dynamic_cast<EnergyoutFilter*>(model::tally_filters[i_eout_filt].get())};
// Note that the score below is weighted by keff. Since the creation of
// fission sites is weighted such that it is expected to create n_particles
// sites, we need to multiply the score by keff to get the true nu-fission
// rate. Otherwise, the sum of all nu-fission rates would be ~1.0.
// loop over number of particles banked
for (auto i = 0; i < p->n_bank; ++i) {
auto i_bank = simulation::n_bank - p->n_bank + i;
const auto& bank = simulation::fission_bank[i_bank];
// get the delayed group
auto g = bank.delayed_group;
// determine score based on bank site weight and keff
double score = simulation::keff * bank.wgt;
// Add derivative information for differential tallies. Note that the
// i_nuclide and atom_density arguments do not matter since this is an
// analog estimator.
if (tally.deriv_ != C_NONE)
apply_derivative_to_score(p, i_tally, 0, 0., SCORE_NU_FISSION, &score);
if (!settings::run_CE && eo_filt.matches_transport_groups_) {
// determine outgoing energy group from fission bank
auto g_out = static_cast<int>(bank.E);
// modify the value so that g_out = 1 corresponds to the highest energy
// bin
g_out = eo_filt.n_bins_ - g_out + 1;
// change outgoing energy bin
simulation::filter_matches[i_eout_filt].bins_[i_bin-1] = g_out;
} else {
double E_out;
if (settings::run_CE) {
E_out = bank.E;
} else {
E_out = data::energy_bin_avg[static_cast<int>(bank.E)];
}
//TODO: do this without the extern "C" function
auto i_match = energy_filter_search(&eo_filt, E_out);
if (i_match == -1) continue;
simulation::filter_matches[i_eout_filt].bins_[i_bin-1] = i_match;
}
// Case for tallying prompt neutrons
if (score_bin == SCORE_NU_FISSION
|| (score_bin == SCORE_PROMPT_NU_FISSION && g == 0)) {
// Find the filter scoring index for this filter combination
//TODO: should this include a weight?
int filter_index = 1;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
auto& match {simulation::filter_matches[i_filt]};
auto i_bin = match.i_bin_;
filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(j);
}
// Update tally results
#pragma omp atomic
results(filter_index-1, i_score-1, RESULT_VALUE) += score;
} else if (score_bin == SCORE_DELAYED_NU_FISSION && g != 0) {
// Get the index of the delayed group filter
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_-1];
// If the delayed group filter is present, tally to corresponding delayed
// group bin if it exists
if (i_dg_filt >= 0) {
const DelayedGroupFilter& dg_filt {*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
// Loop over delayed group bins until the corresponding bin is found
for (auto d_bin = 0; d_bin < dg_filt.n_bins_; ++d_bin) {
if (dg_filt.groups_[d_bin] == g) {
// Find the filter index and weight for this filter combination
int filter_index = 1;
double filter_weight = 1.;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
auto& match {simulation::filter_matches[i_filt]};
auto i_bin = match.i_bin_;
filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(j);
filter_weight *= match.weights_[i_bin-1];
}
score_fission_delayed_dg(i_tally, d_bin+1, score*filter_weight,
i_score);
}
}
// If the delayed group filter is not present, add score to tally
} else {
// Find the filter index and weight for this filter combination
int filter_index = 1;
double filter_weight = 1.;
for (auto j = 0; j < tally.filters().size(); ++j) {
auto i_filt = tally.filters(j);
auto& match {simulation::filter_matches[i_filt]};
auto i_bin = match.i_bin_;
filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(j);
filter_weight *= match.weights_[i_bin-1];
}
// Update tally results
#pragma omp atomic
results(filter_index-1, i_score-1, RESULT_VALUE) += score*filter_weight;
}
}
}
// Reset outgoing energy bin and score index
simulation::filter_matches[i_eout_filt].bins_[i_bin-1] = bin_energyout;
}
//! Tally rates for when the user requests a tally on all nuclides.
void