Move score_analog_tally_mg to C++

This commit is contained in:
Sterling Harper 2019-01-30 22:18:54 -05:00
parent 6c0e28f069
commit abe9bf4a6f
2 changed files with 120 additions and 139 deletions

View file

@ -51,6 +51,11 @@ module tally
type(Particle), intent(in) :: p
end subroutine
subroutine score_analog_tally_mg(p) bind(C)
import Particle
type(Particle), intent(in) :: p
end subroutine
subroutine score_tracklength_tally(p, distance) bind(C)
import Particle, C_DOUBLE
type(Particle) :: p
@ -2077,141 +2082,6 @@ contains
end associate
end subroutine score_all_nuclides
!===============================================================================
! SCORE_ANALOG_TALLY keeps track of how many events occur in a specified cell,
! energy range, etc. Note that since these are "analog" tallies, they are only
! triggered at every collision, not every event
!===============================================================================
subroutine score_analog_tally_mg(p)
type(Particle), intent(in) :: p
integer :: i, j
integer :: i_tally
integer :: i_filt
integer :: i_bin
integer :: k ! loop index for nuclide bins
integer :: filter_index ! single index for single bin
integer :: i_nuclide ! index in nuclides array
real(8) :: filter_weight ! combined weight of all filters
real(8) :: atom_density
logical :: finished ! found all valid bin combinations
type(Material), pointer :: mat
! A loop over all tallies is necessary because we need to simultaneously
! determine different filter bins for the same tally in order to score to it
TALLY_LOOP: do i = 1, active_analog_tallies_size()
! Get index of tally and pointer to tally
i_tally = active_analog_tallies_data(i)
associate (t => tallies(i_tally) % obj)
! Find all valid bins in each filter if they have not already been found
! for a previous tally.
do j = 1, t % n_filters()
i_filt = t % filter(j)
if (.not. filter_matches(i_filt) % bins_present) then
call filter_matches(i_filt) % bins_clear()
call filter_matches(i_filt) % weights_clear()
call filters(i_filt) % obj % get_all_bins(p, t % estimator(), &
filter_matches(i_filt))
filter_matches(i_filt) % bins_present = .true.
end if
! If there are no valid bins for this filter, then there is nothing to
! score and we can move on to the next tally.
if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP
! Set the index of the bin used in the first filter combination
call filter_matches(i_filt) % set_i_bin(1)
end do
! ========================================================================
! Loop until we've covered all valid bins on each of the filters.
FILTER_LOOP: do
! Reset scoring index and weight
filter_index = 1
filter_weight = ONE
! Determine scoring index and weight for this filter combination
do j = 1, t % n_filters()
i_filt = t % filter(j)
i_bin = filter_matches(i_filt) % i_bin()
filter_index = filter_index + (filter_matches(i_filt) &
% bins_data(i_bin) - 1) * t % stride(j)
filter_weight = filter_weight * filter_matches(i_filt) &
% weights_data(i_bin)
end do
! ======================================================================
! Nuclide logic
! Check for nuclide bins
NUCLIDE_LOOP: do k = 1, t % n_nuclide_bins()
! Get index of nuclide in nuclides array
i_nuclide = t % nuclide_bins(k)
if (i_nuclide > 0) then
if (p % material /= MATERIAL_VOID) then
! Get pointer to current material
mat => materials(p % material)
! Determine index of nuclide in Material % atom_density array
j = mat % mat_nuclide_index(i_nuclide)
if (j == 0) cycle NUCLIDE_LOOP
! Copy corresponding atom density
atom_density = mat % atom_density(j)
end if
else
atom_density = ZERO
end if
call score_general(p, i_tally, (k-1)*t % n_score_bins, filter_index, &
i_nuclide, atom_density, filter_weight)
end do NUCLIDE_LOOP
! ======================================================================
! Filter logic
! Increment the filter bins, starting with the last filter to find the
! next valid bin combination
finished = .true.
do j = t % n_filters(), 1, -1
i_filt = t % filter(j)
if (filter_matches(i_filt) % i_bin() < filter_matches(i_filt) % &
bins_size()) then
call filter_matches(i_filt) % set_i_bin(filter_matches(i_filt) % i_bin() + 1)
finished = .false.
exit
else
call filter_matches(i_filt) % set_i_bin(1)
end if
end do
! Once we have finished all valid bins for each of the filters, exit
! the loop.
if (finished) exit FILTER_LOOP
end do FILTER_LOOP
! If the user has specified that we can assume all tallies are spatially
! separate, this implies that once a tally has been scored to, we needn't
! check the others. This cuts down on overhead when there are many
! tallies specified
if (assume_separate) exit TALLY_LOOP
end associate
end do TALLY_LOOP
! Reset filter matches flag
filter_matches(:) % bins_present = .false.
end subroutine score_analog_tally_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

View file

