Merge pull request #1357 from drewejohnson/internal-tally

Create and use internal tallies for depletion
This commit is contained in:
Paul Romano 2019-10-01 08:26:12 -05:00 committed by GitHub
commit db8cdff1c5
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
12 changed files with 164 additions and 63 deletions

View file

@ -109,6 +109,10 @@ The current version of the statepoint file format is 17.0.
**/tallies/tally <uid>/**
:Attributes: - **internal** (*int*) -- Flag indicating the presence of tally
data (0) or absence of tally data (1). All user defined
tallies will have a value of 0 unless otherwise instructed.
:Datasets: - **n_realizations** (*int*) -- Number of realizations.
- **n_filters** (*int*) -- Number of filters used.
- **filters** (*int[]*) -- User-defined unique IDs of the filters on

View file

@ -110,6 +110,7 @@ extern "C" {
int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n);
int openmc_tally_get_scores(int32_t index, int** scores, int* n);
int openmc_tally_get_type(int32_t index, int32_t* type);
int openmc_tally_get_writable(int32_t index, bool* writable);
int openmc_tally_reset(int32_t index);
int openmc_tally_results(int32_t index, double** ptr, size_t shape_[3]);
int openmc_tally_set_active(int32_t index, bool active);
@ -119,6 +120,7 @@ extern "C" {
int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides);
int openmc_tally_set_scores(int32_t index, int n, const char** scores);
int openmc_tally_set_type(int32_t index, const char* type);
int openmc_tally_set_writable(int32_t index, bool writable);
int openmc_zernike_filter_get_order(int32_t index, int* order);
int openmc_zernike_filter_get_params(int32_t index, double* x, double* y, double* r);
int openmc_zernike_filter_set_order(int32_t index, int order);

View file

@ -37,6 +37,8 @@ public:
void set_active(bool active) { active_ = active; }
void set_writable(bool writable) { writable_ = writable; }
void set_scores(pugi::xml_node node);
void set_scores(const std::vector<std::string>& scores);
@ -55,6 +57,8 @@ public:
int32_t n_filter_bins() const {return n_filter_bins_;}
bool writable() const { return writable_;}
//----------------------------------------------------------------------------
// Other methods.
@ -98,6 +102,9 @@ public:
//! (e.g. specific cell, specific energy group, etc.)
xt::xtensor<double, 3> results_;
//! True if this tally should be written to statepoint files
bool writable_ {true};
//----------------------------------------------------------------------------
// Miscellaneous public members.

View file

@ -546,6 +546,7 @@ class TalliedFissionYieldHelper(FissionYieldHelper):
# Tally group-wise fission reaction rates
self._fission_rate_tally = Tally()
self._fission_rate_tally.writable = False
self._fission_rate_tally.scores = ['fission']
self._fission_rate_tally.filters = [MaterialFilter(materials)]

View file

@ -59,6 +59,7 @@ class DirectReactionRateHelper(ReactionRateHelper):
``"(n, gamma)"``, needed for the reaction rate tally.
"""
self._rate_tally = Tally()
self._rate_tally.writable = False
self._rate_tally.scores = scores
self._rate_tally.filters = [MaterialFilter(materials)]
@ -194,6 +195,7 @@ class EnergyScoreHelper(EnergyHelper):
"""
self._tally = Tally()
self._tally.writable = False
self._tally.scores = [self.score]
def reset(self):
@ -570,6 +572,7 @@ class AveragedFissionYieldHelper(TalliedFissionYieldHelper):
func_filter = EnergyFunctionFilter()
func_filter.set_data((0, self._upper_energy), (0, self._upper_energy))
weighted_tally = Tally()
weighted_tally.writable = False
weighted_tally.scores = ['fission']
weighted_tally.filters = filters + [func_filter]
self._weighted_tally = weighted_tally

View file

@ -53,6 +53,9 @@ _dll.openmc_tally_get_scores.errcheck = _error_handler
_dll.openmc_tally_get_type.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_type.restype = c_int
_dll.openmc_tally_get_type.errcheck = _error_handler
_dll.openmc_tally_get_writable.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_tally_get_writable.restype = c_int
_dll.openmc_tally_get_writable.errcheck = _error_handler
_dll.openmc_tally_reset.argtypes = [c_int32]
_dll.openmc_tally_reset.restype = c_int
_dll.openmc_tally_reset.errcheck = _error_handler
@ -81,6 +84,9 @@ _dll.openmc_tally_set_scores.errcheck = _error_handler
_dll.openmc_tally_set_type.argtypes = [c_int32, c_char_p]
_dll.openmc_tally_set_type.restype = c_int
_dll.openmc_tally_set_type.errcheck = _error_handler
_dll.openmc_tally_set_writable.argtypes = [c_int32, c_bool]
_dll.openmc_tally_set_writable.restype = c_int
_dll.openmc_tally_set_writable.errcheck = _error_handler
_dll.tallies_size.restype = c_size_t
@ -344,6 +350,16 @@ class Tally(_FortranObjectWithID):
return std_dev
@property
def writable(self):
writable = c_bool()
_dll.openmc_tally_get_writable(self._index, writable)
return writable.value
@writable.setter
def writable(self, writable):
_dll.openmc_tally_set_writable(self._index, writable)
def reset(self):
"""Reset results and num_realizations of tally"""
_dll.openmc_tally_reset(self._index)

View file

@ -375,13 +375,18 @@ class StatePoint(object):
for tally_id in tally_ids:
group = tallies_group['tally {}'.format(tally_id)]
# Read the number of realizations
n_realizations = group['n_realizations'][()]
# Check if tally is internal and therefore has no data
if group.attrs.get("internal"):
continue
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_id)
tally._sp_filename = self._f.filename
tally.name = group['name'][()].decode() if 'name' in group else ''
# Read the number of realizations
n_realizations = group['n_realizations'][()]
tally.estimator = group['estimator'][()].decode()
tally.num_realizations = n_realizations

View file

@ -638,6 +638,16 @@ write_tallies()
for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
const auto& tally {*model::tallies[i_tally]};
// Write header block.
std::string tally_header("TALLY " + std::to_string(tally.id_));
if (!tally.name_.empty()) tally_header += ": " + tally.name_;
tallies_out << header(tally_header) << "\n\n";
if (!tally.writable_) {
tallies_out << " Internal\n\n";
continue;
}
// Calculate t-value for confidence intervals
double t_value = 1;
if (settings::confidence_intervals) {
@ -645,11 +655,6 @@ write_tallies()
t_value = t_percentile(1 - alpha*0.5, tally.n_realizations_ - 1);
}
// Write header block.
std::string tally_header("TALLY " + std::to_string(tally.id_));
if (!tally.name_.empty()) tally_header += ": " + tally.name_;
tallies_out << header(tally_header) << "\n\n";
// Write derivative information.
if (tally.deriv_ != C_NONE) {
const auto& deriv {model::tally_derivs[tally.deriv_]};

View file

@ -165,36 +165,43 @@ openmc_statepoint_write(const char* filename, bool* write_source)
write_attribute(tallies_group, "ids", tally_ids);
// Write all tally information except results
for (const auto& tally_ptr : model::tallies) {
const auto& tally {*tally_ptr};
for (const auto& tally : model::tallies) {
hid_t tally_group = create_group(tallies_group,
"tally " + std::to_string(tally.id_));
"tally " + std::to_string(tally->id_));
write_dataset(tally_group, "name", tally.name_);
write_dataset(tally_group, "name", tally->name_);
if (tally.estimator_ == ESTIMATOR_ANALOG) {
if (tally->writable_) {
write_attribute(tally_group, "internal", 0);
} else {
write_attribute(tally_group, "internal", 1);
close_group(tally_group);
continue;
}
if (tally->estimator_ == ESTIMATOR_ANALOG) {
write_dataset(tally_group, "estimator", "analog");
} else if (tally.estimator_ == ESTIMATOR_TRACKLENGTH) {
} else if (tally->estimator_ == ESTIMATOR_TRACKLENGTH) {
write_dataset(tally_group, "estimator", "tracklength");
} else if (tally.estimator_ == ESTIMATOR_COLLISION) {
} else if (tally->estimator_ == ESTIMATOR_COLLISION) {
write_dataset(tally_group, "estimator", "collision");
}
write_dataset(tally_group, "n_realizations", tally.n_realizations_);
write_dataset(tally_group, "n_realizations", tally->n_realizations_);
// Write the ID of each filter attached to this tally
write_dataset(tally_group, "n_filters", tally.filters().size());
if (!tally.filters().empty()) {
write_dataset(tally_group, "n_filters", tally->filters().size());
if (!tally->filters().empty()) {
std::vector<int32_t> filter_ids;
filter_ids.reserve(tally.filters().size());
for (auto i_filt : tally.filters())
filter_ids.reserve(tally->filters().size());
for (auto i_filt : tally->filters())
filter_ids.push_back(model::tally_filters[i_filt]->id());
write_dataset(tally_group, "filters", filter_ids);
}
// Write the nuclides this tally scores
std::vector<std::string> nuclides;
for (auto i_nuclide : tally.nuclides_) {
for (auto i_nuclide : tally->nuclides_) {
if (i_nuclide == -1) {
nuclides.push_back("total");
} else {
@ -207,12 +214,12 @@ openmc_statepoint_write(const char* filename, bool* write_source)
}
write_dataset(tally_group, "nuclides", nuclides);
if (tally.deriv_ != C_NONE) write_dataset(tally_group, "derivative",
model::tally_derivs[tally.deriv_].id);
if (tally->deriv_ != C_NONE) write_dataset(tally_group, "derivative",
model::tally_derivs[tally->deriv_].id);
// Write the tally score bins
std::vector<std::string> scores;
for (auto sc : tally.scores_) scores.push_back(reaction_name(sc));
for (auto sc : tally->scores_) scores.push_back(reaction_name(sc));
write_dataset(tally_group, "n_score_bins", scores.size());
write_dataset(tally_group, "score_bins", scores);
@ -232,6 +239,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
// Write all tally results
for (const auto& tally : model::tallies) {
if (!tally->writable_) continue;
// Write sum and sum_sq for each bin
std::string name = "tally " + std::to_string(tally->id_);
hid_t tally_group = open_group(tallies_group, name.c_str());
@ -425,11 +433,21 @@ void load_state_point()
// Read sum, sum_sq, and N for each bin
std::string name = "tally " + std::to_string(tally->id_);
hid_t tally_group = open_group(tallies_group, name.c_str());
auto& results = tally->results_;
read_tally_results(tally_group, results.shape()[0],
results.shape()[1], results.data());
read_dataset(tally_group, "n_realizations", tally->n_realizations_);
close_group(tally_group);
int internal=0;
if (attribute_exists(tally_group, "internal")) {
read_attribute(tally_group, "internal", internal);
}
if (internal) {
tally->writable_ = false;
} else {
auto& results = tally->results_;
read_tally_results(tally_group, results.shape()[0],
results.shape()[1], results.data());
read_dataset(tally_group, "n_realizations", tally->n_realizations_);
close_group(tally_group);
}
}
close_group(tallies_group);
@ -697,6 +715,7 @@ void write_tally_results_nr(hid_t file_id)
for (const auto& t : model::tallies) {
// Skip any tallies that are not active
if (!t->active_) continue;
if (!t->writable_) continue;
if (mpi::master && !object_exists(file_id, "tallies_present")) {
write_attribute(file_id, "tallies_present", 1);

View file

@ -1220,6 +1220,30 @@ openmc_tally_set_active(int32_t index, bool active)
return 0;
}
extern "C" int
openmc_tally_get_writable(int32_t index, bool* writable)
{
if (index < 0 || index >= model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*writable = model::tallies[index]->writable();
return 0;
}
extern "C" int
openmc_tally_set_writable(int32_t index, bool writable)
{
if (index < 0 || index >= model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
model::tallies[index]->set_writable(writable);
return 0;
}
extern "C" int
openmc_tally_get_scores(int32_t index, int** scores, int* n)
{

View file

@ -108,11 +108,17 @@ def test_full(run_in_tmpdir):
t_ref, k_ref = res_ref.get_eigenvalue()
k_state = np.empty_like(k_ref)
n_tallies = np.empty(N + 1, dtype=int)
# Get statepoint files for all BOS points and EOL
for n in range(N + 1):
statepoint = openmc.StatePoint("openmc_simulation_n{}.h5".format(n))
k_n = statepoint.k_combined
k_state[n] = [k_n.nominal_value, k_n.std_dev]
n_tallies[n] = len(statepoint.tallies)
# Look for exact match pulling from statepoint and depletion_results
assert np.all(k_state == k_test)
assert np.allclose(k_test, k_ref)
# Check that no additional tallies are loaded from the files
assert np.all(n_tallies == 0)

View file

@ -50,24 +50,24 @@ def pincell_model():
@pytest.fixture(scope='module')
def capi_init(pincell_model, mpi_intracomm):
def lib_init(pincell_model, mpi_intracomm):
openmc.lib.init(intracomm=mpi_intracomm)
yield
openmc.lib.finalize()
@pytest.fixture(scope='module')
def capi_simulation_init(capi_init):
def lib_simulation_init(lib_init):
openmc.lib.simulation_init()
yield
@pytest.fixture(scope='module')
def capi_run(capi_simulation_init):
def lib_run(lib_simulation_init):
openmc.lib.run()
def test_cell_mapping(capi_init):
def test_cell_mapping(lib_init):
cells = openmc.lib.cells
assert isinstance(cells, Mapping)
assert len(cells) == 3
@ -76,7 +76,7 @@ def test_cell_mapping(capi_init):
assert cell_id == cell.id
def test_cell(capi_init):
def test_cell(lib_init):
cell = openmc.lib.cells[1]
assert isinstance(cell.fill, openmc.lib.Material)
cell.fill = openmc.lib.materials[1]
@ -85,7 +85,7 @@ def test_cell(capi_init):
cell.name = "Not fuel"
assert cell.name == "Not fuel"
def test_cell_temperature(capi_init):
def test_cell_temperature(lib_init):
cell = openmc.lib.cells[1]
cell.set_temperature(100.0, 0)
assert cell.get_temperature(0) == 100.0
@ -93,7 +93,7 @@ def test_cell_temperature(capi_init):
assert cell.get_temperature() == 200.0
def test_new_cell(capi_init):
def test_new_cell(lib_init):
with pytest.raises(exc.AllocationError):
openmc.lib.Cell(1)
new_cell = openmc.lib.Cell()
@ -101,7 +101,7 @@ def test_new_cell(capi_init):
assert len(openmc.lib.cells) == 5
def test_material_mapping(capi_init):
def test_material_mapping(lib_init):
mats = openmc.lib.materials
assert isinstance(mats, Mapping)
assert len(mats) == 3
@ -110,7 +110,7 @@ def test_material_mapping(capi_init):
assert mat_id == mat.id
def test_material(capi_init):
def test_material(lib_init):
m = openmc.lib.materials[3]
assert m.nuclides == ['H1', 'O16', 'B10', 'B11']
@ -136,14 +136,14 @@ def test_material(capi_init):
m.name = "Not hot borated water"
assert m.name == "Not hot borated water"
def test_material_add_nuclide(capi_init):
def test_material_add_nuclide(lib_init):
m = openmc.lib.materials[3]
m.add_nuclide('Xe135', 1e-12)
assert m.nuclides[-1] == 'Xe135'
assert m.densities[-1] == 1e-12
def test_new_material(capi_init):
def test_new_material(lib_init):
with pytest.raises(exc.AllocationError):
openmc.lib.Material(1)
new_mat = openmc.lib.Material()
@ -151,7 +151,7 @@ def test_new_material(capi_init):
assert len(openmc.lib.materials) == 5
def test_nuclide_mapping(capi_init):
def test_nuclide_mapping(lib_init):
nucs = openmc.lib.nuclides
assert isinstance(nucs, Mapping)
assert len(nucs) == 13
@ -160,7 +160,7 @@ def test_nuclide_mapping(capi_init):
assert name == nuc.name
def test_settings(capi_init):
def test_settings(lib_init):
settings = openmc.lib.settings
assert settings.batches == 10
settings.batches = 10
@ -175,7 +175,7 @@ def test_settings(capi_init):
settings.run_mode = 'eigenvalue'
def test_tally_mapping(capi_init):
def test_tally_mapping(lib_init):
tallies = openmc.lib.tallies
assert isinstance(tallies, Mapping)
assert len(tallies) == 3
@ -184,7 +184,7 @@ def test_tally_mapping(capi_init):
assert tally_id == tally.id
def test_energy_function_filter(capi_init):
def test_energy_function_filter(lib_init):
"""Test special __new__ and __init__ for EnergyFunctionFilter"""
efunc = openmc.lib.EnergyFunctionFilter([0.0, 1.0], [0.0, 2.0])
assert len(efunc.energy) == 2
@ -193,7 +193,7 @@ def test_energy_function_filter(capi_init):
assert (efunc.y == [0.0, 2.0]).all()
def test_tally(capi_init):
def test_tally(lib_init):
t = openmc.lib.tallies[1]
assert t.type == 'volume'
assert len(t.filters) == 2
@ -239,7 +239,7 @@ def test_tally(capi_init):
assert len(t3_f.y) == 3
def test_new_tally(capi_init):
def test_new_tally(lib_init):
with pytest.raises(exc.AllocationError):
openmc.lib.Material(1)
new_tally = openmc.lib.Tally()
@ -249,16 +249,25 @@ def test_new_tally(capi_init):
assert len(openmc.lib.tallies) == 5
def test_tally_activate(capi_simulation_init):
def test_tally_activate(lib_simulation_init):
t = openmc.lib.tallies[1]
assert not t.active
t.active = True
assert t.active
def test_tally_results(capi_run):
def test_tally_writable(lib_simulation_init):
t = openmc.lib.tallies[1]
assert t.num_realizations == 10 # t was made active in test_tally
assert t.writable
t.writable = False
assert not t.writable
# Revert tally to writable state for lib_run fixtures
t.writable = True
def test_tally_results(lib_run):
t = openmc.lib.tallies[1]
assert t.num_realizations == 10 # t was made active in test_tally_active
assert np.all(t.mean >= 0)
nonzero = (t.mean > 0.0)
assert np.all(t.std_dev[nonzero] >= 0)
@ -269,26 +278,26 @@ def test_tally_results(capi_run):
assert t2.mean.size == (n + 1) * (n + 2) // 2 * 3 # Number of Zernike coeffs * 3 cells
def test_global_tallies(capi_run):
def test_global_tallies(lib_run):
assert openmc.lib.num_realizations() == 5
gt = openmc.lib.global_tallies()
for mean, std_dev in gt:
assert mean >= 0
def test_statepoint(capi_run):
def test_statepoint(lib_run):
openmc.lib.statepoint_write('test_sp.h5')
assert os.path.exists('test_sp.h5')
def test_source_bank(capi_run):
def test_source_bank(lib_run):
source = openmc.lib.source_bank()
assert np.all(source['E'] > 0.0)
assert np.all(source['wgt'] == 1.0)
assert np.allclose(np.linalg.norm(source['u'], axis=1), 1.0)
def test_by_batch(capi_run):
def test_by_batch(lib_run):
openmc.lib.hard_reset()
# Running next batch before simulation is initialized should raise an
@ -313,7 +322,7 @@ def test_by_batch(capi_run):
openmc.lib.simulation_finalize()
def test_reset(capi_run):
def test_reset(lib_run):
# Init and run 10 batches.
openmc.lib.hard_reset()
openmc.lib.simulation_init()
@ -344,7 +353,7 @@ def test_reset(capi_run):
openmc.lib.simulation_finalize()
def test_reproduce_keff(capi_init):
def test_reproduce_keff(lib_init):
# Get k-effective after run
openmc.lib.hard_reset()
openmc.lib.run()
@ -357,7 +366,7 @@ def test_reproduce_keff(capi_init):
assert keff0 == pytest.approx(keff1)
def test_find_cell(capi_init):
def test_find_cell(lib_init):
cell, instance = openmc.lib.find_cell((0., 0., 0.))
assert cell is openmc.lib.cells[1]
cell, instance = openmc.lib.find_cell((0.4, 0., 0.))
@ -366,14 +375,14 @@ def test_find_cell(capi_init):
openmc.lib.find_cell((100., 100., 100.))
def test_find_material(capi_init):
def test_find_material(lib_init):
mat = openmc.lib.find_material((0., 0., 0.))
assert mat is openmc.lib.materials[1]
mat = openmc.lib.find_material((0.4, 0., 0.))
assert mat is openmc.lib.materials[2]
def test_mesh(capi_init):
def test_mesh(lib_init):
mesh = openmc.lib.RegularMesh()
mesh.dimension = (2, 3, 4)
assert mesh.dimension == (2, 3, 4)
@ -408,7 +417,7 @@ def test_mesh(capi_init):
assert msf.mesh == mesh
def test_restart(capi_init, mpi_intracomm):
def test_restart(lib_init, mpi_intracomm):
# Finalize and re-init to make internal state consistent with XML.
openmc.lib.hard_reset()
openmc.lib.finalize()
@ -440,7 +449,7 @@ def test_restart(capi_init, mpi_intracomm):
assert keff0 == pytest.approx(keff1)
def test_load_nuclide(capi_init):
def test_load_nuclide(lib_init):
# load multiple nuclides
openmc.lib.load_nuclide('H3')
assert 'H3' in openmc.lib.nuclides
@ -451,7 +460,7 @@ def test_load_nuclide(capi_init):
openmc.lib.load_nuclide('Pu3')
def test_id_map(capi_init):
def test_id_map(lib_init):
expected_ids = np.array([[(3, 3), (2, 2), (3, 3)],
[(2, 2), (1, 1), (2, 2)],
[(3, 3), (2, 2), (3, 3)]], dtype='int32')
@ -469,7 +478,7 @@ def test_id_map(capi_init):
ids = openmc.lib.plot.id_map(s)
assert np.array_equal(expected_ids, ids)
def test_property_map(capi_init):
def test_property_map(lib_init):
expected_properties = np.array(
[[(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)],
[ (293.6, 6.55), (293.6, 10.29769), (293.6, 6.55)],
@ -489,7 +498,7 @@ def test_property_map(capi_init):
assert np.allclose(expected_properties, properties, atol=1e-04)
def test_position(capi_init):
def test_position(lib_init):
pos = openmc.lib.plot._Position(1.0, 2.0, 3.0)
@ -502,7 +511,7 @@ def test_position(capi_init):
assert tuple(pos) == (1.3, 2.3, 3.3)
def test_global_bounding_box(capi_init):
def test_global_bounding_box(lib_init):
expected_llc = (-0.63, -0.63, -np.inf)
expected_urc = (0.63, 0.63, np.inf)