Merge pull request #1270 from drewejohnson/dep-timings

Store time required to deplete materials in depletion results

docs/source/io_formats/depletion_results.rst

@@ -25,6 +25,9 @@ The current version of the depletion results file format is 1.0.
              nuclides, number of reactions).
            - **time** (*double[][2]*) -- Time in [s] at beginning/end of each
              step.
+           - **depletion time** (*double[]*) -- Average process time in [s] 
+             spent depleting a material across all burnable materials and,
+             if applicable, MPI processes.
 
 **/materials/<id>/**
 

openmc/deplete/__init__.py

@@ -8,6 +8,8 @@ from sys import exit
 
 from h5py import get_config
 
+from unittest.mock import Mock
+
 from .dummy_comm import DummyCommunicator
 
 try:
@@ -27,6 +29,7 @@ try:
 except ImportError:
     comm = DummyCommunicator()
     have_mpi = False
+    MPI = Mock()
 
 from .nuclide import *
 from .chain import *

openmc/deplete/integrator/cecm.py

@@ -3,7 +3,7 @@
 import copy
 from collections.abc import Iterable
 
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 
 
@@ -94,7 +94,8 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
                 op_results[0].rates *= ratio_power[0]
 
             # Deplete for first half of timestep
-            x_middle = deplete(chain, x[0], op_results[0].rates, dt/2, print_out)
+            proc_time, x_middle = timed_deplete(
+                chain, x[0], op_results[0].rates, dt/2, print_out)
 
             # Get middle-of-timestep reaction rates
             x.append(x_middle)
@@ -102,10 +103,13 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
 
             # Deplete for full timestep using beginning-of-step materials
             # and middle-of-timestep reaction rates
-            x_end = deplete(chain, x[0], op_results[1].rates, dt, print_out)
+            pt_end, x_end = timed_deplete(
+                chain, x[0], op_results[1].rates, dt, print_out)
 
             # Create results, write to disk
-            Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+            Results.save(
+                operator, x, op_results, [t, t + dt], p, i_res + i,
+                proc_time + pt_end)
 
             # Advance time, update vector
             t += dt

openmc/deplete/integrator/celi.py

@@ -3,7 +3,7 @@
 import copy
 from collections.abc import Iterable
 
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 
 
@@ -149,19 +149,19 @@ def celi_inner(operator, vec, p, i, i_res, t, dt, print_out):
         op_results[0].rates *= ratio_power[0]
 
     # Deplete to end
-    x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out)
+    proc_time, x_new = timed_deplete(chain, x[0], op_results[0].rates, dt, print_out)
     x.append(x_new)
     op_results.append(operator(x[1], p))
 
     # Deplete with two matrix exponentials
     rates = list(zip(op_results[0].rates, op_results[1].rates))
-    x_end = deplete(chain, x[0], rates, dt, print_out,
+    time_1, x_end = timed_deplete(chain, x[0], rates, dt, print_out,
                     matrix_func=_celi_f1)
-    x_end = deplete(chain, x_end, rates, dt, print_out,
+    time_2, x_end = timed_deplete(chain, x_end, rates, dt, print_out,
                     matrix_func=_celi_f2)
 
     # Create results, write to disk
-    Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+    Results.save(operator, x, op_results, [t, t + dt], p, i_res + i, proc_time + time_1 + time_2)
 
     # return updated time and vectors
     return x_end, t + dt, op_results[0]

openmc/deplete/integrator/cf4.py

@@ -3,7 +3,7 @@
 import copy
 from collections.abc import Iterable
 
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 
 
@@ -119,34 +119,38 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True):
                 op_results[0].rates *= ratio_power[0]
 
             # Step 1: deplete with matrix 1/2*A(y0)
-            x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out,
-                            matrix_func=_cf4_f1)
+            time_1, x_new = timed_deplete(
+                chain, x[0], op_results[0].rates, dt, print_out,
+                matrix_func=_cf4_f1)
             x.append(x_new)
             op_results.append(operator(x_new, p))
 
             # Step 2: deplete with matrix 1/2*A(y1)
-            x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out,
-                            matrix_func=_cf4_f1)
+            time_2, x_new = timed_deplete(
+                chain, x[0], op_results[1].rates, dt, print_out,
+                matrix_func=_cf4_f1)
             x.append(x_new)
             op_results.append(operator(x_new, p))
 
             # Step 3: deplete with matrix -1/2*A(y0)+A(y2)
             rates = list(zip(op_results[0].rates, op_results[2].rates))
-            x_new = deplete(chain, x[1], rates, dt, print_out,
-                            matrix_func=_cf4_f2)
+            time_3, x_new = timed_deplete(
+                chain, x[1], rates, dt, print_out, matrix_func=_cf4_f2)
             x.append(x_new)
             op_results.append(operator(x_new, p))
 
             # Step 4: deplete with two matrix exponentials
             rates = list(zip(op_results[0].rates, op_results[1].rates,
                              op_results[2].rates, op_results[3].rates))
-            x_end = deplete(chain, x[0], rates, dt, print_out,
-                            matrix_func=_cf4_f3)
-            x_end = deplete(chain, x_end, rates, dt, print_out,
-                            matrix_func=_cf4_f4)
+            time_4, x_end = timed_deplete(
+                chain, x[0], rates, dt, print_out, matrix_func=_cf4_f3)
+            time_5, x_end = timed_deplete(
+                chain, x_end, rates, dt, print_out, matrix_func=_cf4_f4)
 
             # Create results, write to disk
-            Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+            Results.save(
+                operator, x, op_results, [t, t + dt], p, i_res + i,
+                time_1 + time_2 + time_3 + time_4 + time_5)
 
             # Advance time, update vector
             t += dt

openmc/deplete/integrator/cram.py

@@ -53,6 +53,25 @@ def deplete(chain, x, rates, dt, print_out=True, matrix_func=None):
     return x_result
 
 
+def timed_deplete(*args, **kwargs):
+    """Wrapper over :func:`deplete` that also returns process time
+
+    All arguments and keyword arguments are passed onto
+    :func:`deplete` directly.
+
+    Returns
+    -------
+    proc_time: float
+        Process time required to return from deplete
+    results: list of numpy arrays
+        Output from :func:`deplete` call
+    """
+
+    start = time.time()
+    results = deplete(*args, **kwargs)
+    return time.time() - start, results
+
+
 def _cram_wrapper(chain, n0, rates, dt, matrix_func=None):
     """Wraps depletion matrix creation / CRAM solve for multiprocess execution
 

openmc/deplete/integrator/epc_rk4.py

@@ -3,7 +3,7 @@
 import copy
 from collections.abc import Iterable
 
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 
 
@@ -105,30 +105,35 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True)
                 op_results[0].rates *= ratio_power[0]
 
             # Step 1: deplete with matrix 1/2*A(y0)
-            x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out,
-                            matrix_func=_rk4_f1)
+            time_1, x_new = timed_deplete(
+                chain, x[0], op_results[0].rates, dt, print_out,
+                matrix_func=_rk4_f1)
             x.append(x_new)
             op_results.append(operator(x[1], p))
 
             # Step 2: deplete with matrix 1/2*A(y1)
-            x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out,
-                            matrix_func=_rk4_f1)
+            time_2, x_new = timed_deplete(
+                chain, x[0], op_results[1].rates, dt, print_out,
+                matrix_func=_rk4_f1)
             x.append(x_new)
             op_results.append(operator(x[2], p))
 
             # Step 3: deplete with matrix A(y2)
-            x_new = deplete(chain, x[0], op_results[2].rates, dt, print_out)
+            time_3, x_new = timed_deplete(
+                chain, x[0], op_results[2].rates, dt, print_out)
             x.append(x_new)
             op_results.append(operator(x[3], p))
 
             # Step 4: deplete with matrix 1/6*A(y0)+1/3*A(y1)+1/3*A(y2)+1/6*A(y3)
             rates = list(zip(op_results[0].rates, op_results[1].rates,
                              op_results[2].rates, op_results[3].rates))
-            x_end = deplete(chain, x[0], rates, dt, print_out,
-                            matrix_func=_rk4_f4)
+            time_4, x_end = timed_deplete(
+                chain, x[0], rates, dt, print_out, matrix_func=_rk4_f4)
 
             # Create results, write to disk
-            Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+            Results.save(
+                operator, x, op_results, [t, t + dt], p, i_res + i,
+                time_1 + time_2 + time_3 + time_4)
 
             # Advance time, update vector
             t += dt
@@ -139,4 +144,5 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True)
         op_results = [operator(x[0], power[-1])]
 
         # Create results, write to disk
-        Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
+        Results.save(
+            operator, x, op_results, [t, t], p, i_res + len(timesteps))

openmc/deplete/integrator/leqi.py

@@ -5,7 +5,7 @@ from collections.abc import Iterable
 from itertools import repeat
 
 from .celi import celi_inner
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 
 
@@ -137,22 +137,24 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True):
 
             inputs = list(zip(op_res_last.rates, op_results[0].rates,
                               repeat(dt_l), repeat(dt)))
-            x_new = deplete(chain, x[0], inputs, dt, print_out,
-                            matrix_func=_leqi_f1)
-            x_new = deplete(chain, x_new, inputs, dt, print_out,
-                            matrix_func=_leqi_f2)
+            time_1, x_new = timed_deplete(
+                chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f1)
+            time_2, x_new = timed_deplete(
+                chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f2)
             x.append(x_new)
             op_results.append(operator(x[1], p))
 
             inputs = list(zip(op_res_last.rates, op_results[0].rates,
                               op_results[1].rates, repeat(dt_l), repeat(dt)))
-            x_new = deplete(chain, x[0], inputs, dt, print_out,
-                            matrix_func=_leqi_f3)
-            x_new = deplete(chain, x_new, inputs, dt, print_out,
-                            matrix_func=_leqi_f4)
+            time_3, x_new = timed_deplete(
+                chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3)
+            time_4, x_new = timed_deplete(
+                chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4)
 
             # Create results, write to disk
-            Results.save(operator, x, op_results, [t, t+dt], p, i_res+i)
+            Results.save(
+                operator, x, op_results, [t, t+dt], p, i_res+i,
+                time_1 + time_2 + time_3 + time_4)
 
             # update results
             op_res_last = copy.deepcopy(op_results[0])
@@ -164,4 +166,5 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True):
         op_results = [operator(x[0], power[-1])]
 
         # Create results, write to disk
-        Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
+        Results.save(
+            operator, x, op_results, [t, t], p, i_res + len(timesteps))

openmc/deplete/integrator/predictor.py

@@ -3,7 +3,7 @@
 import copy
 from collections.abc import Iterable
 
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 
 
@@ -53,6 +53,8 @@ def predictor(operator, timesteps, power=None, power_density=None,
     if not isinstance(power, Iterable):
         power = [power]*len(timesteps)
 
+    proc_time = None
+
     # Generate initial conditions
     with operator as vec:
         # Initialize time and starting index
@@ -74,7 +76,7 @@ def predictor(operator, timesteps, power=None, power_density=None,
                 op_results = [operator(x[0], p)]
 
                 # Create results, write to disk
-                Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
+                Results.save(operator, x, op_results, [t, t + dt], p, i_res + i, proc_time)
             else:
                 # Get initial concentration
                 x = [operator.prev_res[-1].data[0]]
@@ -89,7 +91,8 @@ def predictor(operator, timesteps, power=None, power_density=None,
                 op_results[0].rates *= ratio_power[0]
 
             # Deplete for full timestep
-            x_end = deplete(chain, x[0], op_results[0].rates, dt, print_out)
+            proc_time, x_end = timed_deplete(
+                chain, x[0], op_results[0].rates, dt, print_out)
 
             # Advance time, update vector
             t += dt
@@ -100,4 +103,4 @@ def predictor(operator, timesteps, power=None, power_density=None,
         op_results = [operator(x[0], power[-1])]
 
         # Create results, write to disk
-        Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
+        Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps), proc_time)

openmc/deplete/integrator/si_celi.py

@@ -3,7 +3,7 @@
 import copy
 from collections.abc import Iterable
 
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 from ..abc import OperatorResult
 from .celi import _celi_f1, _celi_f2
@@ -94,7 +94,8 @@ def si_celi(operator, timesteps, power=None, power_density=None,
                                              i, i_res, t, dt, print_out, m)
 
         # Create results for last point, write to disk
-        Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
+        Results.save(
+            operator, x, op_results, [t, t], p, i_res + len(timesteps))
 
 
 def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out, m=10):
@@ -136,7 +137,8 @@ def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out, m=10):
     chain = operator.chain
 
     # Deplete to end
-    x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out)
+    proc_time, x_new = timed_deplete(
+        chain, x[0], op_results[0].rates, dt, print_out)
     x.append(x_new)
 
     for j in range(m + 1):
@@ -150,14 +152,15 @@ def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out, m=10):
             op_res_bar = OperatorResult(k, rates)
 
         rates = list(zip(op_results[0].rates, op_res_bar.rates))
-        x_new = deplete(chain, x[0], rates, dt, print_out,
-                        matrix_func=_celi_f1)
-        x_new = deplete(chain, x_new, rates, dt, print_out,
-                        matrix_func=_celi_f2)
+        time_1, x_new = timed_deplete(
+            chain, x[0], rates, dt, print_out, matrix_func=_celi_f1)
+        time_2, x_new = timed_deplete(
+            chain, x_new, rates, dt, print_out, matrix_func=_celi_f2)
+        proc_time += time_1 + time_2
 
     # Create results, write to disk
     op_results.append(op_res_bar)
-    Results.save(operator, x, op_results, [t, t+dt], p, i_res+i)
+    Results.save(operator, x, op_results, [t, t+dt], p, i_res+i, proc_time)
 
     # return updated time and vectors
     return [x_new], t + dt, [op_res_bar]

openmc/deplete/integrator/si_leqi.py

@@ -6,7 +6,7 @@ from itertools import repeat
 
 from .si_celi import si_celi_inner
 from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4
-from .cram import deplete
+from .cram import timed_deplete
 from ..results import Results
 from ..abc import OperatorResult
 
@@ -112,12 +112,14 @@ def si_leqi(operator, timesteps, power=None, power_density=None,
             # Perform remaining LE/QI
             inputs = list(zip(op_res_last.rates, op_results[0].rates,
                               repeat(dt_l), repeat(dt)))
-            x_new = deplete(chain, x[0], inputs, dt, print_out,
-                            matrix_func=_leqi_f1)
-            x_new = deplete(chain, x_new, inputs, dt, print_out,
-                            matrix_func=_leqi_f2)
+            proc_time, x_new = timed_deplete(
+                chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f1)
+            time_1, x_new = timed_deplete(
+                chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f2)
             x.append(x_new)
 
+            proc_time += time_1
+
             # Loop on inner
             for j in range(m + 1):
                 op_res = operator(x_new, p)
@@ -131,14 +133,17 @@ def si_leqi(operator, timesteps, power=None, power_density=None,
 
                 inputs = list(zip(op_res_last.rates, op_results[0].rates,
                                   op_res_bar.rates, repeat(dt_l), repeat(dt)))
-                x_new = deplete(chain, x[0], inputs, dt, print_out,
-                                matrix_func=_leqi_f3)
-                x_new = deplete(chain, x_new, inputs, dt, print_out,
-                                matrix_func=_leqi_f4)
+                time_1, x_new = timed_deplete(
+                    chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3)
+                time_2, x_new = timed_deplete(
+                    chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4)
+
+                proc_time += time_1 + time_2
 
             # Create results, write to disk
             op_results.append(op_res_bar)
-            Results.save(operator, x, op_results, [t, t+dt], p, i_res+i)
+            Results.save(
+                operator, x, op_results, [t, t+dt], p, i_res+i, proc_time)
 
             # update results
             x = [x_new]
@@ -148,4 +153,5 @@ def si_leqi(operator, timesteps, power=None, power_density=None,
             dt_l = dt
 
         # Create results for last point, write to disk
-        Results.save(operator, x, op_results, [t, t], p, i_res+len(timesteps))
+        Results.save(
+                operator, x, op_results, [t, t], p, i_res+len(timesteps))

openmc/deplete/results.py

@@ -10,7 +10,7 @@ from warnings import warn
 import numpy as np
 import h5py
 
-from . import comm, have_mpi
+from . import comm, have_mpi, MPI
 from .reaction_rates import ReactionRates
 
 _VERSION_RESULTS = (1, 0)
@@ -47,6 +47,9 @@ class Results(object):
         Number of stages in simulation.
     data : numpy.ndarray
         Atom quantity, stored by stage, mat, then by nuclide.
+    proc_time: int
+        Average time spent depleting a material across all
+        materials and processes
 
     """
     def __init__(self):
@@ -55,6 +58,7 @@ class Results(object):
         self.power = None
         self.rates = None
         self.volume = None
+        self.proc_time = None
 
         self.mat_to_ind = None
         self.nuc_to_ind = None
@@ -245,6 +249,10 @@ class Results(object):
         handle.create_dataset("power", (1, n_stages), maxshape=(None, n_stages),
                               dtype='float64')
 
+        handle.create_dataset(
+            "depletion time", (1,), maxshape=(None,),
+            dtype="float64")
+
     def _to_hdf5(self, handle, index):
         """Converts results object into an hdf5 object.
 
@@ -268,6 +276,7 @@ class Results(object):
         eigenvalues_dset = handle["/eigenvalues"]
         time_dset = handle["/time"]
         power_dset = handle["/power"]
+        proc_time_dset = handle["/depletion time"]
 
         # Get number of results stored
         number_shape = list(number_dset.shape)
@@ -296,6 +305,10 @@ class Results(object):
             power_shape[0] = new_shape
             power_dset.resize(power_shape)
 
+            proc_shape = list(proc_time_dset.shape)
+            proc_shape[0] = new_shape
+            proc_time_dset.resize(proc_shape)
+
         # If nothing to write, just return
         if len(self.mat_to_ind) == 0:
             return
@@ -314,6 +327,10 @@ class Results(object):
         if comm.rank == 0:
             time_dset[index, :] = self.time
             power_dset[index, :] = self.power
+            if self.proc_time is not None:
+                proc_time_dset[index] = (
+                    self.proc_time / (comm.size * self.n_hdf5_mats)
+                )
 
     @classmethod
     def from_hdf5(cls, handle, step):
@@ -324,8 +341,7 @@ class Results(object):
         handle : h5py.File or h5py.Group
             An HDF5 file or group type to load from.
         step : int
-            What step is this?
-
+            Index for depletion step
         """
         results = cls()
 
@@ -340,6 +356,14 @@ class Results(object):
         results.time = time_dset[step, :]
         results.power = power_dset[step, :]
 
+        if "depletion time" in handle:
+            proc_time_dset = handle["/depletion time"]
+            if step < proc_time_dset.shape[0]:
+                results.proc_time = proc_time_dset[step]
+
+        if results.proc_time is None:
+            results.proc_time = np.array([np.nan])
+
         # Reconstruct dictionaries
         results.volume = OrderedDict()
         results.mat_to_ind = OrderedDict()
@@ -375,7 +399,7 @@ class Results(object):
         return results
 
     @staticmethod
-    def save(op, x, op_results, t, power, step_ind):
+    def save(op, x, op_results, t, power, step_ind, proc_time=None):
         """Creates and writes depletion results to disk
 
         Parameters
@@ -392,6 +416,10 @@ class Results(object):
             Power during time step
         step_ind : int
             Step index.
+        proc_time : float or None
+            Total process time spent depleting materials. This may
+            be process-dependent and will be reduced across MPI
+            processes.
 
         """
         # Get indexing terms
@@ -422,6 +450,9 @@ class Results(object):
         results.rates = [r.rates for r in op_results]
         results.time = t
         results.power = power
+        results.proc_time = proc_time
+        if results.proc_time is not None:
+            results.proc_time = comm.reduce(proc_time, op=MPI.SUM)
 
         results.export_to_hdf5("depletion_results.h5", step_ind)
 

openmc/deplete/results_list.py

@@ -105,3 +105,33 @@ class ResultsList(list):
             eigenvalue[i] = result.k[0]
 
         return time, eigenvalue
+
+    def get_depletion_time(self):
+        """Return an array of the average time to deplete a material
+
+        ..note::
+
+            Will have one fewer row than number of other methods,
+            like :meth:`get_eigenvalues`, because no depletion
+            is performed at the final transport stage
+
+        Returns
+        -------
+
+        times : :class:`numpy.ndarray`
+            Vector of average time to deplete a single material
+            across all processes and materials.
+
+        """
+        times = np.empty(len(self) - 1)
+        # Need special logic because the predictor
+        # writes EOS values for step i as BOS values
+        # for step i+1
+        # The first proc_time may be zero
+        if self[0].proc_time > 0.0:
+            items = self[:-1]
+        else:
+            items = self[1:]
+        for ix, res in enumerate(items):
+            times[ix] = res.proc_time
+        return times

tests/unit_tests/test_deplete_cecm.py

@@ -35,3 +35,9 @@ def test_cecm(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)

tests/unit_tests/test_deplete_celi.py

@@ -35,3 +35,9 @@ def test_celi(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)

tests/unit_tests/test_deplete_cf4.py

@@ -35,3 +35,9 @@ def test_cf4(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)

tests/unit_tests/test_deplete_epc_rk4.py

@@ -35,3 +35,9 @@ def test_epc_rk4(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)

tests/unit_tests/test_deplete_leqi.py

@@ -35,3 +35,9 @@ def test_leqi(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)

tests/unit_tests/test_deplete_predictor.py

@@ -35,3 +35,9 @@ def test_predictor(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)

tests/unit_tests/test_deplete_si_celi.py

@@ -35,3 +35,9 @@ def test_si_celi(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)

tests/unit_tests/test_deplete_si_leqi.py

@@ -35,3 +35,9 @@ def test_si_leqi(run_in_tmpdir):
 
     assert y1[2] == approx(s2[0])
     assert y2[2] == approx(s2[1])
+
+    # Test structure of depletion time dataset
+
+    dep_time = res.get_depletion_time()
+    assert dep_time.shape == (len(dt), )
+    assert all(dep_time > 0)