@ -143,7 +143,7 @@ adaptor_type<3> tally_results(int idx)
return xt::adapt(results, size, xt::no_ownership(), shape);
}
//! Score tallies based on a simple count of events.
//! Score tallies based on a simple count of events (for continuous energy).
//
//! Analog tallies ar etriggered at every collision, not every event.
@ -210,7 +210,9 @@ score_analog_tally_ce(Particle* p)
auto i_nuclide = tally.nuclides_[i];
// Tally this event in the present nuclide bin if that bin represents
// the event nuclide or the total material.
// the event nuclide or the total material. Note that the i_nuclide
// and flux arguments for score_general are not used for analog
// tallies.
if (i_nuclide == p->event_nuclide || i_nuclide == -1)
score_general_ce(p, i_tally, i*n_score_bins, filter_index,
-1, -1., filter_weight);
@ -219,8 +221,6 @@ score_analog_tally_ce(Particle* p)
// In the case that the user has requested to tally all nuclides, we
// can take advantage of the fact that we know exactly how nuclide
// bins correspond to nuclide indices. First, tally the nuclide.
// Note that the i_nuclide and flux arguments for score_general are
// not used for analog tallies.
auto i = p->event_nuclide;
score_general_ce(p, i_tally, i*n_score_bins, filter_index,
-1, -1., filter_weight);
@ -266,6 +266,117 @@ next_tally:
match.bins_present_ = false;
}
//! Score tallies based on a simple count of events (for multigroup).
//
//! Analog tallies ar etriggered at every collision, not every event.
extern "C" void
score_analog_tally_mg(Particle* p)
{
for (auto i_tally : model::active_analog_tallies) {
//TODO: off-by-one
const Tally& tally {*model::tallies[i_tally-1]};
auto n_score_bins = tally_get_n_score_bins(i_tally);
//--------------------------------------------------------------------------
// Loop through all relevant filters and find the filter bins and weights
// for this event.
// Find all valid bins in each filter if they have not already been found
// for a previous tally.
for (auto i_filt : tally.filters()) {
//TODO: off-by-one
auto& match {simulation::filter_matches[i_filt-1]};
if (!match.bins_present_) {
match.bins_.clear();
match.weights_.clear();
//TODO: off-by-one
model::tally_filters[i_filt-1]
->get_all_bins(p, tally.estimator_, match);
match.bins_present_ = true;
}
// If there are no valid bins for this filter, then there is nothing to
// score so we can move on to the next tally.
if (match.bins_.size() == 0) goto next_tally;
// Set the index of the bin used in the first filter combination
match.i_bin_ = 1;
}
//--------------------------------------------------------------------------
// Loop over filter bins and nuclide bins
for (bool filter_loop_done = false; !filter_loop_done; ) {
//------------------------------------------------------------------------
// Filter logic
// Determine scoring index and weight for the current filter combination
int filter_index = 1;
double filter_weight = 1.;
for (auto i = 0; i < tally.filters().size(); ++i) {
auto i_filt = tally.filters(i);
//TODO: off-by-one
auto& match {simulation::filter_matches[i_filt-1]};
auto i_bin = match.i_bin_;
//TODO: off-by-one
filter_index += (match.bins_[i_bin-1] - 1) * tally.strides(i);
filter_weight *= match.weights_[i_bin-1];
}
//------------------------------------------------------------------------
// Nuclide logic
for (auto i = 0; i < tally.nuclides_.size(); ++i) {
auto i_nuclide = tally.nuclides_[i];
double atom_density = 0.;
if (i_nuclide > 0) {
auto j = material_nuclide_index(p->material, i_nuclide);
if (j == 0) continue;
atom_density = material_atom_density(p->material, j);
}
score_general_mg(p, i_tally, i*n_score_bins, filter_index,
i_nuclide, atom_density, filter_weight);
}
//------------------------------------------------------------------------
// Further filter logic
// Increment the filter bins, starting with the last filter to find the
// next valid bin combination
filter_loop_done = true;
for (int i = tally.filters().size()-1; i >= 0; --i) {
auto i_filt = tally.filters(i);
//TODO: off-by-one
auto& match {simulation::filter_matches[i_filt-1]};
if (match.i_bin_ < match.bins_.size()) {
++match.i_bin_;
filter_loop_done = false;
break;
} else {
match.i_bin_ = 1;
}
}
}
// If the user has specified that we can assume all tallies are spatially
// separate, this implies that once a tally has been scored to, we needn't
// check the others. This cuts down on overhead when there are many
// tallies specified
if (settings::assume_separate) break;
next_tally:
;
}
// Reset all the filter matches for the next tally event.
for (auto& match : simulation::filter_matches)
match.bins_present_ = false;
}
//! Score tallies using a tracklength estimate of the flux.
//
//! This is triggered at every event (surface crossing, lattice crossing, or