Make capi.keff() usable during inactive and active batches

openmc/capi/core.py

@@ -166,9 +166,17 @@ def keff():
         Mean k-eigenvalue and standard deviation of the mean
 
     """
-    k = (c_double*2)()
-    _dll.openmc_get_keff(k)
-    return tuple(k)
+    n = openmc.capi.num_realizations()
+    if n > 3:
+        # Use the combined estimator if there are enough realizations
+        k = (c_double*2)()
+        _dll.openmc_get_keff(k)
+        return tuple(k)
+    else:
+        # Otherwise, return the tracklength estimator
+        mean = c_double.in_dll(_dll, 'keff').value
+        std_dev = c_double.in_dll(_dll, 'keff_std').value if n > 1 else np.inf
+        return (mean, std_dev)
 
 
 def next_batch():

openmc/capi/error.py

@@ -4,39 +4,39 @@ from warnings import warn
 from . import _dll
 
 
-class Error(Exception):
+class OpenMCError(Exception):
     """Root exception class for OpenMC."""
 
 
-class GeometryError(Error):
+class GeometryError(OpenMCError):
     """Geometry-related error"""
 
 
-class InvalidIDError(Error):
+class InvalidIDError(OpenMCError):
     """Use of an ID that is invalid."""
 
 
-class AllocationError(Error):
+class AllocationError(OpenMCError):
     """Error related to memory allocation."""
 
 
-class OutOfBoundsError(Error):
+class OutOfBoundsError(OpenMCError):
     """Index in array out of bounds."""
 
 
-class DataError(Error):
+class DataError(OpenMCError):
     """Error relating to nuclear data."""
 
 
-class PhysicsError(Error):
+class PhysicsError(OpenMCError):
     """Error relating to performing physics."""
 
 
-class InvalidArgumentError(Error):
+class InvalidArgumentError(OpenMCError):
     """Argument passed was invalid."""
 
 
-class InvalidTypeError(Error):
+class InvalidTypeError(OpenMCError):
     """Tried to perform an operation on the wrong type."""
 
 
@@ -71,4 +71,4 @@ def _error_handler(err, func, args):
     elif err == errcode('e_warning'):
         warn(msg)
     elif err < 0:
-        raise Exception("Unknown error encountered (code {}).".format(err))
+        raise OpenMCError("Unknown error encountered (code {}).".format(err))

openmc/capi/tally.py

@@ -92,12 +92,17 @@ def global_tallies():
     n = num_realizations()
 
     # Determine mean
-    mean = sum_ / n
+    if n > 0:
+        mean = sum_ / n
+    else:
+        mean = sum_.copy()
 
     # 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))
+    stdev = np.empty_like(mean)
+    stdev.fill(np.inf)
+    if n > 1:
+        stdev[nonzero] = np.sqrt((sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1))
 
     return list(zip(mean, stdev))
 
@@ -265,7 +270,7 @@ class Tally(_FortranObjectWithID):
     def std_dev(self):
         results = self.results
         std_dev = np.empty(results.shape[:2])
-        std_dev.fill(np.nan)
+        std_dev.fill(np.inf)
 
         n = self.num_realizations
         if n > 1:

src/simulation_header.F90

@@ -1,5 +1,7 @@
 module simulation_header
 
+  use, intrinsic :: ISO_C_BINDING
+
   use bank_header
   use constants
   use settings, only: gen_per_batch
@@ -33,8 +35,8 @@ module simulation_header
 
   ! Temporary k-effective values
   type(VectorReal) :: k_generation ! single-generation estimates of k
-  real(8) :: keff = ONE       ! average k over active batches
-  real(8) :: keff_std         ! standard deviation of average k
+  real(C_DOUBLE), bind(C) :: keff = ONE  ! average k over active batches
+  real(C_DOUBLE), bind(C) :: keff_std    ! standard deviation of average k
   real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
   real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
   real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength