From a07ea809f296b289a373811cbd2afa9bf323fd72 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 31 Jan 2019 13:39:33 -0500 Subject: [PATCH] Remove some tally-related C++/Fortran glue-code This commit also removes some tally trigger special-case code for meshsurface current tallies. That code was broken and is also no longer needed. --- CMakeLists.txt | 3 +- include/openmc/tallies/tally.h | 3 - include/openmc/tallies/trigger.h | 36 ++++ src/input_xml.F90 | 4 +- src/output.F90 | 2 - src/settings.cpp | 8 +- src/tallies/filter.cpp | 24 --- src/tallies/tally.cpp | 24 +-- src/tallies/tally_filter.F90 | 2 - src/tallies/tally_filter_header.F90 | 90 +------- src/tallies/tally_filter_mesh.F90 | 38 ---- src/tallies/tally_filter_meshsurface.F90 | 22 -- src/tallies/tally_header.F90 | 48 ----- src/tallies/trigger.F90 | 250 ++--------------------- src/tallies/trigger.cpp | 32 +++ src/tallies/trigger_header.F90 | 14 +- 16 files changed, 106 insertions(+), 494 deletions(-) create mode 100644 include/openmc/tallies/trigger.h delete mode 100644 src/tallies/tally_filter_mesh.F90 delete mode 100644 src/tallies/tally_filter_meshsurface.F90 create mode 100644 src/tallies/trigger.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 4671b75c0..8fb3f9c0d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -359,8 +359,6 @@ add_library(libopenmc SHARED src/tallies/tally_filter_delayedgroup.F90 src/tallies/tally_filter_distribcell.F90 src/tallies/tally_filter_energy.F90 - src/tallies/tally_filter_mesh.F90 - src/tallies/tally_filter_meshsurface.F90 src/tallies/tally_filter_particle.F90 src/tallies/tally_filter_sph_harm.F90 src/tallies/tally_header.F90 @@ -442,6 +440,7 @@ add_library(libopenmc SHARED src/tallies/filter_universe.cpp src/tallies/filter_zernike.cpp src/tallies/tally.cpp + src/tallies/trigger.cpp src/timer.cpp src/thermal.cpp src/wmp.cpp diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index a41c82b61..d4bbfbaa1 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -56,11 +56,8 @@ public: // We need to have quick access to some filters. The following gives indices // for various filters that could be in the tally or C_NONE if they are not // present. - int energyin_filter_ {C_NONE}; int energyout_filter_ {C_NONE}; int delayedgroup_filter_ {C_NONE}; - int surface_filter_ {C_NONE}; - int mesh_filter_ {C_NONE}; int deriv_ {C_NONE}; //!< Index of a TallyDerivative object for diff tallies. diff --git a/include/openmc/tallies/trigger.h b/include/openmc/tallies/trigger.h new file mode 100644 index 000000000..56c15f7bc --- /dev/null +++ b/include/openmc/tallies/trigger.h @@ -0,0 +1,36 @@ +#ifndef OPENMC_TALLIES_TRIGGER_H +#define OPENMC_TALLIES_TRIGGER_H + +#include + +#include "pugixml.hpp" + +namespace openmc { + +//! Stops the simulation early if a desired tally uncertainty is reached. + +extern "C" struct Trigger +{ + int type; //!< variance, std_dev, or rel_err + double threshold; //!< uncertainty value below which trigger is satisfied + int score_index; //!< index of the relevant score in the tally's arrays + double variance {0.}; + double std_dev {0.}; + double rel_err {0.}; +}; + +//! Stops the simulation early if a desired k-effective uncertainty is reached. + +extern "C" struct KTrigger +{ + int type{0}; + double threshold {0.}; +}; + +//TODO: consider a different namespace +namespace settings { + extern "C" KTrigger keff_trigger; +} + +} // namespace openmc +#endif // OPENMC_TALLIES_TRIGGER_H diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c657b9722..160acd263 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1152,9 +1152,9 @@ contains has_energyout = (t % energyout_filter() > 0) has_delayedgroup = (t % delayedgroup_filter() > 0) has_legendre = .false. - has_surface = (t % surface_filter() > 0) has_cell = .false. has_cellfrom = .false. + has_surface = .false. has_meshsurface = .false. particle_filter_index = 0 do j = 1, t % n_filters() @@ -1165,6 +1165,8 @@ contains has_cell = .true. type is (CellfromFilter) has_cellfrom = .true. + type is (SurfaceFilter) + has_surface = .true. type is (MeshSurfaceFilter) has_meshsurface = .true. type is (ParticleFilter) diff --git a/src/output.F90 b/src/output.F90 index c4e0c3477..5895cd59c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -20,8 +20,6 @@ module output use tally_header use tally_derivative_header use tally_filter - use tally_filter_mesh, only: MeshFilter - use tally_filter_header, only: TallyFilterMatch use timer_header implicit none diff --git a/src/settings.cpp b/src/settings.cpp index dbbe9c6da..a7543d5e6 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -25,6 +25,7 @@ #include "openmc/simulation.h" #include "openmc/source.h" #include "openmc/string_utils.h" +#include "openmc/tallies/trigger.h" #include "openmc/xml_interface.h" namespace openmc { @@ -106,13 +107,6 @@ int verbosity {7}; double weight_cutoff {0.25}; double weight_survive {1.0}; -// TODO: Move to separate file -struct KTrigger { - int type; - double threshold; -}; -extern "C" KTrigger keff_trigger; - } // namespace settings //============================================================================== diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 7c641f854..f331e6312 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -54,30 +54,6 @@ extern "C" { filter_match_get_i_bin(FilterMatch* match) {return match->i_bin_;} - void - filter_match_set_i_bin(FilterMatch* match, int i) - {match->i_bin_ = i;} - - void - filter_match_bins_push_back(FilterMatch* match, int val) - {match->bins_.push_back(val);} - - void - filter_match_weights_push_back(FilterMatch* match, double val) - {match->weights_.push_back(val);} - - void - filter_match_bins_clear(FilterMatch* match) - {match->bins_.clear();} - - void - filter_match_weights_clear(FilterMatch* match) - {match->weights_.clear();} - - int - filter_match_bins_size(FilterMatch* match) - {return match->bins_.size();} - int filter_match_bins_data(FilterMatch* match, int indx) {return match->bins_[indx-1];} diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 104885edd..835747253 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -86,18 +86,10 @@ Tally::set_filters(const int32_t filter_indices[], int n) const auto* filt = model::tally_filters[i_filt-1].get(); //TODO: off-by-one on each index - if (dynamic_cast(filt)) { - if (dynamic_cast(filt)) { - energyout_filter_ = i + 1; - } else { - energyin_filter_ = i + 1; - } + if (dynamic_cast(filt)) { + energyout_filter_ = i + 1; } else if (dynamic_cast(filt)) { delayedgroup_filter_ = i + 1; - } else if (dynamic_cast(filt)) { - surface_filter_ = i + 1; - } else if (dynamic_cast(filt)) { - mesh_filter_ = i + 1; } } @@ -1021,9 +1013,6 @@ extern "C" { int active_tallies_size() {return model::active_tallies.size();} - int active_analog_tallies_data(int i) - {return model::active_analog_tallies[i-1];} - int active_analog_tallies_size() {return model::active_analog_tallies.size();} @@ -1075,21 +1064,12 @@ extern "C" { bool tally_get_all_nuclides_c(Tally* tally) {return tally->all_nuclides_;} - int tally_get_energyin_filter_c(Tally* tally) - {return tally->energyin_filter_;} - int tally_get_energyout_filter_c(Tally* tally) {return tally->energyout_filter_;} int tally_get_delayedgroup_filter_c(Tally* tally) {return tally->delayedgroup_filter_;} - int tally_get_surface_filter_c(Tally* tally) - {return tally->surface_filter_;} - - int tally_get_mesh_filter_c(Tally* tally) - {return tally->mesh_filter_;} - int tally_get_deriv_c(Tally* tally) {return tally->deriv_;} int tally_set_deriv_c(Tally* tally, int deriv) {tally->deriv_ = deriv;} diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index a1fc01dab..231fe5525 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -10,8 +10,6 @@ module tally_filter use tally_filter_delayedgroup use tally_filter_distribcell use tally_filter_energy - use tally_filter_mesh - use tally_filter_meshsurface use tally_filter_particle use tally_filter_sph_harm diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index 5a8784ae0..2e1e7d660 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -35,18 +35,8 @@ module tally_filter_header type, public :: TallyFilterMatch type(C_PTR) :: ptr - - ! Indicates whether all valid bins for this filter have been found - logical :: bins_present = .false. - contains procedure :: i_bin - procedure :: set_i_bin - procedure :: bins_push_back - procedure :: weights_push_back - procedure :: bins_clear - procedure :: weights_clear - procedure :: bins_size procedure :: bins_data procedure :: weights_data procedure :: bins_set_data @@ -98,6 +88,12 @@ module tally_filter_header type, public, extends(TallyFilter) :: MaterialFilter end type + type, public, extends(TallyFilter) :: MeshFilter + end type + + type, public, extends(TallyFilter) :: MeshSurfaceFilter + end type + type, public, extends(TallyFilter) :: MuFilter end type @@ -162,80 +158,6 @@ contains i = filter_match_get_i_bin(this % ptr) end function - subroutine set_i_bin(this, i) - class(TallyFilterMatch) :: this - integer :: i - interface - subroutine filter_match_set_i_bin(ptr, i) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT), value :: i - end subroutine - end interface - call filter_match_set_i_bin(this % ptr, i) - end subroutine - - subroutine bins_push_back(this, val) - class(TallyFilterMatch), intent(inout) :: this - integer, intent(in) :: val - interface - subroutine filter_match_bins_push_back(ptr, val) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT), intent(in), value :: val - end subroutine - end interface - call filter_match_bins_push_back(this % ptr, val) - end subroutine bins_push_back - - subroutine weights_push_back(this, val) - class(TallyFilterMatch), intent(inout) :: this - real(8), intent(in) :: val - interface - subroutine filter_match_weights_push_back(ptr, val) bind(C) - import C_PTR, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE), intent(in), value :: val - end subroutine - end interface - call filter_match_weights_push_back(this % ptr, val) - end subroutine weights_push_back - - subroutine bins_clear(this) - class(TallyFilterMatch), intent(inout) :: this - interface - subroutine filter_match_bins_clear(ptr) bind(C) - import C_PTR - type(C_PTR), value :: ptr - end subroutine - end interface - call filter_match_bins_clear(this % ptr) - end subroutine bins_clear - - subroutine weights_clear(this) - class(TallyFilterMatch), intent(inout) :: this - interface - subroutine filter_match_weights_clear(ptr) bind(C) - import C_PTR - type(C_PTR), value :: ptr - end subroutine - end interface - call filter_match_weights_clear(this % ptr) - end subroutine weights_clear - - function bins_size(this) result(len) - class(TallyFilterMatch), intent(inout) :: this - integer :: len - interface - function filter_match_bins_size(ptr) bind(C) result(len) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT) :: len - end function - end interface - len = filter_match_bins_size(this % ptr) - end function bins_size - function bins_data(this, indx) result(val) class(TallyFilterMatch), intent(inout) :: this integer, intent(in) :: indx diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 deleted file mode 100644 index 626ed4b15..000000000 --- a/src/tallies/tally_filter_mesh.F90 +++ /dev/null @@ -1,38 +0,0 @@ -module tally_filter_mesh - - use, intrinsic :: ISO_C_BINDING - - use tally_filter_header - - implicit none - - interface - function openmc_mesh_filter_set_mesh(index, index_mesh) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: index_mesh - integer(C_INT) :: err - end function - end interface - - type, public, extends(TallyFilter) :: MeshFilter - contains - procedure :: mesh => get_mesh - end type MeshFilter - -contains - - function get_mesh(this) result(mesh) - class(MeshFilter) :: this - integer :: mesh - interface - function mesh_filter_get_mesh(filt) result(index_mesh) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: filt - integer(C_INT) :: index_mesh - end function mesh_filter_get_mesh - end interface - mesh = mesh_filter_get_mesh(this % ptr) - end function get_mesh - -end module tally_filter_mesh diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 deleted file mode 100644 index 574d2c8b7..000000000 --- a/src/tallies/tally_filter_meshsurface.F90 +++ /dev/null @@ -1,22 +0,0 @@ -module tally_filter_meshsurface - - use, intrinsic :: ISO_C_BINDING - - use tally_filter_header - - implicit none - - interface - function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) & - bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: index_mesh - integer(C_INT) :: err - end function - end interface - - type, extends(TallyFilter) :: MeshSurfaceFilter - end type MeshSurfaceFilter - -end module tally_filter_meshsurface diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index b8be99424..3e7610027 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -64,12 +64,6 @@ module tally_header integer(C_INT) :: size end function - function active_analog_tallies_data(i) result(tally) bind(C) - import C_INT - integer(C_INT), value :: i - integer(C_INT) :: tally - end function - function active_tracklength_tallies_size() result(size) bind(C) import C_INT integer(C_INT) :: size @@ -143,11 +137,8 @@ module tally_header procedure :: nuclide_bins => tally_get_nuclide_bins procedure :: set_nuclide_bins => tally_set_nuclide_bins procedure :: all_nuclides => tally_get_all_nuclides - procedure :: energyin_filter => tally_get_energyin_filter procedure :: energyout_filter => tally_get_energyout_filter procedure :: delayedgroup_filter => tally_get_delayedgroup_filter - procedure :: surface_filter => tally_get_surface_filter - procedure :: mesh_filter => tally_get_mesh_filter procedure :: deriv => tally_get_deriv procedure :: set_deriv => tally_set_deriv end type TallyObject @@ -484,19 +475,6 @@ contains all_nuc = tally_get_all_nuclides_c(this % ptr) end function - function tally_get_energyin_filter(this) result(filt) - class(TallyObject) :: this - integer(C_INT) :: filt - interface - function tally_get_energyin_filter_c(tally) result(filt) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: tally - integer(C_INT) :: filt - end function - end interface - filt = tally_get_energyin_filter_c(this % ptr) - end function - function tally_get_energyout_filter(this) result(filt) class(TallyObject) :: this integer(C_INT) :: filt @@ -523,32 +501,6 @@ contains filt = tally_get_delayedgroup_filter_c(this % ptr) end function - function tally_get_surface_filter(this) result(filt) - class(TallyObject) :: this - integer(C_INT) :: filt - interface - function tally_get_surface_filter_c(tally) result(filt) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: tally - integer(C_INT) :: filt - end function - end interface - filt = tally_get_surface_filter_c(this % ptr) - end function - - function tally_get_mesh_filter(this) result(filt) - class(TallyObject) :: this - integer(C_INT) :: filt - interface - function tally_get_mesh_filter_c(tally) result(filt) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: tally - integer(C_INT) :: filt - end function - end interface - filt = tally_get_mesh_filter_c(this % ptr) - end function - function tally_get_deriv(this) result(deriv) class(TallyObject) :: this integer(C_INT) :: deriv diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 785e7b97b..93068f4f5 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -2,27 +2,32 @@ module trigger use, intrinsic :: ISO_C_BINDING -#ifdef OPENMC_MPI - use message_passing -#endif - use constants use eigenvalue, only: openmc_get_keff use endf, only: reaction_name use error, only: warning, write_message use string, only: to_str - use mesh_header, only: RegularMesh, meshes use message_passing, only: master use settings use simulation_header use trigger_header use tally, only: TallyObject - use tally_filter_mesh, only: MeshFilter - use tally_filter_header, only: filter_matches, filters use tally_header, only: tallies, n_tallies implicit none + interface + subroutine get_tally_uncertainty(i_tally, score_index, filter_index, & + std_dev, rel_err) bind(C) + import C_INT, C_DOUBLE + integer(C_INT), value :: i_tally + integer(C_INT), value :: score_index + integer(C_INT), value :: filter_index + real(C_DOUBLE) :: std_dev + real(C_DOUBLE) :: rel_err + end subroutine + end interface + contains !=============================================================================== @@ -171,22 +176,10 @@ contains TRIGGER_LOOP: do s = 1, t % n_triggers associate (trigger => t % triggers(s)) - ! Initialize trigger uncertainties to zero - trigger % std_dev = ZERO - trigger % rel_err = ZERO - trigger % variance = ZERO - - ! Mesh current tally triggers require special treatment - if (t % type() == TALLY_MESH_SURFACE) then - call compute_tally_current(t, trigger) - - else - - ! Initialize bins, filter level - do j = 1, t % n_filters() - call filter_matches(t % filter(j)) % bins_clear() - call filter_matches(t % filter(j)) % bins_push_back(0) - end do + ! Initialize trigger uncertainties to zero + trigger % std_dev = ZERO + trigger % rel_err = ZERO + trigger % variance = ZERO FILTER_LOOP: do filter_index = 1, t % n_filter_bins() @@ -196,8 +189,8 @@ contains ! Initialize score bin index NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins() - call get_trigger_uncertainty(std_dev, rel_err, & - score_index, filter_index, t) + call get_tally_uncertainty(i, score_index, filter_index, & + std_dev, rel_err) if (trigger % variance < variance) then trigger % variance = std_dev ** 2 @@ -236,217 +229,10 @@ contains end do NUCLIDE_LOOP if (t % n_filters() == 0) exit FILTER_LOOP end do FILTER_LOOP - end if end associate end do TRIGGER_LOOP end associate end do TALLY_LOOP end subroutine check_tally_triggers -!=============================================================================== -! COMPUTE_TALLY_CURRENT computes the current for a mesh current tally with -! precision trigger(s). -!=============================================================================== - - subroutine compute_tally_current(t, trigger) - type(TallyObject), intent(in) :: t ! mesh current tally - type(TriggerObject), intent(inout) :: trigger ! mesh current tally trigger - - integer :: i ! mesh index - integer :: j ! loop index for tally filters - integer :: ijk(3) ! indices of mesh cells - integer :: n_dim ! number of mesh dimensions - integer :: n_cells ! number of mesh cells - integer :: l ! index for energy - integer :: i_filter_mesh ! index for mesh filter - integer :: i_filter_ein ! index for incoming energy filter - integer :: i_filter_surf ! index for surface filter - integer :: n ! number of incoming energy bins - integer :: filter_index ! index in results array for filters - logical :: print_ebin ! should incoming energy bin be displayed? - real(8) :: rel_err = ZERO ! temporary relative error of result - real(8) :: std_dev = ZERO ! temporary standard deviration of result - type(RegularMesh) :: m ! surface current mesh - - ! Get pointer to mesh - i_filter_mesh = t % mesh_filter() - i_filter_surf = t % surface_filter() - select type(filt => filters(i_filter_mesh) % obj) - type is (MeshFilter) - m = meshes(filt % mesh()) - end select - - ! initialize bins array - do j = 1, t % n_filters() - call filter_matches(t % filter(j)) % bins_clear() - call filter_matches(t % filter(j)) % bins_push_back(1) - end do - - ! determine how many energyin bins there are - i_filter_ein = t % energyin_filter() - if (i_filter_ein > 0) then - print_ebin = .true. - n = filters(t % filter(i_filter_ein)) % obj % n_bins - i_filter_ein = t % filter(i_filter_ein) - else - print_ebin = .false. - n = 1 - end if - - ! Get the dimensions and number of cells in the mesh - n_dim = m % n_dimension() - n_cells = 1 - do j = 1, n_dim - n_cells = n_cells * m % dimension(j) - end do - - ! Loop over all the mesh cells - do i = 1, n_cells - - ! Get the indices for this cell - call m % get_indices_from_bin(i, ijk) - call filter_matches(i_filter_mesh) % bins_set_data(1, i) - - do l = 1, n - - if (print_ebin) then - call filter_matches(i_filter_ein) % bins_set_data(1, l) - end if - - ! Left Surface - call filter_matches(i_filter_surf) % bins_set_data(1, OUT_LEFT) - filter_index = 1 - do j = 1, t % n_filters() - filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins_data(1) - 1) * t % stride(j) - end do - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = std_dev**2 - - ! Right Surface - call filter_matches(i_filter_surf) % bins_set_data(1, OUT_RIGHT) - filter_index = 1 - do j = 1, t % n_filters() - filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins_data(1) - 1) * t % stride(j) - end do - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - - ! Back Surface - call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BACK) - filter_index = 1 - do j = 1, t % n_filters() - filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins_data(1) - 1) * t % stride(j) - end do - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - - ! Front Surface - call filter_matches(i_filter_surf) % bins_set_data(1, OUT_FRONT) - filter_index = 1 - do j = 1, t % n_filters() - filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins_data(1) - 1) * t % stride(j) - end do - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - - ! Bottom Surface - call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BOTTOM) - filter_index = 1 - do j = 1, t % n_filters() - filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins_data(1) - 1) * t % stride(j) - end do - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - - ! Top Surface - call filter_matches(i_filter_surf) % bins_set_data(1, OUT_TOP) - filter_index = 1 - do j = 1, t % n_filters() - filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins_data(1) - 1) * t % stride(j) - end do - call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) - if (trigger % std_dev < std_dev) then - trigger % std_dev = std_dev - end if - if (trigger % rel_err < rel_err) then - trigger % rel_err = rel_err - end if - trigger % variance = trigger % std_dev**2 - - end do - end do - - end subroutine compute_tally_current - -!=============================================================================== -! GET_TRIGGER_UNCERTAINTY computes the standard deviation and relative error -! for a single tally bin for CHECK_TALLY_TRIGGERS. -!=============================================================================== - - subroutine get_trigger_uncertainty(std_dev, rel_err, score_index, & - filter_index, t) - - real(8), intent(inout) :: std_dev ! tally standard deviation - real(8), intent(inout) :: rel_err ! tally relative error - integer, intent(in) :: score_index ! tally results score index - integer, intent(in) :: filter_index ! tally results filter index - type(TallyObject), intent(in) :: t ! tally - - integer :: n ! number of realizations - real(8) :: mean ! tally mean - real(8) :: tally_sum, tally_sum_sq ! results for a single tally bin - - n = t % n_realizations - tally_sum = t % results(RESULT_SUM, score_index, filter_index) - tally_sum_sq = t % results(RESULT_SUM_SQ, score_index, filter_index) - - ! Compute the tally mean and standard deviation - mean = tally_sum / n - std_dev = sqrt((tally_sum_sq / n - mean * mean) / (n - 1)) - - ! Compute the relative error if the mean is non-zero - if (mean == ZERO) then - rel_err = ZERO - else - rel_err = std_dev / mean - end if - - end subroutine get_trigger_uncertainty - end module trigger diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp new file mode 100644 index 000000000..d3069d9a0 --- /dev/null +++ b/src/tallies/trigger.cpp @@ -0,0 +1,32 @@ +#include "openmc/tallies/trigger.h" + +#include + +#include "openmc/capi.h" +#include "openmc/tallies/tally.h" + +namespace openmc { + +extern "C" void +get_tally_uncertainty(int i_tally, int score_index, int filter_index, + double* std_dev, double* rel_err) +{ + int n; + int err = openmc_tally_get_n_realizations(i_tally, &n); + + auto results = tally_results(i_tally); + //TODO: off-by-one + auto sum = results(filter_index-1, score_index-1, RESULT_SUM); + auto sum_sq = results(filter_index-1, score_index-1, RESULT_SUM_SQ); + + auto mean = sum / n; + *std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1)); + + if (mean > 0) { + *rel_err = *std_dev / mean; + } else { + *rel_err = 0.; + } +} + +} // namespace openmc diff --git a/src/tallies/trigger_header.F90 b/src/tallies/trigger_header.F90 index ccf497d1d..c370483ce 100644 --- a/src/tallies/trigger_header.F90 +++ b/src/tallies/trigger_header.F90 @@ -11,13 +11,13 @@ module trigger_header ! TRIGGEROBJECT stores the variance, relative error and standard deviation ! for some user-specified trigger. !=============================================================================== - type, public :: TriggerObject - integer :: type ! "variance", "std_dev" or "rel_err" - real(8) :: threshold ! a convergence threshold - integer :: score_index ! the index of the score - real(8) :: variance = ZERO ! temp variance container - real(8) :: std_dev = ZERO ! temp std. dev. container - real(8) :: rel_err = ZERO ! temp rel. err. container + type, public, bind(C) :: TriggerObject + integer(C_INT) :: type ! "variance", "std_dev" or "rel_err" + real(C_DOUBLE) :: threshold ! a convergence threshold + integer(C_INT) :: score_index ! the index of the score + real(C_DOUBLE) :: variance = ZERO ! temp variance container + real(C_DOUBLE) :: std_dev = ZERO ! temp std. dev. container + real(C_DOUBLE) :: rel_err = ZERO ! temp rel. err. container end type TriggerObject !===============================================================================