diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 9bfe4aa9e..ebc5126df 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -42,9 +42,11 @@ extern "C" {
int openmc_global_tallies(double** ptr);
int openmc_hard_reset();
int openmc_init(int argc, char* argv[], const void* intracomm);
+ bool openmc_is_statepoint_batch();
int openmc_legendre_filter_get_order(int32_t index, int* order);
int openmc_legendre_filter_set_order(int32_t index, int order);
int openmc_load_nuclide(const char* name);
+ bool openmc_master();
int openmc_material_add_nuclide(int32_t index, const char name[], double density);
int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n);
int openmc_material_get_id(int32_t index, int32_t* id);
diff --git a/openmc/capi/core.py b/openmc/capi/core.py
index 922877c8e..b85adb32f 100644
--- a/openmc/capi/core.py
+++ b/openmc/capi/core.py
@@ -47,6 +47,7 @@ _init_linsolver_argtypes = [_array_1d_int, c_int, _array_1d_int, c_int, c_int,
c_double, _array_1d_int, _array_1d_int]
_dll.openmc_initialize_linsolver.argtypes = _init_linsolver_argtypes
_dll.openmc_initialize_linsolver.restype = None
+_dll.openmc_is_statepoint_batch.restype = c_bool
_dll.openmc_master.restype = c_bool
_dll.openmc_next_batch.argtypes = [POINTER(c_int)]
_dll.openmc_next_batch.restype = c_int
@@ -186,6 +187,18 @@ def init(args=None, intracomm=None):
_dll.openmc_init(argc, argv, intracomm)
+def is_statepoint_batch():
+ """Return whether statepoint will be written in current batch or not.
+
+ Returns
+ -------
+ bool
+ Whether is statepoint batch or not
+
+ """
+ return _dll.openmc_is_statepoint_batch()
+
+
def iter_batches():
"""Iterator over batches.
diff --git a/openmc/cmfd.py b/openmc/cmfd.py
index 0a5ba7cd1..60156c7f5 100644
--- a/openmc/cmfd.py
+++ b/openmc/cmfd.py
@@ -17,9 +17,11 @@ import sys
import time
from ctypes import c_int, c_double
import warnings
+import os.path
import numpy as np
from scipy import sparse
+import h5py
import openmc.capi
from openmc.checkvalue import (check_type, check_length, check_value,
@@ -767,6 +769,10 @@ class CMFDRun(object):
# Write CMFD output if CMFD on for current batch
if openmc.capi.master():
self._write_cmfd_output()
+
+ # Write CMFD data to statepoint
+ if openmc.capi.is_statepoint_batch():
+ self.statepoint_write()
return status
def finalize(self):
@@ -777,6 +783,64 @@ class CMFDRun(object):
# Print out CMFD timing statistics
self._write_cmfd_timing_stats()
+ def statepoint_write(self, filename=None):
+ if filename is None:
+ batch_str_len = len(str(openmc.capi.settings.batches))
+ batch_str = str(openmc.capi.current_batch()).zfill(batch_str_len)
+ filename = 'statepoint.{}.h5'.format(batch_str)
+
+ # Call C API statepoint_write if statepoint file doesn't exist
+ if not os.path.isfile(filename):
+ openmc.capi.statepoint_write(filename=filename)
+
+ # Append CMFD data to statepoint file using h5py
+ self._write_cmfd_statepoint(filename)
+
+ def _write_cmfd_statepoint(self, filename):
+ if openmc.capi.master():
+ if openmc.capi.settings.verbosity >= 5:
+ print(' Writing CMFD data to {}'.format(filename))
+ sys.stdout.flush()
+ # TODO write CMFD data to statepoint file
+ # TODO check if "cmfd" group already exists
+ with h5py.File(filename, 'a') as h5f:
+ cmfd_group = h5f.create_group("cmfd")
+ cmfd_group.attrs['cmfd_on'] = self._cmfd_on
+ cmfd_group.attrs['feedback_begin'] = self._feedback_begin
+ cmfd_group.attrs['mesh_id'] = self._mesh_id
+ cmfd_group.attrs['tally_begin'] = self._tally_begin
+ cmfd_group.attrs['time_cmfd'] = self._time_cmfd
+ cmfd_group.attrs['time_cmfdbuild'] = self._time_cmfdbuild
+ cmfd_group.attrs['time_cmfdsolve'] = self._time_cmfdsolve
+ cmfd_group.attrs['window_size'] = self._window_size
+ cmfd_group.attrs['window_type'] = self._window_type
+ cmfd_group.create_dataset('k_cmfd', data=self._k_cmfd)
+ cmfd_group.create_dataset('dom', data=self._dom)
+ cmfd_group.create_dataset('src_cmp', data=self._src_cmp)
+ cmfd_group.create_dataset('balance', data=self._balance)
+ cmfd_group.create_dataset('entropy', data=self._entropy)
+ cmfd_group.create_dataset('albedo', data=self._albedo)
+ cmfd_group.create_dataset('coremap', data=self._coremap)
+ cmfd_group.create_dataset('egrid', data=self._egrid)
+ cmfd_group.create_dataset('indices', data=self._indices)
+ cmfd_group.create_dataset('tally_ids', data=self._tally_ids)
+
+ if self._cmfd_on:
+ cmfd_group.create_dataset('current_rate',
+ data=self._current_rate)
+ cmfd_group.create_dataset('flux_rate',
+ data=self._flux_rate)
+ cmfd_group.create_dataset('nfiss_rate',
+ data=self._nfiss_rate)
+ cmfd_group.create_dataset('openmc_src_rate',
+ data=self._openmc_src_rate)
+ cmfd_group.create_dataset('p1scatt_rate',
+ data=self._p1scatt_rate)
+ cmfd_group.create_dataset('scatt_rate',
+ data=self._scatt_rate)
+ cmfd_group.create_dataset('total_rate',
+ data=self._total_rate)
+
def _initialize_linsolver(self):
# Determine number of rows in CMFD matrix
ng = self._indices[3]
@@ -841,6 +905,21 @@ class CMFDRun(object):
# Set up CMFD coremap
self._set_coremap()
+ # Extract spatial and energy indices
+ nx = self._indices[0]
+ ny = self._indices[1]
+ nz = self._indices[2]
+ ng = self._indices[3]
+
+ # Allocate parameters that need to stored for rolling window
+ self._openmc_src_rate = np.zeros((nx, ny, nz, ng, 0))
+ self._flux_rate = np.zeros((nx, ny, nz, ng, 0))
+ self._total_rate = np.zeros((nx, ny, nz, ng, 0))
+ self._p1scatt_rate = np.zeros((nx, ny, nz, ng, 0))
+ self._scatt_rate = np.zeros((nx, ny, nz, ng, ng, 0))
+ self._nfiss_rate = np.zeros((nx, ny, nz, ng, ng, 0))
+ self._current_rate = np.zeros((nx, ny, nz, 12, ng, 0))
+
# Initialize timers
self._time_cmfd = 0.0
self._time_cmfdbuild = 0.0
@@ -904,8 +983,7 @@ class CMFDRun(object):
'feedback begin')
# Set number of batches where cmfd tallies should be reset
- if self._reset is not None:
- self._n_resets = len(self._reset)
+ self._n_resets = len(self._reset)
# Create tally objects
self._create_cmfd_tally()
@@ -922,15 +1000,6 @@ class CMFDRun(object):
self._hxyz = np.zeros((nx, ny, nz, 3))
self._hxyz[:,:,:,:] = openmc.capi.meshes[self._mesh_id].width
- # Allocate parameters that need to stored for rolling window
- self._openmc_src_rate = np.zeros((nx, ny, nz, ng, 0))
- self._flux_rate = np.zeros((nx, ny, nz, ng, 0))
- self._total_rate = np.zeros((nx, ny, nz, ng, 0))
- self._p1scatt_rate = np.zeros((nx, ny, nz, ng, 0))
- self._scatt_rate = np.zeros((nx, ny, nz, ng, ng, 0))
- self._nfiss_rate = np.zeros((nx, ny, nz, ng, ng, 0))
- self._current_rate = np.zeros((nx, ny, nz, 12, ng, 0))
-
# Allocate flux, cross sections and diffusion coefficient
self._flux = np.zeros((nx, ny, nz, ng))
self._totalxs = np.zeros((nx, ny, nz, ng))
diff --git a/src/message_passing.cpp b/src/message_passing.cpp
index 278bf924d..e53c65825 100644
--- a/src/message_passing.cpp
+++ b/src/message_passing.cpp
@@ -1,4 +1,5 @@
#include "openmc/message_passing.h"
+#include "openmc/capi.h"
namespace openmc {
namespace mpi {
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 7d30ccea9..c0e534fc3 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -233,6 +233,16 @@ int openmc_next_batch(int* status)
return 0;
}
+bool openmc_is_statepoint_batch() {
+ using namespace openmc;
+ using openmc::simulation::current_gen;
+
+ if (!simulation::initialized)
+ return false;
+ else
+ return contains(settings::statepoint_batch, simulation::current_batch);
+}
+
namespace openmc {
//==============================================================================
diff --git a/tests/regression_tests/cmfd_nofeed/settings.xml b/tests/regression_tests/cmfd_nofeed/settings.xml
index 24b0b6ab5..5c89af8d2 100644
--- a/tests/regression_tests/cmfd_nofeed/settings.xml
+++ b/tests/regression_tests/cmfd_nofeed/settings.xml
@@ -23,4 +23,6 @@
10
+
+