crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Introduce ResultsList class to replace scattered functions

Browse source
This commit is contained in:
Paul Romano 2018-02-21 15:36:27 -06:00
parent 7ef5174274
commit 92697ec06e
15 changed files with 220 additions and 262 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/pythonapi/deplete.rst
View file

@ -59,6 +59,7 @@ data, such as number densities and reaction rates for each material.
OperatorResult
ReactionRates
Results
ResultsList
TransportOperator
Each of the integrator functions also relies on a number of "helper" functions
@ -71,4 +72,3 @@ as follows:
integrator.CRAM16
integrator.CRAM48
integrator.save_results

2
openmc/deplete/__init__.py
View file

@ -20,5 +20,5 @@ from .operator import *
from .reaction_rates import *
from .abc import *
from .results import *
from .results_list import *
from .integrator import *
from .utilities import *

1
openmc/deplete/integrator/__init__.py
View file

@ -8,4 +8,3 @@ The integrator subcomponents.
from .cecm import *
from .cram import *
from .predictor import *
from .save_results import *

16
openmc/deplete/integrator/cecm.py
View file

@ -4,7 +4,7 @@ import copy
from collections.abc import Iterable
from .cram import deplete
from .save_results import save_results
from ..results import Results
def cecm(operator, timesteps, power, print_out=True):
@ -51,20 +51,20 @@ def cecm(operator, timesteps, power, print_out=True):
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep reaction rates
x = [copy.deepcopy(vec)]
results = [operator(x[0], p)]
op_results = [operator(x[0], p)]
# Deplete for first half of timestep
x_middle = deplete(chain, x[0], results[0], dt/2, print_out)
x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out)
# Get middle-of-timestep reaction rates
x.append(x_middle)
results.append(operator(x_middle, p))
op_results.append(operator(x_middle, p))
# Deplete for full timestep using beginning-of-step materials
x_end = deplete(chain, x[0], results[1], dt, print_out)
x_end = deplete(chain, x[0], op_results[1], dt, print_out)
# Create results, write to disk
save_results(operator, x, results, [t, t + dt], i)
Results.save(operator, x, op_results, [t, t + dt], i)
# Advance time, update vector
t += dt
@ -72,7 +72,7 @@ def cecm(operator, timesteps, power, print_out=True):
# Perform one last simulation
x = [copy.deepcopy(vec)]
results = [operator(x[0], power[-1])]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
save_results(operator, x, results, [t, t], len(timesteps))
Results.save(operator, x, op_results, [t, t], len(timesteps))

12
openmc/deplete/integrator/predictor.py
View file

@ -4,7 +4,7 @@ import copy
from collections.abc import Iterable
from .cram import deplete
from .save_results import save_results
from ..results import Results
def predictor(operator, timesteps, power, print_out=True):
@ -46,13 +46,13 @@ def predictor(operator, timesteps, power, print_out=True):
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep reaction rates
x = [copy.deepcopy(vec)]
results = [operator(x[0], p)]
op_results = [operator(x[0], p)]
# Create results, write to disk
save_results(operator, x, results, [t, t + dt], i)
Results.save(operator, x, op_results, [t, t + dt], i)
# Deplete for full timestep
x_end = deplete(chain, x[0], results[0], dt, print_out)
x_end = deplete(chain, x[0], op_results[0], dt, print_out)
# Advance time, update vector
t += dt
@ -60,7 +60,7 @@ def predictor(operator, timesteps, power, print_out=True):
# Perform one last simulation
x = [copy.deepcopy(vec)]
results = [operator(x[0], power[-1])]
op_results = [operator(x[0], power[-1])]
# Create results, write to disk
save_results(operator, x, results, [t, t], len(timesteps))
Results.save(operator, x, op_results, [t, t], len(timesteps))

42
openmc/deplete/integrator/save_results.py
View file

@ -1,42 +0,0 @@
""" Generic result saving code for integrators.
"""
from ..results import Results, write_results
def save_results(op, x, op_results, t, step_ind):
"""Creates and writes depletion results to disk
Parameters
----------
op : openmc.deplete.TransportOperator
The operator used to generate these results.
x : list of list of numpy.array
The prior x vectors. Indexed [i][cell] using the above equation.
op_results : list of openmc.deplete.OperatorResult
Results of applying transport operator
t : list of float
Time indices.
step_ind : int
Step index.
"""
# Get indexing terms
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
# Create results
stages = len(x)
results = Results()
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
n_mat = len(burn_list)
for i in range(stages):
for mat_i in range(n_mat):
results[i, mat_i, :] = x[i][mat_i][:]
results.k = [r.k for r in op_results]
results.rates = [r.rates for r in op_results]
results.time = t
write_results(results, "depletion_results.h5", step_ind)

93
openmc/deplete/results.py
View file

@ -11,7 +11,6 @@ import h5py
from . import comm, have_mpi
from .reaction_rates import ReactionRates
from openmc.checkvalue import check_filetype_version
_VERSION_RESULTS = (1, 0)
@ -146,6 +145,27 @@ class Results(object):
# Create storage array
self.data = np.zeros((stages, self.n_mat, self.n_nuc))
def export_to_hdf5(self, filename, step):
"""Export results to an HDF5 file
Parameters
----------
filename : str
The filename to write to
step : int
What step is this?
"""
if have_mpi and h5py.get_config().mpi:
kwargs = {'driver': 'mpio', 'comm': comm}
else:
kwargs = {}
kwargs['mode'] = "w" if step == 0 else "a"
with h5py.File(filename, **kwargs) as handle:
self._to_hdf5(handle, step)
def _write_hdf5_metadata(self, handle):
"""Writes result metadata in HDF5 file
@ -217,7 +237,7 @@ class Results(object):
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
def to_hdf5(self, handle, index):
def _to_hdf5(self, handle, index):
"""Converts results object into an hdf5 object.
Parameters
@ -338,49 +358,40 @@ class Results(object):
return results
@staticmethod
def save(op, x, op_results, t, step_ind):
"""Creates and writes depletion results to disk
def write_results(result, filename, step):
"""Outputs result to an HDF5 file.
Parameters
----------
op : openmc.deplete.TransportOperator
The operator used to generate these results.
x : list of list of numpy.array
The prior x vectors. Indexed [i][cell] using the above equation.
op_results : list of openmc.deplete.OperatorResult
Results of applying transport operator
t : list of float
Time indices.
step_ind : int
Step index.
Parameters
----------
result : Results
Object to be stored in a file.
filename : String
Target filename.
step : int
What step is this?
"""
# Get indexing terms
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
"""
if have_mpi and h5py.get_config().mpi:
kwargs = {'driver': 'mpio', 'comm': comm}
else:
kwargs = {}
# Create results
stages = len(x)
results = Results()
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
kwargs['mode'] = "w" if step == 0 else "a"
n_mat = len(burn_list)
with h5py.File(filename, **kwargs) as handle:
result.to_hdf5(handle, step)
for i in range(stages):
for mat_i in range(n_mat):
results[i, mat_i, :] = x[i][mat_i][:]
results.k = [r.k for r in op_results]
results.rates = [r.rates for r in op_results]
results.time = t
def read_results(filename):
"""Return a list of Results objects from an HDF5 file.
Parameters
----------
filename : str
The filename to read from.
Returns
-------
results : list of openmc.deplete.Results
The result objects.
"""
with h5py.File(str(filename), "r") as fh:
check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
# Get number of results stored
n = fh["number"].value.shape[0]
return [Results.from_hdf5(fh, i) for i in range(n)]
results.export_to_hdf5("depletion_results.h5", step_ind)

105
openmc/deplete/results_list.py Normal file
View file

@ -0,0 +1,105 @@
import h5py
import numpy as np
from .results import Results, _VERSION_RESULTS
from openmc.checkvalue import check_filetype_version
class ResultsList(list):
"""A list of openmc.deplete.Results objects
Parameters
----------
filename : str
The filename to read from.
"""
def __init__(self, filename):
super().__init__()
with h5py.File(str(filename), "r") as fh:
check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
# Get number of results stored
n = fh["number"].value.shape[0]
for i in range(n):
self.append(Results.from_hdf5(fh, i))
def get_atoms(self, mat, nuc):
"""Get nuclide concentration over time from a single material
Parameters
----------
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
Returns
-------
time : numpy.ndarray
Array of times in [s]
concentration : numpy.ndarray
Total number of atoms for specified nuclide
"""
time = np.empty_like(self)
concentration = np.empty_like(self)
# Evaluate value in each region
for i, result in enumerate(self):
time[i] = result.time[0]
concentration[i] = result[0, mat, nuc]
return time, concentration
def get_reaction_rate(self, mat, nuc, rx):
"""Get reaction rate in a single material/nuclide over time
Parameters
----------
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
rx : str
Reaction rate to evaluate
Returns
-------
time : numpy.ndarray
Array of times in [s]
rate : numpy.ndarray
Array of reaction rates
"""
time = np.empty_like(self)
rate = np.empty_like(self)
# Evaluate value in each region
for i, result in enumerate(self):
time[i] = result.time[0]
rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
return time, rate
def get_eigenvalue(self):
"""Evaluates the eigenvalue from a results list.
Returns
-------
time : numpy.ndarray
Array of times in [s]
eigenvalue : numpy.ndarray
k-eigenvalue at each time
"""
time = np.empty_like(self)
eigenvalue = np.empty_like(self)
# Get time/eigenvalue at each point
for i, result in enumerate(self):
time[i] = result.time[0]
eigenvalue[i] = result.k[0]
return time, eigenvalue

101
openmc/deplete/utilities.py
View file

@ -1,101 +0,0 @@
"""The utilities module.
Contains functions that can be used to post-process objects that come out of
the results module.
"""
import numpy as np
def evaluate_single_nuclide(results, mat, nuc):
"""Evaluates a single nuclide in a single material from a results list.
Parameters
----------
results : list of results
The results to extract data from. Must be sorted and continuous.
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
Returns
-------
time : numpy.ndarray
Time vector
concentration : numpy.ndarray
Total number of atoms in the material
"""
n_points = len(results)
time = np.zeros(n_points)
concentration = np.zeros(n_points)
# Evaluate value in each region
for i, result in enumerate(results):
time[i] = result.time[0]
concentration[i] = result[0, mat, nuc]
return time, concentration
def evaluate_reaction_rate(results, mat, nuc, rx):
"""Return reaction rate in a single material/nuclide from a results list.
Parameters
----------
results : list of openmc.deplete.Results
The results to extract data from. Must be sorted and continuous.
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
rx : str
Reaction rate to evaluate
Returns
-------
time : numpy.ndarray
Time vector.
rate : numpy.ndarray
Reaction rate.
"""
n_points = len(results)
time = np.zeros(n_points)
rate = np.zeros(n_points)
# Evaluate value in each region
for i, result in enumerate(results):
time[i] = result.time[0]
rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
return time, rate
def evaluate_eigenvalue(results):
"""Evaluates the eigenvalue from a results list.
Parameters
----------
results : list of openmc.deplete.Results
The results to extract data from. Must be sorted and continuous.
Returns
-------
time : numpy.ndarray
Time vector.
eigenvalue : numpy.ndarray
Eigenvalue.
"""
n_points = len(results)
time = np.zeros(n_points)
eigenvalue = np.zeros(n_points)
# Evaluate value in each region
for i, result in enumerate(results):
time[i] = result.time[0]
eigenvalue[i] = result.k[0]
return time, eigenvalue

42
tests/dummy_operator.py
View file

