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

Remove some tally-related C++/Fortran glue-code

Browse source 
This commit also removes some tally trigger special-case code for meshsurface current tallies. That code was broken and is also no longer needed.
This commit is contained in:
Sterling Harper 2019-01-31 13:39:33 -05:00
parent a90d8f3e86
commit a07ea809f2
16 changed files with 106 additions and 494 deletions
Download patch file Download diff file Expand all files Collapse all files

3
CMakeLists.txt
View file

@ -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

3
include/openmc/tallies/tally.h
View file

@ -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.

36
include/openmc/tallies/trigger.h Normal file
View file

@ -0,0 +1,36 @@
#ifndef OPENMC_TALLIES_TRIGGER_H
#define OPENMC_TALLIES_TRIGGER_H
#include <string>
#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

4
src/input_xml.F90
View file

@ -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)

2
src/output.F90
View file

@ -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

8
src/settings.cpp
View file

@ -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
//==============================================================================

24
src/tallies/filter.cpp
View file

@ -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];}

24
src/tallies/tally.cpp
View file

@ -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<const EnergyFilter*>(filt)) {
if (dynamic_cast<const EnergyoutFilter*>(filt)) {
energyout_filter_ = i + 1;
} else {
energyin_filter_ = i + 1;
}
if (dynamic_cast<const EnergyoutFilter*>(filt)) {
energyout_filter_ = i + 1;
} else if (dynamic_cast<const DelayedGroupFilter*>(filt)) {
delayedgroup_filter_ = i + 1;
} else if (dynamic_cast<const SurfaceFilter*>(filt)) {
surface_filter_ = i + 1;
} else if (dynamic_cast<const MeshFilter*>(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;}

2
src/tallies/tally_filter.F90
View file

@ -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

90
src/tallies/tally_filter_header.F90
View file

@ -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

38
src/tallies/tally_filter_mesh.F90
View file

@ -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

22
src/tallies/tally_filter_meshsurface.F90
View file

@ -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

48
src/tallies/tally_header.F90
View file

@ -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

250
src/tallies/trigger.F90
View file

@ -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

32
src/tallies/trigger.cpp Normal file
View file

@ -0,0 +1,32 @@
#include "openmc/tallies/trigger.h"
#include <cmath>
#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

14
src/tallies/trigger_header.F90
View file

@ -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
!===============================================================================