From f14957940dad27da3258ce0b11cdfa3e2ba666c2 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 3 Nov 2017 21:45:49 -0500
Subject: [PATCH] Add openmc.capi.global_tallies()

---
 openmc/capi/tally.py         | 33 +++++++++++++++++++++++++++++----
 src/tallies/tally_header.F90 | 15 ++++++++++++++-
 2 files changed, 43 insertions(+), 5 deletions(-)

diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py
index 7037484ae..14f9d1c1c 100644
--- a/openmc/capi/tally.py
+++ b/openmc/capi/tally.py
@@ -2,6 +2,7 @@ from collections import Mapping
 from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
 from weakref import WeakValueDictionary
 
+import numpy as np
 from numpy.ctypeslib import as_array
 
 from openmc.data.reaction import REACTION_NAME
@@ -11,15 +12,18 @@ from .error import _error_handler, AllocationError, InvalidIDError
 from .filter import _get_filter
 
 
-__all__ = ['Tally', 'tallies']
+__all__ = ['Tally', 'tallies', 'global_tallies']
 
 # Tally functions
-_dll.openmc_get_tally_index.argtypes = [c_int32, POINTER(c_int32)]
-_dll.openmc_get_tally_index.restype = c_int
-_dll.openmc_get_tally_index.errcheck = _error_handler
 _dll.openmc_extend_tallies.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)]
 _dll.openmc_extend_tallies.restype = c_int
 _dll.openmc_extend_tallies.errcheck = _error_handler
+_dll.openmc_get_tally_index.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_get_tally_index.restype = c_int
+_dll.openmc_get_tally_index.errcheck = _error_handler
+_dll.openmc_global_tallies.argtypes = [POINTER(POINTER(c_double))]
+_dll.openmc_global_tallies.restype = c_int
+_dll.openmc_global_tallies.errcheck = _error_handler
 _dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_tally_get_id.restype = c_int
 _dll.openmc_tally_get_id.errcheck = _error_handler
@@ -65,6 +69,27 @@ _SCORES = {
 }
 
 
+def global_tallies():
+    ptr = POINTER(c_double)()
+    _dll.openmc_global_tallies(ptr)
+    array = as_array(ptr, (4, 3))
+
+    # Get sum, sum-of-squares, and number of realizations
+    sum_ = array[:, 1]
+    sum_sq = array[:, 2]
+    n = c_int32.in_dll(_dll, 'n_realizations').value
+
+    # Determine mean
+    mean = sum_ / n
+
+    # Determine standard deviation
+    nonzero = np.abs(mean) > 0
+    stdev = np.zeros_like(mean)
+    stdev[nonzero] = np.sqrt((sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1))
+
+    return list(zip(mean, stdev))
+
+
 class Tally(_FortranObjectWithID):
     """Tally stored internally.
 
diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90
index e137a1327..5aa16c2d0 100644
--- a/src/tallies/tally_header.F90
+++ b/src/tallies/tally_header.F90
@@ -146,7 +146,7 @@ module tally_header
   type(VectorInt), public :: active_surface_tallies
 
   ! Normalization for statistics
-  integer, public :: n_realizations = 0 ! # of independent realizations
+  integer(C_INT32_T), public, bind(C) :: n_realizations = 0 ! # of independent realizations
   real(8), public :: total_weight       ! total starting particle weight in realization
 
 contains
@@ -470,6 +470,19 @@ contains
   end function openmc_get_tally_index
 
 
+  function openmc_global_tallies(ptr) result(err) bind(C)
+    type(C_PTR), intent(out) :: ptr
+    integer(C_INT) :: err
+
+    if (.not. allocated(global_tallies)) then
+      err = E_ALLOCATE
+    else
+      err = 0
+      ptr = C_LOC(global_tallies)
+    end if
+  end function openmc_global_tallies
+
+
   function openmc_tally_get_id(index, id) result(err) bind(C)
     ! Return the ID of a tally
     integer(C_INT32_T), value       :: index