@ -34,19 +34,15 @@ class DummyOperator(TransportOperator):
Returns
-------
k : float
Zero.
rates : ReactionRates
Reaction rates from this simulation.
seed : int
Zero.
openmc.deplete.OperatorResult
Result of transport operator
"""
mats = ["1"]
nuclides = ["1", "2"]
reactions = ["1"]
cell_to_ind = {"1" : 0}
nuc_to_ind = {"1" : 0, "2" : 1}
react_to_ind = {"1" : 0}
reaction_rates = ReactionRates(cell_to_ind, nuc_to_ind, react_to_ind)
reaction_rates = ReactionRates(mats, nuclides, reactions)
reaction_rates[0, 0, 0] = vec[0][0]
reaction_rates[0, 1, 0] = vec[0][1]
@ -105,18 +101,18 @@ class DummyOperator(TransportOperator):
return ["1", "2"]
@property
def burn_list(self):
def local_mats(self):
"""
burn_list : list of str
A list of all cell IDs to be burned. Used for sorting the simulation.
local_mats : list of str
A list of all material IDs to be burned. Used for sorting the simulation.
"""
return ["1"]
@property
def mat_tally_ind(self):
def burnable_mats(self):
"""Maps cell name to index in global geometry."""
return {"1": 0}
return self.local_mats
@property
@ -125,11 +121,11 @@ class DummyOperator(TransportOperator):
reaction_rates : ReactionRates
Reaction rates from the last operator step.
"""
cell_to_ind = {"1" : 0}
nuc_to_ind = {"1" : 0, "2" : 1}
react_to_ind = {"1" : 0}
mats = ["1"]
nuclides = ["1", "2"]
reactions = ["1"]
return ReactionRates(cell_to_ind, nuc_to_ind, react_to_ind)
return ReactionRates(mats, nuclides, reactions)
def initial_condition(self):
"""Returns initial vector.
@ -153,8 +149,8 @@ class DummyOperator(TransportOperator):
A list of all nuclide names. Used for sorting the simulation.
burn_list : list of int
A list of all cell IDs to be burned. Used for sorting the simulation.
full_burn_dict : OrderedDict of str to int
full_burn_list : OrderedDict of str to int
Maps cell name to index in global geometry.
"""
return self.volume, self.nuc_list, self.burn_list, self.mat_tally_ind
"""
return self.volume, self.nuc_list, self.local_mats, self.burnable_mats

10
tests/regression_tests/test_deplete_full.py
View file

@ -8,8 +8,6 @@ import numpy as np
import openmc
from openmc.data import JOULE_PER_EV
import openmc.deplete
from openmc.deplete import results
from openmc.deplete import utilities
from tests.regression_tests import config
from .example_geometry import generate_problem
@ -67,8 +65,8 @@ def test_full(run_in_tmpdir):
return
# Load the reference/test results
res_test = results.read_results(path_test)
res_ref = results.read_results(path_reference)
res_test = openmc.deplete.ResultsList(path_test)
res_ref = openmc.deplete.ResultsList(path_reference)
# Assert same mats
for mat in res_ref[0].mat_to_ind:
@ -88,8 +86,8 @@ def test_full(run_in_tmpdir):
tol = 1.0e-6
for mat in res_test[0].mat_to_ind:
for nuc in res_test[0].nuc_to_ind:
_, y_test = utilities.evaluate_single_nuclide(res_test, mat, nuc)
_, y_old = utilities.evaluate_single_nuclide(res_ref, mat, nuc)
_, y_test = res_test.get_atoms(mat, nuc)
_, y_old = res_ref.get_atoms(mat, nuc)
# Test each point
correct = True

8
tests/unit_tests/test_deplete_cecm.py
View file

@ -5,8 +5,6 @@ These tests integrate a simple test problem described in dummy_geometry.py.
from pytest import approx
import openmc.deplete
from openmc.deplete import results
from openmc.deplete import utilities
from tests import dummy_operator
@ -23,10 +21,10 @@ def test_cecm(run_in_tmpdir):
openmc.deplete.cecm(op, dt, power, print_out=False)
# Load the files
res = results.read_results(op.output_dir / "depletion_results.h5")
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
_, y1 = utilities.evaluate_single_nuclide(res, "1", "1")
_, y2 = utilities.evaluate_single_nuclide(res, "1", "2")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")
# Mathematica solution
s1 = [1.86872629872102, 1.395525772416039]

12
tests/unit_tests/test_deplete_integrator.py
View file

@ -1,4 +1,4 @@
"""Tests for integrator.py
"""Tests for saving results
It is worth noting that openmc.deplete.integrate is extremely complex, to the
point I am unsure if it can be reasonably unit-tested. For the time being, it
@ -11,11 +11,11 @@ import os
from unittest.mock import MagicMock
import numpy as np
from openmc.deplete import (integrator, ReactionRates, results, comm,
from openmc.deplete import (ReactionRates, Results, ResultsList, comm,
OperatorResult)
def test_save_results(run_in_tmpdir):
def test_results_save(run_in_tmpdir):
"""Test data save module"""
stages = 3
@ -72,11 +72,11 @@ def test_save_results(run_in_tmpdir):
op_result1 = [OperatorResult(k, rates) for k, rates in zip(eigvl1, rate1)]
op_result2 = [OperatorResult(k, rates) for k, rates in zip(eigvl2, rate2)]
integrator.save_results(op, x1, op_result1, t1, 0)
integrator.save_results(op, x2, op_result2, t2, 1)
Results.save(op, x1, op_result1, t1, 0)
Results.save(op, x2, op_result2, t2, 1)
# Load the files
res = results.read_results("depletion_results.h5")
res = ResultsList("depletion_results.h5")
for i in range(stages):
for mat_i, mat in enumerate(burn_list):

8
tests/unit_tests/test_deplete_predictor.py
View file

@ -5,8 +5,6 @@ These tests integrate a simple test problem described in dummy_geometry.py.
from pytest import approx
import openmc.deplete
from openmc.deplete import results
from openmc.deplete import utilities
from tests import dummy_operator
@ -23,10 +21,10 @@ def test_predictor(run_in_tmpdir):
openmc.deplete.predictor(op, dt, power, print_out=False)
# Load the files
res = results.read_results(op.output_dir / "depletion_results.h5")
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
_, y1 = utilities.evaluate_single_nuclide(res, "1", "1")
_, y2 = utilities.evaluate_single_nuclide(res, "1", "2")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")
# Mathematica solution
s1 = [2.46847546272295, 0.986431226850467]

28
tests/regression_tests/test_deplete_utilities.py → tests/unit_tests/test_deplete_resultslist.py
View file

@ -1,26 +1,22 @@
""" Tests the utilities classes.
This also tests the results read/write code.
"""
"""Tests the ResultsList class"""
from pathlib import Path
import numpy as np
import pytest
from openmc.deplete import results
from openmc.deplete import utilities
import openmc.deplete
@pytest.fixture
def res():
"""Load the reference results"""
filename = Path(__file__).with_name('test_reference.h5')
return results.read_results(filename)
filename = Path(__file__).parents[1] / 'regression_tests' / 'test_reference.h5'
return openmc.deplete.ResultsList(filename)
def test_evaluate_single_nuclide(res):
"""Tests evaluating single nuclide utility code."""
t, n = utilities.evaluate_single_nuclide(res, "1", "Xe135")
def test_get_atoms(res):
"""Tests evaluating single nuclide concentration."""
t, n = res.get_atoms("1", "Xe135")
t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
n_ref = [6.6747328233649218e+08, 3.5421791038348462e+14,
@ -29,9 +25,9 @@ def test_evaluate_single_nuclide(res):
np.testing.assert_array_equal(t, t_ref)
np.testing.assert_array_equal(n, n_ref)
def test_evaluate_reaction_rate(res):
"""Tests evaluating reaction rate utility code."""
t, r = utilities.evaluate_reaction_rate(res, "1", "Xe135", "(n,gamma)")
def test_get_reaction_rate(res):
"""Tests evaluating reaction rate."""
t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)")
t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
n_ref = np.array([6.6747328233649218e+08, 3.5421791038348462e+14,
@ -43,9 +39,9 @@ def test_evaluate_reaction_rate(res):
np.testing.assert_array_equal(r, n_ref * xs_ref)
def test_evaluate_eigenvalue(res):
def test_get_eigenvalue(res):
"""Tests evaluating eigenvalue."""
t, k = utilities.evaluate_eigenvalue(res)
t, k = res.get_eigenvalue()
t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
k_ref = [1.181281798790367, 1.1798750921988739, 1.1965943696058159,