diff --git a/.gitignore b/.gitignore
index c61074d..001bbd2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,4 +2,7 @@ build
misc
commit_msg
__pycache__
+
+*.egg-info
*.swp
+*.json
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d0455b7..a5892cb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -120,11 +120,20 @@ FIND_PACKAGE(COIN REQUIRED)
SET(STUB_INCLUDE_DIRS ${STUB_INCLUDE_DIRS} ${COIN_INCLUDE_DIR})
set(LIBS ${LIBS} ${COIN_LIBRARIES})
+# Find Python and python.h
+FIND_PACKAGE(PythonLibs REQUIRED)
+SET(STUB_INCLUDE_DIRS ${STUB_INCLUDE_DIRS} ${PYTHON_INCLUDE_DIRS})
+
+# Find Niels Lohmann's 'JSON for Modern C++' package
+FIND_PACKAGE(nlohmann_json REQUIRED)
+# The line below does not work!
+# For unknown reasons, CMake does not find the include directories
+# and nlohmann_json_INCLUDE_DIRS is an unset variable.
+SET(STUB_INCLUDE_DIRS ${STUB_INCLUDE_DIRS} ${nlohmann_json_INCLUDE_DIRS})
# include all the directories we just found
INCLUDE_DIRECTORIES(${STUB_INCLUDE_DIRS})
-
# add the agents
ADD_SUBDIRECTORY(src)
diff --git a/LICENSE b/LICENSE
index 4a914fa..dfa5f20 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
BSD 3-Clause License
-Copyright (c) 2020, RWTH Aachen University Nuclear Verification and Disarmament Group
+Copyright (c) 2020-2021, RWTH Aachen University Nuclear Verification and Disarmament Group
All rights reserved.
Redistribution and use in source and binary forms, with or without
diff --git a/README.md b/README.md
index b8b381a..2e5acba 100644
--- a/README.md
+++ b/README.md
@@ -1,19 +1,34 @@
# misoenrichment: a Cyclus multi component isotope enrichment module
`misoenrichment` is a module developed at the [Nuclear Verification and Disarmament Group](https://www.nvd.rwth-aachen.de/) at RWTH Aachen University for the nuclear fuel cycle simulator
-[Cyclus](http://fuelcycle.org). It provides a Cyclus facility that enriches
-uranium streams composed of two or more isotopes taking into account the
+[Cyclus](http://fuelcycle.org). It currently provides two Cyclus facilities.
+
+`MIsoEnrich` is an enrichment facility that enriches
+streams composed of two or more uranium isotopes taking into account the
different enrichment behaviour of minor isotopes such as 234U (present in
natural as well as in reprocessed uranium) or 236U (present in
reprocessed uranium from spent nuclear fuel). The tracking of minor
isotopes makes this module suitable for nuclear archaeology, see, e.g., Ref 3.
+`GprReactor` is a Cyclus reactor facility that uses Gaussian Process
+Regression (GPR) to calculate the composition of the irradiated fuel depending
+on various input parameters. Generally, this implementation works for any
+reactor type and any input parameters. However, one needs the appropriate
+GPR model (which needs to be generated using training data) and depending
+on which input parameters are chosen, the source code of `GprReactor` may
+need minor tweaking. Additional information on this issue will be given
+in future commits.
+
Table of Contents
-- [Getting Started](#getting-started)
-- [Theoretical background](#theoretical-background)
+- [MIsoEnrich](#misoenrich)
+ - [Getting Started](#getting-started)
+ - [Theoretical background](#theoretical-background)
+- [GprReactor](#gprreactor)
+ - [Requirements](#requirements)
- [References](#references)
-## Getting started
+## MIsoEnrich
+### Getting started
An example input file is found in `input/main.py` featuring a
`cycamore::Source` source agent, a `MIsoEnrich` enrichment facility and two
`cycamore::Sink` agents, one for enriched and one for depleted uranium.
@@ -34,7 +49,7 @@ level, as explained above) and then blending the product with uranium from
the feed. This procedure is only performed if the `use_downblending`
variable is set to `True` in the input file.
-## Theoretical background
+### Theoretical background
The implementation of the facility itself and the interaction with Cyclus'
Dynamic Resource Exchange is based on the binary enrichment facility from
the [Cycamore](https://github.com/cyclus/cycamore) package.
@@ -44,6 +59,31 @@ Ref 1 derives the mathematical formalism of a matched abundance ratio cascade
using constant overall stage separation factors while Ref 2 gives a new
physically founded approach to calculating said separation factors.
+## GprReactor
+### Requirements
+This facility needs Niels Lohmann's [JSON for Modern C++](https://json.nlohmann.me/)
+library. It can be downloaded from his [GitHub repository](https://github.com/nlohmann/json)
+or using one of the many package managers, see [here](https://github.com/nlohmann/json#package-managers).
+Successfully tested using `conda install -c conda-forge nlohmann_json` and using `CMake`.
+When using `CMake`, then the package needs to be installed globally
+(i.e., under `/usr/local`), as shown in the following:
+```
+$ git clone https://github.com/nlohmann/json
+$ cd json
+$ mkdir build
+$cd build
+$ cmake ..
+$ make
+$ sudo make install
+```
+While it _should_ be possible to install `JSON for Modern C++` locally,
+i.e. in `~/.local`, this results in `CMake` not finding `nlohmann/json.hpp`
+during the `misoenrichment` installation. I will hopefully manage to
+fix this in future versions.
+
+Additionally, one needs `Python3` in combination with the `NumPy` and
+`SciPy` packages.
+
## References
1. E. von Halle, _Multicomponent Isotope Separation in Matched Abundance
diff --git a/input/archetypes.py b/input/archetypes.py
index dd97548..fde0ac1 100644
--- a/input/archetypes.py
+++ b/input/archetypes.py
@@ -5,7 +5,9 @@ def archetypes():
{"lib": "agents", "name": "NullRegion"},
{"lib": "cycamore", "name": "Sink"},
{"lib": "cycamore", "name": "Source"},
- {"lib": "misoenrichment", "name":"MIsoEnrich"}
+ {"lib": "cycamore", "name": "Storage"},
+ {"lib": "misoenrichment", "name": "GprReactor"},
+ {"lib": "misoenrichment", "name": "MIsoEnrich"}
]
}}
diff --git a/input/control.py b/input/control.py
index 5afe7ff..79d096a 100644
--- a/input/control.py
+++ b/input/control.py
@@ -2,8 +2,8 @@ def control():
d = {"control": {
"startyear": 2020,
"startmonth": 1,
- "duration": 12,
- "dt": 2629846, # duration of a time step in seconds, here: 1 month
+ "duration": 365,
+ "dt": 86400, # duration of a time step in seconds, here: 1 day
"simhandle": "misoenrichment tutorial file",
}
}
diff --git a/input/facility.py b/input/facility.py
index 0558f03..eb31b1b 100644
--- a/input/facility.py
+++ b/input/facility.py
@@ -5,14 +5,13 @@ def facility():
"config": {"Source": {
"outcommod": "NaturalU",
"outrecipe": "NaturalURecipe",
- "throughput": 2000
+ "throughput": 10000
}}
},
{
- "name": "EnrichedUSink",
+ "name": "SpentFuelSink",
"config": {"Sink": {
- "in_commods": {"val": ["EnrichedU"]},
- "recipe_name": "EnrichedURecipe"
+ "in_commods": {"val": ["SpentFuel"]}
}}
},
{
@@ -22,6 +21,15 @@ def facility():
"recipe_name": "DepletedURecipe"
}}
},
+ {
+ "name": "FreshFuelStorage",
+ "config": {"Storage": {
+ "in_commods": {"val": ["EnrichedU"]},
+ "out_commods": {"val": ["FreshFuel"]},
+ "in_recipe": ["EnrichedURecipe"],
+ "residence_time": 0
+ }}
+ },
{
"name": "EnrichmentFacility",
"config": {"MIsoEnrich": {
@@ -31,13 +39,28 @@ def facility():
"tails_commod": "DepletedU",
"tails_assay": 0.003,
"initial_feed": 0,
- "max_feed_inventory": 10000,
+ "max_feed_inventory": 1e299,
"gamma_235": 1.35,
"swu_capacity": 1e299,
- "swu_capacity_vals": {"val": [1e5, 1000, 2000]},
+ "swu_capacity_vals": {"val": [1e5, 5e4, 5e5]},
"swu_capacity_times": {"val": [0, 5, 6]},
"use_downblending": True
}}
+ },
+ {
+ "name": "SavannahRiverReactor",
+ "config": {"GprReactor": {
+ "in_commods": {"val": ["FreshFuel"]},
+ "out_commods": {"val": ["SpentFuel"]},
+ "in_recipes": {"val": ["EnrichedURecipe"]},
+ "n_assem_core": 1,
+ "n_assem_batch": 1,
+ "assem_size": 110820,
+ "cycle_time": 88,
+ "refuel_time": 6,
+ "power_output": 2400,
+ "temperature": 350
+ }}
}
]}
return d
diff --git a/input/recipe.py b/input/recipe.py
index 36a8005..ecf1645 100644
--- a/input/recipe.py
+++ b/input/recipe.py
@@ -4,17 +4,19 @@ def recipe():
"name": "NaturalURecipe",
"basis": "mass",
"nuclide": [
- {"id": "U234", "comp": 5.5e-3},
+ # At the moment (May 2021), U234 is not taken into
+ # account by the Gpr model.
+ #{"id": "U234", "comp": 5.5e-3},
{"id": "U235", "comp": 0.711},
- {"id": "U238", "comp": 100 - 0.711 - 5.5e-3}
+ {"id": "U238", "comp": 100 - 0.711}
]
},
{
"name": "EnrichedURecipe",
"basis": "mass",
"nuclide": [
- {"id": "U235", "comp": 90},
- {"id": "U238", "comp": 10}
+ {"id": "U235", "comp": 1.1},
+ {"id": "U238", "comp": 98.9}
]
},
{
diff --git a/install.py b/install.py
index dba1439..83df77f 100644
--- a/install.py
+++ b/install.py
@@ -37,6 +37,11 @@ def install(args):
root_dir = os.path.split(__file__)[0]
makefile = os.path.join(args.build_dir, 'Makefile')
+ print("Installing Python (sub)module 'spentfuelgpr'...")
+ spent_fuel_gpr_dir = os.path.join(root_dir, "spentfuelgpr")
+ subprocess.check_call(["pip3", "install", "-e", spent_fuel_gpr_dir])
+ print()
+
if not os.path.exists(makefile):
rtn = subprocess.call(['which', 'cmake'], shell=(os.name == 'nt'))
if rtn != 0:
diff --git a/spentfuelgpr/data/x_trainingset.npy b/spentfuelgpr/data/x_trainingset.npy
new file mode 100644
index 0000000..7f1db3c
Binary files /dev/null and b/spentfuelgpr/data/x_trainingset.npy differ
diff --git a/spentfuelgpr/data/y_trainingset_reduced.npy b/spentfuelgpr/data/y_trainingset_reduced.npy
new file mode 100644
index 0000000..0a8a6cb
Binary files /dev/null and b/spentfuelgpr/data/y_trainingset_reduced.npy differ
diff --git a/spentfuelgpr/setup.py b/spentfuelgpr/setup.py
new file mode 100644
index 0000000..409c277
--- /dev/null
+++ b/spentfuelgpr/setup.py
@@ -0,0 +1,22 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from setuptools import setup
+
+def main():
+ setup(
+ name="spentfuelgpr",
+ version="0.1",
+ description="Python part of the MIsoEnrichment module",
+ author="Nuclear Verification and Disarmament Group, RWTH Aachen University",
+ url="https://github.com/maxschalz/miso_enrichment/",
+ license="BSD-3-Clause",
+ packages=["spentfuelgpr"],
+ classifiers=["License :: OSI Approved :: BSD-3-Clause License",
+ "Programming Language :: Python :: 3"],
+ install_requires=["numpy", "scipy"]
+ )
+ return
+
+if __name__=="__main__":
+ main()
diff --git a/spentfuelgpr/spentfuelgpr/__init__.py b/spentfuelgpr/spentfuelgpr/__init__.py
new file mode 100644
index 0000000..d19b56c
--- /dev/null
+++ b/spentfuelgpr/spentfuelgpr/__init__.py
@@ -0,0 +1,24 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Calculate the composition of spent fuel based on a set of trained GPRs.
+
+This module is part of the misoenrichment module for Cyclus, see
+https://github.com/maxschalz/miso_enrichment/ .
+
+The module is strongly based on Antonio Figueroa's work, who kindly
+provided the original code. For more information, please see
+https://doi.org/10.1016/j.anucene.2020.108085
+or https://arxiv.org/abs/2006.12921
+or https://github.com/FigueroaAC/GPs-for-SpentFuel
+"""
+
+__author__ = "Nuclear Verification and Disarmament Group, RWTH Aachen University"
+__copyright__ = "Copyright 2020-2021, Nuclear Verification and Disarmament Group, RWTH Aachen University"
+__credits__ = ["Antonio Figueroa", "Max Schalz"]
+__license__ = "BSD-3-Clause"
+__version__ = "2.0"
+__maintainer__ = "Max Schalz"
+
+from .kernel import *
+from .predict_posterior import *
diff --git a/spentfuelgpr/spentfuelgpr/kernel.py b/spentfuelgpr/spentfuelgpr/kernel.py
new file mode 100644
index 0000000..76e3361
--- /dev/null
+++ b/spentfuelgpr/spentfuelgpr/kernel.py
@@ -0,0 +1,327 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Define the kernel function used to calculate the predictions.
+
+This module is strongly based on Antonio Figueroa's work, who kindly
+provided the original code. For more information, please see
+https://doi.org/10.1016/j.anucene.2020.108085
+or https://arxiv.org/abs/2006.12921
+or https://github.com/FigueroaAC/GPs-for-SpentFuel
+"""
+
+__all__ = ["Kernel"]
+
+import numpy as np
+from scipy.spatial.distance import cdist
+
+def Kernel(X1, X2, Type, *params, gradient=False):
+ '''
+ All the kernels have a bit of noise added in order to prevent the covariance
+ matrix becoming singular. The amount of noise to be added is also a
+ parameter to reconstruct
+ '''
+
+ def dif(array1, array2):
+ Output = [array1[i] - array2[i] for i in range(len(array1))]
+ Output = np.array(Output).ravel()
+ return sum(Output)
+
+ if Type == 'SQE':
+ sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
+ K = (params[0][0]**2) * np.exp(-0.5 * sqdist / params[0][1]**2 ) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ gradients = [2*params[0][0]*K,np.multiply(sqdist/(params[0][1]**3),K),
+ 2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'ASQE':
+ LAMBDA = np.eye(len(X2[0]))
+ length_scales = 1/params[0][1:-1]
+ np.fill_diagonal(LAMBDA,length_scales)
+
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ sqdist = cdist(X1 , X2, metric='sqeuclidean').T
+ K = (params[0][0]**2) * np.exp(-0.5 * sqdist) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ g = [cdist(np.expand_dims(X1[:,i],-1),
+ np.expand_dims(X2[:,i],-1),
+ metric='sqeuclidean')/(params[0][i+1]) \
+ for i in range(X1.shape[1]) ]
+ gradients = [np.multiply(g[i],K) for i in range(X1.shape[1])]
+ gradients = [2*params[0][0]*K] + gradients + [2*params[0][-1]*np.eye(X1.shape[0])]
+
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'LAP':
+ dist = np.sqrt(np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T))
+ K = (params[0][0]**2) * np.exp(-0.5 * dist / params[0][1] ) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ gradients = [2*params[0][0]*K,np.multiply(dist/(params[0][1]**2),K),
+ 2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'ALAP':
+ dist = np.sqrt(cdist(X1/params[0][1:-1],X2/params[0][1:-1], metric='euclidean').T)
+ K = (params[0][0]**2) * np.exp(-0.5 * dist) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ g = [cdist(np.expand_dims(X1[:,i]/params[0][i+1],-1),
+ np.expand_dims(X2[:,i]/params[0][i+1],-1),
+ metric='euclidean')/(params[0][i+1]) \
+ for i in range(X1.shape[1]) ]
+ gradients = [np.multiply(g[i],K) for i in range(X1.shape[1])]
+ gradients = [2*params[0][0]*K] + gradients + [2*params[0][-1]*np.eye(X1.shape[0])]
+
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'Linear':
+ # This is a version of Linear modified to include lengthscales for each
+ # Input Variable
+# =============================================================================
+# K = a*X.Y+b
+# =============================================================================
+ LAMBDA = np.eye(len(params[0][2:-1]))
+ length_scales = 1/params[0][2:-1]
+ np.fill_diagonal(LAMBDA,length_scales)
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ K = (params[0][0]*np.dot(X1,X2.T).T)+params[0][1] +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+
+ if gradient:
+ gradients = [np.dot(X1,X2.T).T] + [np.eye(K.shape[0])] + [2*params[0][0]*\
+ np.dot(np.expand_dims(X1[:,i],-1)
+ ,np.expand_dims(X2[:,i],-1).T).T/\
+ params[0][i]\
+ for i in range(X1.shape[1])]
+ gradients = gradients + [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'Poly':
+# =============================================================================
+# K = (a*(X.Y)+b)**c
+# =============================================================================
+ LAMBDA = np.eye(len(params[0][3:-1]))
+ length_scales = 1/params[0][3:-1]
+ np.fill_diagonal(LAMBDA,length_scales)
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ K = ((params[0][0]*np.dot(X1,X2.T).T)+params[0][1])**params[0][2] +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+
+ if gradient:
+ db = params[0][2]*((params[0][0]*np.dot(X1,X2.T))+params[0][1])**(params[0][2]-1)
+ da = np.multiply(np.dot(X1,X2.T),db)
+ dc = np.multiply(K,np.log((params[0][0]*np.dot(X1,X2.T)) + params[0][1]))
+ gradients = [da,db,dc]+\
+ [2*params[0][0]*np.multiply(np.dot(np.expand_dims(X1[:,i],-1),
+ np.expand_dims(X2[:,i],-1).T)/\
+ params[0][i],db) \
+ for i in range(X1.shape[1])]
+ gradients = gradients + [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+ return
+ if Type == 'Anova':
+ # This isn't working, probably some parameters for lengthscales are needed for each dimension
+ K = np.exp((-params[0][0]*cdist(X1,X2, metric='sqeuclidean'))**params[0][1]) +\
+ np.exp((-params[0][0]*cdist(X1**2,X2**2, metric='sqeuclidean'))**params[0][1])+ \
+ np.exp((-params[0][0]*cdist(X1**3,X2**3, metric='sqeuclidean'))**params[0][1]) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))# Not working
+ if gradient:
+ dd = params[0][1] * (np.exp((-params[0][0]*cdist(X1 , X2, metric='sqeuclidean'))**(params[0][1]-1)) +\
+ np.exp((-params[0][0] *cdist(X1**2 , X2**2, metric='sqeuclidean'))**(params[0][1]-1)) + \
+ np.exp((-params[0][0]*cdist(X1**3 , X2**3, metric='sqeuclidean'))**(params[0][1]-1)))
+ gradients = [-np.multiply(cdist(X1 , X2, metric='sqeuclidean'),dd),
+ -np.multiply(cdist(X1**2 , X2**2, metric='sqeuclidean'),dd),
+ -np.multiply(cdist(X1**3 , X2**3, metric='sqeuclidean'),dd), dd,
+ 2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+
+ else:
+ return K
+
+ if Type == 'Sigmoid':
+ # This is not a positive semidefinite Kernel. Some tweaking has to be done
+ # to make this work. Probably along the lines of KdotK.T
+ K = np.tanh(params[0][0]*np.dot(X1,X2.T) + params[0][1]) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ sech2 = 1/(np.cosh(params[0][0]*np.dot(X1,X2.T) + params[0][1])**2)
+ gradients = [np.multiply(np.dot(X1,X2.T),sech2),sech2,
+ 2*params[0][-1]*np.eye(X1.shape[0]) ]
+ return K, gradients
+ # This isn't working, probably some parameters for lengthscales are needed for each dimension
+ else:
+ return K
+
+ if Type == 'RQ':
+ LAMBDA = np.eye(len(X1[0]))
+ length_scales = 1/params[0][1:-1]
+ np.fill_diagonal(LAMBDA,length_scales)
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ sqdist = cdist(X1 , X2, metric='sqeuclidean').T
+ K = 1 - (sqdist/(sqdist+(params[0][0])**2)) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+
+ if gradient:
+ denom = 1/((sqdist+(params[0][0])**2))**2
+ g = [cdist(np.expand_dims(X1[:,i],-1),
+ np.expand_dims(X2[:,i],-1),
+ metric='sqeuclidean')/(params[0][i+1]) \
+ for i in range(X1.shape[1]) ]
+ gradients = [-2*denom*g[i]*(params[0][0])**2 for i in range(X1.shape[1])]
+ gradients = [2*params[0][0]*sqdist*denom] + gradients +\
+ [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+
+ else:
+ return K
+
+
+ if Type == 'SRQ':
+ LAMBDA = np.eye(len(X1[0]))
+ length_scales = 1/params[0][:-1]
+ np.fill_diagonal(LAMBDA,length_scales)
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ sqdist = cdist(X1 , X2, metric='sqeuclidean').T
+ return (1/(1 + (sqdist/params[0][-1])) )**params[0][-1]
+
+ if Type == 'MultiQuad':
+ # This Kernel is not positive semidefinite so a lot of tweaking would
+ # have to be done to make this work. Maybe take some product of KdotK.T
+ LAMBDA = np.eye(len(X1[0]))
+ length_scales = 1/params[0][1:-1]
+ np.fill_diagonal(LAMBDA,length_scales)
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ sqdist = cdist(X1 , X2, metric='sqeuclidean').T
+ K = np.sqrt(sqdist + (params[0][0]**2)) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ g = [cdist(np.expand_dims(X1[:,i],-1),
+ np.expand_dims(X2[:,i],-1),
+ metric='sqeuclidean')/(params[0][i+1]) \
+ for i in range(X1.shape[1]) ]
+ gradients = [np.multiply(g[i],1/K) for i in range(X1.shape[1])]
+ gradients = [params[0][0]/K] + gradients +\
+ [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+
+ else:
+ return K
+
+ if Type == 'InvMultiQuad':
+ LAMBDA = np.eye(len(X1[0]))
+ length_scales = 1/params[0][1:-1]
+ np.fill_diagonal(LAMBDA,length_scales)
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ sqdist = np.sqrt(cdist(X1 , X2, metric='sqeuclidean').T + (params[0][0]**2))
+ K = 1/sqdist + ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ g = [cdist(np.expand_dims(X1[:,i],-1),
+ np.expand_dims(X2[:,i],-1),
+ metric='sqeuclidean')/(params[0][i+1]) \
+ for i in range(X1.shape[1]) ]
+ gradients = [-np.multiply(g[i],K**3) for i in range(X1.shape[1])]
+ gradients = [-params[0][0]*(K**3)] + gradients +\
+ [2*params[0][-1]*np.eye(X1.shape[0])]
+
+ return K, gradients
+ else:
+ return K
+ if Type == 'Wave':
+ dist = cdist(X1 , X2, metric='euclidean')
+ K = ((params[0][0]/dist) * np.sin(dist/params[0][0])) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ arg = dist/params[0][0]
+ gradients = (1/dist) * (np.sin(arg) - (np.cos(arg)/(params[0][0])**2)) +\
+ [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'Power':
+ K = cdist(X1/params[0][1:-1],X2/params[0][1:-1], metric='euclidean')**params[0][0] +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ dd = params[0][0] * cdist(X1/params[0][1:-1],X2/params[0][1:-1],
+ metric='euclidean')**(params[0][0]-1)
+ g = [cdist(np.expand_dims(X1[:,i]/params[0][i+1],-1),
+ np.expand_dims(X2[:,i]/params[0][i+1],-1),
+ metric='sqeuclidean')/(params[0][i+1]) \
+ for i in range(X1.shape[1]) ]
+ gradients = [np.multiply(g[i],dd) for i in range(X1.shape[1])]
+ gradients = [dd] + gradients + [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'Log':
+ # This is a non positive definite Kernel
+ dist = cdist(X1 , X2, metric='euclidean').T
+ K = -np.log((dist**params[0][0]) + 1) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ arg = dist**params[0][0]
+ gradients = arg * np.log(dist)/ (arg+1) +\
+ [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'Cauchy':
+ # This works very nice and fast
+ sqdist = cdist(X1 , X2, metric='sqeuclidean').T
+ K = 1/(1+(sqdist/(params[0][0]**2))) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ gradients = [sqdist/(((params[0][0]**2)+sqdist)**2)] +\
+ [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
+ if Type == 'Tstudent':
+ # This isn't working at the moment
+ LAMBDA = np.eye(len(X1[0]))
+ length_scales = 1/params[0][1:]
+ np.fill_diagonal(LAMBDA,length_scales)
+ X1 = np.dot(X1,LAMBDA)
+ X2 = np.dot(X2,LAMBDA)
+ sqdist = cdist(X1 , X2, metric='euclidean').T
+ K = 1/(1+(sqdist**params[0][0])) +\
+ ((params[0][-1]**2)*np.eye(X1.shape[0]))
+ if gradient:
+ da = -(sqdist**params[0][0])*np.log(sqdist)*(K**2)
+ arg = (sqdist**(params[0][0]-2))*(K**2)
+ gradients = [da] + [params[0][0]*np.multiply(
+ cdist(np.expand_dims(X1[:,i],-1),
+ np.expand_dims(X2[:,i],-1),
+ metric='sqeuclidean')/(params[0][i+1]),arg)
+ for i in range(X1.shape[1])] +\
+ [2*params[0][-1]*np.eye(X1.shape[0])]
+ return K, gradients
+ else:
+ return K
+
diff --git a/spentfuelgpr/spentfuelgpr/predict_posterior.py b/spentfuelgpr/spentfuelgpr/predict_posterior.py
new file mode 100644
index 0000000..53bb029
--- /dev/null
+++ b/spentfuelgpr/spentfuelgpr/predict_posterior.py
@@ -0,0 +1,275 @@
+#!/usr/bin/env python3]
+# -*- coding: utf-8 -*-
+"""Interface used by Cyclus to calculate the spent fuel composition."""
+
+__all__ = ["predict", "run_kernel"]
+
+import json
+import numpy as np
+import os
+
+from . import kernel
+
+
+#TODO
+# - Find a workaround for the ugly calculations that are currently
+# performed in 'get_input_params' for the 'power_output' parameter
+# and that are only relevant (?) for the SRS reactor.
+
+def main():
+ return
+
+def predict():
+ """Calculate the spent fuel composition."""
+ # Rationale for the choice of these isotopes:
+ # - isotope fraction > 1e-15
+ # - stable U isotopes, i.e., 1/2 time > days (omit U230, 231, 237)
+ # - stable Pu isotopes
+ # - decay into 'interesting' material (relevant for U->Np->Pu)
+ # - stable Pu isotopes
+
+ isotopes = ('U232', 'U233', 'U234', 'U235', 'U235m', 'U236', 'U238',
+ 'U239', 'U240', 'Pu238', 'Pu239', 'Pu240', 'Pu241',
+ 'Pu242', 'Pu243', 'Pu244', 'Np239', 'Np240', 'Np240m',
+ 'Np241')
+
+ # Check if the needed kernels and parameter information exist.
+ data_dir = os.path.join(os.path.split(__file__)[0], "..", "data")
+ kernel_dir = os.path.join(data_dir, "trained_kernels")
+ if not os.path.isdir(kernel_dir):
+ raise OSError("'trained_kernels' directory not found!")
+ for iso in isotopes:
+ fname = f"{iso}.npy"
+ if not os.path.isfile(os.path.join(kernel_dir, fname)):
+ msg = f"Trained kernel '{fname}' not found!"
+ raise FileNotFoundError(msg)
+
+ fname = f"training_params_{iso}.json"
+ if not os.path.isfile(os.path.join(kernel_dir, fname)):
+ msg = f"Training parameters '{fname}' not found!"
+ raise FileNotFoundError(msg)
+
+ # Load input parameters, training data and the kernel type.
+ training_data = np.load(os.path.join(data_dir, "x_trainingset.npy"),
+ allow_pickle=True)
+ if not os.path.isfile(os.path.join(data_dir,
+ "y_trainingset_reduced.npy")):
+ y_data = np.load(os.path.join(data_dir, "y_trainingset.npy"),
+ allow_pickle=True).item()
+ shrink_dictionary(
+ y_data, isotopes,
+ os.path.join(data_dir, "y_trainingset_reduced.npy"))
+
+ y_data = np.load(
+ os.path.join(data_dir, "y_trainingset_reduced.npy"),
+ allow_pickle=True).item()
+ par = ("enrichment", "temperature", "power_output", "burnup")
+ reactor_input_params = get_input_params(par)
+ reactor_input_params = np.expand_dims(reactor_input_params, axis=0)
+
+ # Calculate the spent fuel composition for all isotopes.
+ spent_fuel_composition = {"spent_fuel_composition": {}}
+ for iso in isotopes:
+ kernel_fname = os.path.join(kernel_dir, f"{iso}.npy")
+ params_fname = os.path.join(kernel_dir,
+ f"training_params_{iso}.json")
+ with open(params_fname, "r") as f:
+ data = json.load(f)
+ kernel_type = data["kernel_type"]
+ size = data["size"]
+ x_train = training_data[:size]
+ y_train = np.array(y_data[iso])[:size]
+ check_input_params(reactor_input_params, x_train)
+ trained_kernel = np.load(kernel_fname,
+ allow_pickle=True).item()
+ mass = run_kernel(reactor_input_params, x_train, y_train,
+ trained_kernel, kernel_type)
+
+ if mass < 0:
+ msg = (f"Calculated mass of {iso} in spent fuel is {mass} "
+ + "kg.\nHowever, the mass cannot be negative!\n"
+ + f"Parameters used: {reactor_input_params}")
+ raise RuntimeError(msg)
+
+ iso = iso[:-1]+"M" if iso[-1] == "m" else iso
+ spent_fuel_composition["spent_fuel_composition"][iso] = mass
+ store_results(spent_fuel_composition)
+
+def check_input_params(params, training_data):
+ """Check the validity of the input parameters.
+
+ The function raises a ValueError if the checks fail, else it has no
+ return value.
+ """
+ min_vals = np.min(training_data, axis=0)
+ max_vals = np.max(training_data, axis=0)
+ is_valid = (np.all(min_vals < params[0])
+ and np.all(params[0] < max_vals))
+
+ if not is_valid:
+ msg = ("[spentfuelgpr] One or more parameters exceed the bounds.\n"
+ + f"Minimum parameter values: {min_vals}\n"
+ + f"Actual parameter values: {params[0]}\n"
+ + f"Maximum parameter values: {max_vals}")
+ raise ValueError(msg)
+
+def shrink_dictionary(data, isotopes, fname):
+ """Remove unnecessary data from dictionary to reduce runtime
+
+ Parameters
+ ----------
+ data : dict
+ A dictionary containing all the data with the keys being the
+ isotopes.
+ isotopes : iterable
+ All the isotopes that remain in the final dictionary
+ fname : string
+ Filename of the final dictionary
+ """
+ d = {}
+ for iso in isotopes:
+ d[iso] = data[iso]
+
+ np.save(fname, d, allow_pickle=True)
+
+ return
+
+def get_input_params(pnames):
+ """Read in the Gpr input parameters.
+
+ Currently, the input filename is hardcoded. Maybe this will be
+ changed in the future but at the moment I have no idea how this
+ could be done.
+
+ Parameters
+ ----------
+ pnames : tuple
+ The names of the parameters to be extracted from the JSON
+ file.
+
+ Returns
+ -------
+ input_params : array
+ An array containing the extracted values in the same order as
+ specified in `pnames`.
+ """
+ fname = "gpr_reactor_input_params.json"
+ with open(fname, "r") as f:
+ data = json.load(f)
+ input_params = []
+
+ for param in pnames:
+ param = param.lower()
+ if param == "enrichment":
+ enrich = data["fresh_fuel_composition"]["922350000"]
+ input_params.append(enrich)
+ elif param == "power_output":
+ # These calculations below have to be performed for the
+ # Savannah River Site reactor. Currently, they are
+ # hardcoded but this is hopefully subject to change.
+ # TODO update this implementation
+ n_assemblies_tot = 500;
+ n_assemblies_model = 18;
+ feet_to_cm = 30.48;
+ assembly_length = 12; # in feet
+ power = (data[param] * n_assemblies_model
+ / n_assemblies_tot / assembly_length
+ / feet_to_cm)
+ power *= 1e6 # conversion MW to W
+ input_params.append(power)
+ else:
+ input_params.append(data[param])
+
+ return np.array(input_params)
+
+def store_results(composition):
+ """Store the composition in a .json file.
+
+ Parameters
+ ----------
+ composition : dict
+ A dictionary with the keys being the isotopes and the values
+ being the corresponding masses.
+ """
+ fname = "gpr_reactor_spent_fuel_composition.json"
+ with open(fname, "w") as f:
+ json.dump(composition, f, indent=2)
+
+def run_kernel(reactor_input_params, x_train, y_train, trained_kernel,
+ kernel_params, kernel_type='ASQE'):
+ """Calculate the mass of one isotope in the spent fuel.
+
+ Parameters
+ ----------
+ reactor_input_params
+ Input parameters used in the calculations and as defined above
+ in 'predict()' and in 'get_input_params'
+ x_train
+ Input parameters used during training
+ y_train
+ Output used during training
+ trained_kernel
+ Trained kernel, specific to the isotope in question
+ kernel_type
+ The type of the trained kernel
+
+ Returns
+ -------
+ mu_s
+ The (mean predicted) mass of the isotope in the spent fuel
+ after one irradiation period.
+ """
+ if (kernel_type != 'ASQE'):
+ msg = "Currently, only the 'ASQE' kernel type is supported."
+ raise ValueError(msg)
+
+ kernel_params = trained_kernel["Params"]
+ alpha = trained_kernel["alpha_"]
+ k_s = kernel.Kernel(reactor_input_params, x_train, kernel_type,
+ kernel_params, gradient=False)
+ mu_s = np.dot(k_s.T, alpha)[0]
+
+ # Revert the normalisation of the output which is used during
+ # training of the Gpr.
+ mu_s = mu_s * np.std(y_train) + np.mean(y_train)
+
+ return mu_s
+
+"""
+Part of the spent fuel compositions obtained from SERPENT simulations
+of one Savannah River Site reactors. Taken from Max Schalz' master
+thesis, see https://github.com/maxschalz/studious_potato/blob/main/data/SERPENT_outputs_NatU_percentages.npy
+Simulations kindly provided by Antonio Figueroa.
+
+BU: 0.5MWd 2MWd
+ U230: 2.701e-21 2.701e-21
+ U231: 1.422e-22 1.422e-22
+ U232: 2.850e-14 2.850e-14
+ U233: 3.099e-11 3.099e-11
+ U234: 8.513e-05 8.513e-05
+ U235: 1.030e-02 1.030e-02
+ U236: 9.760e-05 9.760e-05
+ U237: 8.535e-07 8.535e-07
+ U238: 9.885e-01 9.885e-01
+ U239: 5.391e-07 5.391e-07
+ U240: 1.744e-10 1.744e-10
+ U241: 1.100e-15 1.100e-15
+ U242: 3.389e-20 3.389e-20
+ Pu238: 9.167e-09 9.167e-09
+ Pu239: 4.052e-04 4.052e-04
+ Pu240: 9.193e-06 9.193e-06
+ Pu241: 5.679e-07 5.679e-07
+ Pu242: 5.122e-09 5.122e-09
+ Pu243: 1.835e-12 1.835e-12
+ Pu244: 4.362e-15 4.362e-15
+ Np239: 7.775e-05 7.775e-05
+ Np240: 1.160e-09 1.160e-09
+ Np241: 3.059e-15 3.059e-15
+ Np242: 4.438e-21 4.438e-21
+ U235m: 8.643e-10 8.643e-10
+Np240m: 2.402e-10 2.402e-10
+
+"""
+
+if __name__=="__main__":
+ main()
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b70d386..1035e03 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,6 +4,8 @@ USE_CYCLUS("misoenrichment" "enrichment_calculator")
USE_CYCLUS("misoenrichment" "miso_helper")
USE_CYCLUS("misoenrichment" "flexible_input")
+USE_CYCLUS("misoenrichment" "gpr_reactor")
+
INSTALL_CYCLUS_MODULE("misoenrichment" "./")
# install header files
diff --git a/src/gpr_reactor.cc b/src/gpr_reactor.cc
new file mode 100644
index 0000000..1983cea
--- /dev/null
+++ b/src/gpr_reactor.cc
@@ -0,0 +1,698 @@
+#define PY_SSIZE_T_CLEAN
+#include
+
+#include "gpr_reactor.h"
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "pyne.h"
+#include
+
+// Future changes relating to the implementation of Antonio's GPRs are marked
+// with the following comment:
+// TODO ANTONIO GPR
+
+namespace misoenrichment {
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+GprReactor::GprReactor(cyclus::Context* ctx)
+ : cyclus::Facility(ctx),
+ in_commods(std::vector()),
+ out_commods(std::vector()),
+ in_recipes(std::vector()),
+ fuel_prefs(std::vector()),
+ n_assem_core(0),
+ n_assem_batch(0),
+ assem_size(0),
+ n_assem_fresh(0),
+ n_assem_spent(0),
+ latitude(0.),
+ longitude(0.),
+ decom_transmute_all(false),
+ cycle_time(0),
+ refuel_time(0),
+ cycle_step(0),
+ discharged(false),
+ power_output(0.),
+ temperature(0.),
+ res_indexes(std::map()),
+ is_hybrid(true),
+ side_products(std::vector()),
+ side_product_quantity(std::vector()),
+ unique_out_commods(std::set()),
+ permitted_fresh_fuel_comps(std::set({922350000, 922380000})),
+ relevant_spent_fuel_comps(std::set(
+ {922320000, 922330000, 922340000, 922350000, 922350001, 922360000,
+ 922380000, 922390000, 922400000, 932390000, 932400000, 932400001,
+ 932410000, 942380000, 942390000, 942400000, 942410000, 942420000,
+ 942430000, 942440000}
+ )),
+ out_fname("gpr_reactor_input_params.json"),
+ in_fname("gpr_reactor_spent_fuel_composition.json") {
+ // TODO check, e.g., runtime performance to determine if calling PyStart here
+ // and doing the imports here (i.e., once) is actually faster or if this is
+ // all optimised anyway.
+ //cyclus::PyStart();
+ //PyRun_SimpleString("import setuptest");
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+GprReactor::~GprReactor() {
+ // TODO see comment in constructor
+ //cyclus::PyStop();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+std::set::Ptr> GprReactor::GetMatlBids(
+ cyclus::CommodMap::type& commod_requests) {
+ using cyclus::BidPortfolio;
+ using cyclus::Material;
+
+ if (unique_out_commods.empty()) {
+ for (int i = 0; i < out_commods.size(); ++i) {
+ unique_out_commods.insert(out_commods[i]);
+ }
+ }
+
+ std::set::Ptr> ports;
+ std::map all_mats;
+ std::set::iterator it;
+ // Loop over all out commodities.
+ for (it = unique_out_commods.begin(); it != unique_out_commods.end(); ++it) {
+ std::string commod = *it;
+ std::vector*>& reqs = commod_requests[commod];
+ if (reqs.size() == 0) {
+ continue;
+ } else {
+ all_mats = PeekSpent_();
+ }
+ cyclus::toolkit::MatVec mats = all_mats[commod];
+ if (mats.size() == 0) {
+ continue;
+ }
+ BidPortfolio::Ptr port(new BidPortfolio());
+ // Loop over all requests for the given out commodity.
+ for (int j = 0; j < reqs.size(); ++j) {
+ cyclus::Request* req = reqs[j];
+ double total_bid = 0;
+ // Loop over all available material of the given commidity.
+ for (int k = 0; k < mats.size(); ++k) {
+ Material::Ptr m = mats[k];
+ total_bid += m->quantity();
+ port->AddBid(req, m, this, true);
+ if (total_bid >= req->target()->quantity()) {
+ break;
+ }
+ }
+ }
+ double total_qty = 0;
+ for (int j = 0; j < mats.size(); ++j) {
+ total_qty += mats[j]->quantity();
+ }
+ cyclus::CapacityConstraint cap_constraint(total_qty);
+ port->AddConstraint(cap_constraint);
+ ports.insert(port);
+ }
+ return ports;
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+std::set::Ptr>
+GprReactor::GetMatlRequests() {
+ using cyclus::Material;
+ using cyclus::Request;
+ using cyclus::RequestPortfolio;
+
+ std::set::Ptr> ports;
+ Material::Ptr m;
+
+ int n_assem_order = n_assem_core - core.count()
+ + n_assem_fresh - fresh_inv.count();
+ // This if clause accounts for the fact that less assemblies may be needed if
+ // the reactor retirement is near.
+ if (exit_time() != -1) {
+ // The '+ 1' accounts for the fact that the reactor is online and gets to
+ // operate during its exit_time timestep.
+ int time_left = exit_time() - context()->time() + 1;
+ int time_left_cycle = cycle_time + refuel_time - cycle_step;
+ double n_cycles_left = static_cast(time_left - time_left_cycle)
+ / static_cast(cycle_time + refuel_time);
+ n_cycles_left = std::ceil(n_cycles_left);
+ int n_needed = std::max(0.0, n_cycles_left*n_assem_batch - n_assem_fresh
+ + n_assem_core - core.count());
+ n_assem_order = std::min(n_assem_order, n_needed);
+ }
+
+ if (n_assem_order == 0 || Retired_()) {
+ return ports;
+ }
+
+ // Make a request for each assembly.
+ for (int i = 0; i < n_assem_order; ++i) {
+ RequestPortfolio::Ptr port(new RequestPortfolio());
+ std::vector*> mutual_reqs;
+ // Make mutual requests for each fuel incommod.
+ for (int j = 0; j < in_commods.size(); ++j) {
+ std::string commod = in_commods[j];
+ double pref = fuel_prefs[j];
+ cyclus::Composition::Ptr recipe = context()->GetRecipe(in_recipes[j]);
+ m = Material::CreateUntracked(assem_size, recipe);
+
+ Request* r = port->AddRequest(m, this, commod, pref, true);
+ mutual_reqs.push_back(r);
+ }
+ std::vector::iterator result;
+ result = std::max_element(fuel_prefs.begin(), fuel_prefs.end());
+ int max_index = std::distance(fuel_prefs.begin(), result);
+
+ cyclus::toolkit::RecordTimeSeries("demand" + in_commods[max_index],
+ this, assem_size*n_assem_order);
+ port->AddMutualReqs(mutual_reqs);
+ ports.insert(port);
+ }
+ return ports;
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+std::string GprReactor::str() { return cyclus::Facility::str(); }
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::AcceptMatlTrades(
+ const std::vector,
+ cyclus::Material::Ptr> >& responses) {
+ std::vector,
+ cyclus::Material::Ptr> >::const_iterator trade;
+ std::stringstream ss;
+ // Number of assemblies that are loaded directly into the core.
+ int n_load = std::min((int)responses.size(), n_assem_core - core.count());
+ if (n_load > 0) {
+ ss << n_load << " assemblies";
+ Record_("LOAD", ss.str());
+ }
+
+ // Accept trades and push material to core or fresh fuel inventory.
+ for (trade = responses.begin(); trade != responses.end(); ++trade) {
+ std::string commod = trade->first.request->commodity();
+ cyclus::Material::Ptr m = trade->second;
+ IndexRes_(m, commod);
+
+ if (core.count() < n_assem_core) {
+ core.Push(m);
+ } else {
+ // TODO check if it is assured that the fresh inventory does not obtain
+ // more fuel than it can accept.
+ fresh_inv.Push(m);
+ }
+ }
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::EnterNotify() {
+ cyclus::Facility::EnterNotify();
+
+ if (fuel_prefs.size() == 0) {
+ for (int i = 0; i < out_commods.size(); i++) {
+ fuel_prefs.push_back(cyclus::kDefaultPref);
+ }
+ }
+
+ // Check if side products have been defined.
+ if (side_products.size() == 0) {
+ is_hybrid = false;
+ }
+ RecordPosition_();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::GetMatlTrades(
+ const std::vector >& trades,
+ std::vector,
+ cyclus::Material::Ptr> >& responses) {
+ std::map mats = PopSpent_();
+ for (int i = 0; i < trades.size(); i++) {
+ std::string commod = trades[i].request->commodity();
+ cyclus::Material::Ptr m = mats[commod].back();
+ mats[commod].pop_back();
+ responses.push_back(std::make_pair(trades[i], m));
+ res_indexes.erase(m->obj_id());
+ }
+ PushSpent_(mats); // return leftovers back to spent buffer
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::Tick() {
+ // Check if the reactor is already retired.
+ if (Retired_()) {
+ Record_("RETIRED", "");
+ // Transmute remaining fuel exactly once.
+ if (context()->time() == exit_time()+1) {
+ if (decom_transmute_all == true) {
+ Transmute_(std::ceil(static_cast(n_assem_core)));
+ }
+ else {
+ Transmute_(std::ceil(static_cast(n_assem_core)/2.));
+ }
+ }
+ // Empty the reactor core if this has not yet been done.
+ while (core.count() > 0) {
+ if (!Discharge_()) {
+ break;
+ }
+ }
+
+ while (fresh_inv.count() > 0 && spent_inv.space() >= assem_size) {
+ spent_inv.Push(fresh_inv.Pop());
+ }
+ if (CheckDecommissionCondition()) {
+ Decommission();
+ }
+ return;
+ }
+
+ // "Burn" the fuel, i.e., change its composition from fresh to spent fuel.
+ if (cycle_step == cycle_time) {
+ Transmute_();
+ Record_("CYCLE_END", "");
+ }
+
+ // If the irradiation period is over and the fuel has not yet been discharged,
+ // e.g., because of a full spent fuel inventory, then discharge it now if
+ // possible.
+ if (cycle_step >= cycle_time && !discharged) {
+ discharged = Discharge_();
+ }
+
+ // If the irradiation period is over, try to load fresh fuel into the reactor
+ // core.
+ if (cycle_step >= cycle_time) {
+ Load_();
+ }
+
+ // In cycamore's Reactor implementation, the changing of preferences and
+ // recipes would take place here.
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::Tock() {
+ using cyclus::toolkit::RecordTimeSeries;
+
+ if (Retired_()) {
+ return;
+ }
+
+ // Check that the irradiation and refuelling period is over, that the core is
+ // full and that the spent fuel was successfully discharged in this refuelling
+ // period. If all of this is the case, start a new cycle.
+ if (cycle_step >= cycle_time+refuel_time && core.count() == n_assem_core
+ && discharged == true) {
+ discharged = false;
+ cycle_step = 0;
+ }
+
+ if (cycle_step == 0 && core.count() == n_assem_core) {
+ Record_("CYCLE_START", "");
+ }
+
+ // Normal reactor operation where power is produced
+ if (cycle_step >= 0 && cycle_step < cycle_time
+ && core.count() == n_assem_core) {
+ RecordTimeSeries(this, power_output);
+ RecordTimeSeries("supplyPOWER", this, power_output);
+ RecordSideProduct_(true);
+ } else {
+ RecordTimeSeries(this, 0);
+ RecordTimeSeries("supplyPOWER", this, 0);
+ RecordSideProduct_(false);
+ }
+
+ // This statement prevents that a newly deployed reactor (`cycle_step = 0`)
+ // increments `cycle_step` although the core might not have been filled yet.
+ if (cycle_step > 0 || core.count() == n_assem_core) {
+ cycle_step++;
+ }
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Below are the private helper functions, not directly interacting with CYCLUS.
+bool GprReactor::CheckDecommissionCondition() {
+ return core.count()==0 && spent_inv.count()==0;
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+bool GprReactor::Discharge_() {
+ int n_pop = std::min(n_assem_batch, core.count());
+ if (n_assem_spent-spent_inv.count() < n_pop) {
+ // Not enough room in spent fuel buffer.
+ Record_("DISCHARGE", "failed");
+ return false;
+ }
+ std::stringstream ss;
+ ss << n_pop << " assemblies";
+ Record_("DISCHARGE", ss.str());
+ spent_inv.Push(core.PopN(n_pop));
+
+ std::map spent_mats;
+ for (std::string commod : out_commods) {
+ spent_mats = PeekSpent_();
+ cyclus::toolkit::MatVec mats = spent_mats[commod];
+ double total_spent = 0.;
+ for (cyclus::Material::Ptr mat_ptr : mats) {
+ total_spent += mat_ptr->quantity();
+ }
+ cyclus::toolkit::RecordTimeSeries("supply"+commod, this,
+ total_spent);
+ }
+ return true;
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+bool GprReactor::Retired_() {
+ return exit_time() != -1 && context()->time() > exit_time();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+std::map GprReactor::PeekSpent_() {
+ std::map mapped;
+ cyclus::toolkit::MatVec mats = spent_inv.PopN(spent_inv.count());
+ spent_inv.Push(mats);
+ for (cyclus::Material::Ptr mat_ptr : mats) {
+ std::string commod = OutCommod_(mat_ptr);
+ mapped[commod].push_back(mat_ptr);
+ }
+ return mapped;
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+std::map GprReactor::PopSpent_() {
+ cyclus::toolkit::MatVec mats = spent_inv.PopN(spent_inv.count());
+ std::map mapped;
+ for (int i = 0; i < mats.size(); ++i) {
+ std::string commod = OutCommod_(mats[i]);
+ mapped[commod].push_back(mats[i]);
+ }
+
+ // Reverse to ensure oldest assemblies are traded away first.
+ // TODO check this in cycamore Reactor, unsure why this is needed.
+ std::map::iterator it;
+ for (it = mapped.begin(); it != mapped.end(); ++it) {
+ std::reverse(it->second.begin(), it->second.end());
+ }
+ return mapped;
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::IndexRes_(cyclus::Resource::Ptr m, std::string incommod) {
+ for (int i = 0; i < in_commods.size(); ++i) {
+ if (in_commods[i] == incommod) {
+ res_indexes[m->obj_id()] = i;
+ return;
+ }
+ }
+ throw cyclus::ValueError(
+ "misoenrichment::GprReactor - received unsupported incommod material.");
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::Load_() {
+ int n_load = std::min(n_assem_core - core.count(), fresh_inv.count());
+ if (n_load == 0) {
+ return;
+ }
+ std::stringstream ss;
+ ss << n_load << " assemblies";
+ Record_("LOAD", ss.str());
+ core.Push(fresh_inv.PopN(n_load));
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::PushSpent_(std::map mats) {
+ std::map::iterator it;
+ for (it = mats.begin(); it != mats.end(); ++it) {
+ // Undo reverse in PopSpent to ensure oldest assemblies come out first.
+ // TODO check this in cycamore Reactor, unsure why this is needed.
+ std::reverse(it->second.begin(), it->second.end());
+ spent_inv.Push(it->second);
+ }
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::Record_(std::string name, std::string val) {
+ context()->NewDatum("ReactorEvents")
+ ->AddVal("AgentId", id())
+ ->AddVal("Time", context()->time())
+ ->AddVal("Event", name)
+ ->AddVal("Value", val)
+ ->Record();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::RecordPosition_() {
+ context()->NewDatum("AgentPosition")
+ ->AddVal("Spec", this->spec())
+ ->AddVal("Prototype", this->prototype())
+ ->AddVal("AgentId", id())
+ ->AddVal("Latitude", latitude)
+ ->AddVal("Longitude", longitude)
+ ->Record();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::RecordSideProduct_(bool is_producing) {
+ if (is_hybrid) {
+ double value;
+ for (int i = 0; i < side_products.size(); ++i) {
+ if (is_producing) {
+ value = side_product_quantity[i];
+ }
+ else {
+ value = 0;
+ }
+ context()->NewDatum("ReactorSideProducts")
+ ->AddVal("AgentId", id())
+ ->AddVal("Time", context()->time())
+ ->AddVal("Product", side_products[i])
+ ->AddVal("Value", value)
+ ->Record();
+ }
+ }
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::Transmute_() { Transmute_(n_assem_batch); }
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::Transmute_(int n_assem) {
+ cyclus::toolkit::MatVec old = core.PopN(std::min(n_assem, core.count()));
+ core.Push(old);
+ if (core.count() > old.size()) {
+ // Rotate the untransmuted materials back to the back of the buffer.
+ core.Push(core.PopN(core.count() - old.size()));
+ }
+ std::stringstream ss;
+ ss << old.size() << " assemblies";
+ Record_("TRANSMUTE", ss.str());
+
+ cyclus::PyStart();
+ int python_exit_code = 0;
+ python_exit_code += PyRun_SimpleString("import spentfuelgpr");
+ for (int i = 0; i < old.size(); ++i) {
+ CompositionToOutFile_(old[i]->comp(), false);
+ python_exit_code += PyRun_SimpleString("spentfuelgpr.predict()");
+ cyclus::Composition::Ptr spent_fuel_comp = ImportSpentFuelComposition_(
+ old[i]->quantity());
+ old[i]->Transmute(spent_fuel_comp);
+ }
+ if (python_exit_code != 0) {
+ throw cyclus::Error("Execution of Python code in Gpr::ReactorTick "
+ "unsuccessful!");
+ }
+
+ // Having finished the GPR calculations, delete the .json files to prevent
+ // cluttering up the working directory.
+ if (std::remove(in_fname.c_str()) != 0) {
+ std::stringstream msg;
+ msg << "Error deleting file '" << in_fname << "'.";
+ throw cyclus::IOError(msg.str());
+ }
+ if (std::remove(out_fname.c_str()) != 0) {
+ std::stringstream msg;
+ msg << "Error deleting file '" << out_fname << "'.";
+ throw cyclus::IOError(msg.str());
+ }
+ // TODO the code below can be deleted UNLESS it turns out that the GPR
+ // predictions are computationally significant and that they are causing a
+ // bottleneck. If this is the case, then do something along the lines shown
+ // below.
+ /*
+ bool same_composition = true;
+ cyclus::Composition::Ptr previous_comp = old[0]->comp();
+ for (cyclus::Material::Ptr mat : old) {
+ it (!cyclus::compmath::AlmostEq(mat->comp(), previous_comp)) {
+ same_composition = false;
+ break;
+ }
+ }
+ */
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactor::CompositionToOutFile_(cyclus::Composition::Ptr comp,
+ bool delete_outfile) {
+ nlohmann::json json_object;
+
+ // TODO check if GPRs use mass or atom percent
+ cyclus::CompMap cm = comp->atom();
+ cyclus::compmath::Normalize(&cm);
+ cyclus::CompMap::iterator compmap_it;
+ // Loop over each isotope in composition.
+ for (compmap_it = cm.begin(); compmap_it != cm.end(); ++compmap_it) {
+ std::set::iterator isotope_it = permitted_fresh_fuel_comps.find(
+ compmap_it->first);
+ // Ensure that the fresh fuel is solely composed of isotopes taken into
+ // account by the GPR, if not, throw an error.
+ if (isotope_it == permitted_fresh_fuel_comps.end()) {
+ std::stringstream msg;
+ msg << "GprReactor fuel must be composed of (some or all of) the "
+ "following isotopes: ";
+ for (const int& iso : permitted_fresh_fuel_comps) {
+ msg << iso << " ";
+ }
+ msg << "\n";
+ throw cyclus::ValueError(msg.str());
+ }
+ std::string nuc_id = std::to_string(compmap_it->first);
+ double fraction = compmap_it->second;
+ json_object["fresh_fuel_composition"][nuc_id] = fraction;
+ }
+
+ // TODO also pass `relevant_spent_fuel_comps` to tell GPR which isotopes to
+ // reconstruct?
+ uint64_t irradiation_time;
+ if (context() != NULL) {
+ const long seconds_per_day = 60 * 60 * 24;
+ irradiation_time = cycle_time * context()->sim_info().dt / seconds_per_day;
+ } else {
+ // kDefaultTimeStepDur defined in cyclus/src/context.h as duration in
+ // seconds of 1/12 of a year, like so:
+ // const uint64_t kDefaultTimeStepDur = 2629846;
+ irradiation_time = cycle_time * kDefaultTimeStepDur;
+ }
+
+ // TODO Update burnup calculation below.
+ // TODO This is strictly speaking not correct in the case of a reactor using
+ // for example uranium dioxide (UO2) as fuel. In that case, the denonimator
+ // would consist only of the mass of uranium, excluding the oxygen mass.
+ double burnup = power_output * irradiation_time / n_assem_core
+ / assem_size; // in MWd/kg
+
+ json_object["temperature"] = temperature; // in K
+ json_object["power_output"] = power_output; // in MWth
+ json_object["irradiation_time"] = irradiation_time; // in days
+ json_object["burnup"] = burnup; // Save JSON output in output file.
+ std::ofstream file(out_fname, std::ofstream::out | std::ofstream::trunc);
+ file << std::setw(2) << json_object << "\n";
+ file.close();
+
+ // Deleting the output file does not make sense except for unit tests to
+ // prevent cluttering up the working directory where the unit tests are
+ // performed.
+ if (delete_outfile && (std::remove(out_fname.c_str()) != 0)) {
+ std::stringstream msg;
+ msg << "Error deleting file '" << out_fname << "'.\n";
+ throw cyclus::IOError(msg.str());
+ }
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+cyclus::Composition::Ptr GprReactor::ImportSpentFuelComposition_(double qty) {
+ using cyclus::Composition;
+
+ nlohmann::json json_object;
+
+ // Dump file contents into json object.
+ std::ifstream file(in_fname, std::ifstream::in);
+ if (!file.is_open()) {
+ std::stringstream msg;
+ msg << "Cannot find file '" << in_fname << "'";
+ throw cyclus::IOError(msg.str());
+ }
+ file >> json_object;
+ file.close();
+
+ try {
+ json_object.at("spent_fuel_composition");
+ } catch (const nlohmann::detail::out_of_range& e) {
+ std::stringstream msg;
+ msg << "Cannot find key 'spent_fuel_composition' in '" << in_fname << "'."
+ << "\nException id: " << e.id;
+ throw cyclus::IOError(msg.str());
+ }
+
+ cyclus::CompMap cm;
+ // This variable is the ratio of the material to be transmuted ('qty' kg) to
+ // the mass of the full core.
+ const double fraction_of_core = qty / (n_assem_core * assem_size);
+ // 'mass' is the absolute mass of the isotope in question in a full reactor
+ // core after one irradiation period.
+ double mass;
+ double sum = 0;
+
+ for (const int& nuc_id : relevant_spent_fuel_comps) {
+ try {
+ // Convert NucID to a human-readable string as used in the json file, for
+ // example: '922350001' is converted to 'U235M'.
+ std::string key = pyne::nucname::name(nuc_id);
+ mass = json_object.at("spent_fuel_composition").at(key);
+ } catch (const nlohmann::detail::out_of_range& e) {
+ continue;
+ }
+ cm[nuc_id] = mass * fraction_of_core;
+ sum += mass * fraction_of_core;
+ }
+ if (cyclus::AlmostEq(sum, 0.)) {
+ // This error is thrown if no isotopes contained in
+ // `relevant_spent_fuel_comps` are found in the spent fuel composition file.
+ // Notably, this prevents the program to continue if the file were empty.
+ throw cyclus::ValueError("No relevant isotopes found in the spent fuel!\n");
+ }
+ // All isotopes part of the spent fuel but not calculated by the Gpr (i.e.,
+ // all isotopes not part of 'relevant_spent_fuel_comps') are `represented' by
+ // hydrogen (H1). This is obviously not correct, but we are not interested in
+ // this part of the spent fuel so it should be alright.
+ if (!cyclus::AlmostEq(qty, sum)) {
+ cm[10010000] = qty - sum;
+ }
+
+ // TODO check if Gpr uses atom or mass fraction
+ Composition::Ptr spent_fuel_comp = Composition::CreateFromMass(cm);
+ return spent_fuel_comp;
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+std::string GprReactor::OutCommod_(cyclus::Material::Ptr m) {
+ int i = res_indexes[m->obj_id()];
+ if (i >= out_commods.size()) {
+ throw cyclus::KeyError("misoenrichment::GprReactor - no outcommod for material object");
+ }
+ return out_commods[i];
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// WARNING! Do not change the following this function!!! This enables your
+// archetype to be dynamically loaded and any alterations will cause your
+// archetype to fail.
+extern "C" cyclus::Agent* ConstructGprReactor(cyclus::Context* ctx) {
+ return new GprReactor(ctx);
+}
+
+} // namespace misoenrichment
diff --git a/src/gpr_reactor.h b/src/gpr_reactor.h
new file mode 100644
index 0000000..517d6f8
--- /dev/null
+++ b/src/gpr_reactor.h
@@ -0,0 +1,295 @@
+#ifndef MISOENRICHMENT_SRC_GPR_REACTOR_H_
+#define MISOENRICHMENT_SRC_GPR_REACTOR_H_
+
+#include
+#include
+
+#include "cyclus.h"
+
+// Future changes relating to the implementation of Antonio's GPRs are marked
+// with the following comment:
+// TODO ANTONIO GPR
+//
+// TODO list:
+// - check line 49 in .cc file. Ensure that `unique_out_commods.empty()`
+// evaluates to `true`, else, the set gets not filled!
+// - no implementation of `Reactor::fuel_pref(cyclus::Material::Ptr)` because
+// it is not used in the reactor class.
+// - think about the weird (?) decommissioning behaviour of the cycamore
+// archetype and whether or not I should use it as well
+// - check if GPRs use mass or molar fractions?
+
+namespace misoenrichment {
+
+class GprReactor : public cyclus::Facility, public cyclus::toolkit::Position {
+ public:
+ /// Constructor for GprReactor Class
+ /// @param ctx the cyclus context for access to simulation-wide parameters
+ explicit GprReactor(cyclus::Context* ctx);
+
+ ~GprReactor();
+
+ friend class GprReactorTest;
+
+ /// The Prime Directive
+ /// Generates code that handles all input file reading and restart operations
+ /// (e.g., reading from the database, instantiating a new object, etc.).
+ /// @warning The Prime Directive must have a space before it!
+
+ #pragma cyclus
+
+ #pragma cyclus note {"doc": "A reactor facility that calculates its spent " \
+ "fuel composition using Gaussian process " \
+ "regression."}
+
+ bool CheckDecommissionCondition();
+ std::set::Ptr> GetMatlBids(
+ cyclus::CommodMap::type& commod_requests);
+ std::set::Ptr> GetMatlRequests();
+ std::string str();
+ void AcceptMatlTrades(
+ const std::vector,
+ cyclus::Material::Ptr> >& responses);
+ void EnterNotify();
+ void GetMatlTrades(
+ const std::vector >& trades,
+ std::vector,
+ cyclus::Material::Ptr> >& responses);
+ void Tick();
+ void Tock();
+
+ private:
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // Fuel commodities, recipes and preferences
+ #pragma cyclus var { \
+ "uitype": ["oneormore", "incommodity"], \
+ "doc": "An ordered list of the input commodities (fresh fuel)." \
+ }
+ std::vector in_commods;
+
+ #pragma cyclus var { \
+ "uitype": ["oneormore", "outcommodity"], \
+ "doc": "An ordered list of the (spent) output commodities." \
+ }
+ std::vector out_commods;
+
+ #pragma cyclus var { \
+ "uitype": ["oneormore", "inrecipe"], \
+ "doc": "An ordered list of the input (fresh fuel) recipes." \
+ }
+ std::vector in_recipes;
+
+ #pragma cyclus var { \
+ "default": [], \
+ "doc": "The preference for each input fuel type, in the same order as " \
+ "the input commodities. If no preferences are specified, then the " \
+ "same default value is used for all fuel requests." \
+ }
+ std::vector fuel_prefs;
+
+ std::set unique_out_commods;
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // Core specifics
+ #pragma cyclus var { \
+ "default": 3, \
+ "doc": "Number of assemblies in a full core. This value is used in " \
+ "combination with the assembly mass, the irradiation time and the " \
+ "thermal power output to calculate the specific burnup needed for " \
+ "the Gaussian process regression." \
+ }
+ int n_assem_core;
+
+ #pragma cyclus var { \
+ "doc": "Number of assemblies constituting a single batch. This is the " \
+ "number of assemblies that gets discharged from the core after " \
+ "one complete irradiation cycle." \
+ }
+ int n_assem_batch;
+
+ #pragma cyclus var { \
+ "uitype": "range", \
+ "range": [1., 1e5], \
+ "doc": "The mass of one assembly in kg. This value is used in " \
+ "combination with the number of assemblies in a full core, the " \
+ "irradiation time and the thermal power output to calculate the " \
+ "specific burnup needed for the Gaussian process regression." \
+ }
+ double assem_size;
+
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // Irradiation cycle parameters
+ #pragma cyclus var { \
+ "default": 0, \
+ "doc": "If true, the archetype transmutes all assemblies upon " \
+ "decommissioning. If false, the archetype only transmutes half." \
+ }
+ bool decom_transmute_all;
+
+ #pragma cyclus var { \
+ "default": 12, \
+ "doc": "The duration of one complete irradiation cycle excluding the " \
+ "refuelling, in units of simulation time steps. This value is " \
+ "used in combination with the core mass and the thermal power " \
+ "output to calculate the specific burnup needed for the Gaussian " \
+ "process regression." \
+ }
+ int cycle_time;
+
+ #pragma cyclus var { \
+ "default": 1, \
+ "doc": "The duration of an entire refuelling period, i.e., the minimum " \
+ "time between the end of a cycle and the start of the next cycle." \
+ "In units of simulation time steps." \
+ }
+ int refuel_time;
+
+ #pragma cyclus var { \
+ "default": 0, \
+ "doc": "The number of time steps since the start of the last cycle. Only " \
+ "set this variable if you know what you are doing." \
+ }
+ int cycle_step;
+
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // Side product settings
+ #pragma cyclus var { \
+ "default": 1, \
+ "internal": True, \
+ "doc": "If the reactor produces side products, then this variable is set " \
+ "to true." \
+ }
+ bool is_hybrid;
+
+ #pragma cyclus var { \
+ "default": [], \
+ "doc": "An ordered vector containing the side products produced during " \
+ "normal reactor operation." \
+ }
+ std::vector side_products;
+
+ #pragma cyclus var { \
+ "default": [], \
+ "doc": "An ordered vector containing the quantities of the produced side " \
+ "products. The vector entries must be of type 'double'." \
+ }
+ std::vector side_product_quantity;
+
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // Other variables
+ #pragma cyclus var { \
+ "default": 0, \
+ "uitype": "range", \
+ "range": [0, 10000], \
+ "units": "MWth", \
+ "doc": "The amount of thermal power that the reactor produces during " \
+ "operation. This value is used in combination with the core mass " \
+ "and the irradiation time to calculate the specific burnup needed " \
+ "for the Gaussian process regression." \
+ }
+ double power_output;
+
+ #pragma cyclus var { \
+ "default": 350, \
+ "units": "K", \
+ "doc": "The reactor moderator temperature." \
+ }
+ double temperature;
+
+ // This variable is internal only and true if fuel has already been discharged
+ // this cycle.
+ #pragma cyclus var { \
+ "default": 0, \
+ "internal": True, \
+ "doc": "This variable should NEVER be set manually." \
+ }
+ bool discharged;
+
+ #pragma cyclus var { \
+ "default": 0, \
+ "uitype": "range", \
+ "range": [0, 1000000000], \
+ "doc": "Maximum number of fresh fuel assemblies that are stored in " \
+ "the facility if available." \
+ }
+ int n_assem_fresh;
+
+ #pragma cyclus var { \
+ "default": 1000000000, \
+ "uitype": "range", \
+ "range": [0, 1000000000], \
+ "doc": "Maximum number of spent fuel assemblies that can be stored in " \
+ "the facility before reactor operation stalls." \
+ }
+ int n_assem_spent;
+
+ // This variable should be hidden/unavailable in ui. Maps resource object
+ // id's to the index for the incommod through which they were received.
+ #pragma cyclus var { \
+ "default": {}, \
+ "internal": True, \
+ "doc": "This variable should NEVER be set manually." \
+ }
+ std::map res_indexes;
+
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // Material buffers
+ # pragma cyclus var {"capacity": "n_assem_fresh * assem_size"}
+ cyclus::toolkit::ResBuf fresh_inv;
+ # pragma cyclus var {"capacity": "n_assem_core * assem_size"}
+ cyclus::toolkit::ResBuf core;
+ # pragma cyclus var {"capacity": "n_assem_spent * assem_size"}
+ cyclus::toolkit::ResBuf spent_inv;
+
+ // This set contains all nuc ids that may be part of fresh fuel. So far,
+ // they are limited to U235 and U238, however this list will expand in the
+ // future. Notably, other U isotopes will be included (probably U232 up to
+ // U236) and possibly Pu as well.
+ const std::set permitted_fresh_fuel_comps;
+ // This set contains all nuc ids of isotopes in spent fuel that we are
+ // interested in and that the GPRs calculate.
+ const std::set relevant_spent_fuel_comps;
+
+ // Filenames of files used to pass arguments and results between the Python
+ // file and the C++ Cyclus archetype.
+ const std::string out_fname;
+ const std::string in_fname;
+
+ // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ // Coordinates
+ #pragma cyclus var { \
+ "default": 0.0, \
+ "uilabel": "Geographical latitude in degrees as a double", \
+ "doc": "Latitude of the agent's geographical position. The value " \
+ "should be expressed in degrees as a double." \
+ }
+ double latitude;
+
+ #pragma cyclus var { \
+ "default": 0.0, \
+ "uilabel": "Geographical longitude in degrees as a double", \
+ "doc": "Longitude of the agent's geographical position. The value " \
+ "should be expressed in degrees as a double." \
+ }
+ double longitude;
+
+ bool Discharge_();
+ bool Retired_();
+ cyclus::Composition::Ptr ImportSpentFuelComposition_(double qty);
+ std::map PeekSpent_();
+ std::map PopSpent_();
+ std::string OutCommod_(cyclus::Material::Ptr m);
+ void CompositionToOutFile_(cyclus::Composition::Ptr comp,
+ bool delete_outfile = false);
+ void IndexRes_(cyclus::Resource::Ptr m, std::string incommod);
+ void Load_();
+ void Record_(std::string name, std::string val);
+ void RecordPosition_();
+ void RecordSideProduct_(bool produce);
+ void PushSpent_(std::map mats);
+ void Transmute_();
+ void Transmute_(int n_assem);
+};
+
+} // namespace misoenrichment
+
+#endif // MISOENRICHMENT_SRC_GPR_REACTOR_H_
diff --git a/src/gpr_reactor_tests.cc b/src/gpr_reactor_tests.cc
new file mode 100644
index 0000000..3bcf291
--- /dev/null
+++ b/src/gpr_reactor_tests.cc
@@ -0,0 +1,1010 @@
+#include "gpr_reactor_tests.h"
+
+#include
+
+#include "agent_tests.h"
+#include "context.h"
+#include "facility_tests.h"
+#include "pyhooks.h"
+
+#include
+
+#include "miso_helper.h"
+
+namespace misoenrichment {
+namespace gpr_reactor_test {
+
+cyclus::Composition::Ptr valid_composition() {
+ cyclus::CompMap valid_cm;
+ valid_cm[922350000] = 1.5;
+ valid_cm[922380000] = 98.5;
+
+ return cyclus::Composition::CreateFromMass(valid_cm);
+}
+
+cyclus::Composition::Ptr invalid_composition() {
+ cyclus::CompMap invalid_cm;
+ invalid_cm[922350000] = 1.5;
+ invalid_cm[922380000] = 98.5;
+ invalid_cm[10010000] = 1;
+
+ return cyclus::Composition::CreateFromMass(invalid_cm);
+}
+
+// The compositions below are taken from CNERG's cycamore module, see
+// https://github.com/cyclus/cycamore/blob/master/src/reactor_tests.cc
+// Copyright under BSD-3 License, University of Wisconsin Computational Nuclear
+// Engineering Research Group, 2010-2016.
+cyclus::Composition::Ptr c_uox() {
+ cyclus::CompMap m;
+ m[pyne::nucname::id("u235")] = 0.04;
+ m[pyne::nucname::id("u238")] = 0.96;
+ return cyclus::Composition::CreateFromMass(m);
+};
+
+cyclus::Composition::Ptr c_mox() {
+ cyclus::CompMap m;
+ m[pyne::nucname::id("u235")] = .7;
+ m[pyne::nucname::id("u238")] = 100;
+ m[pyne::nucname::id("pu239")] = 3.3;
+ return cyclus::Composition::CreateFromMass(m);
+};
+
+cyclus::Composition::Ptr c_spentuox() {
+ cyclus::CompMap m;
+ m[pyne::nucname::id("u235")] = .8;
+ m[pyne::nucname::id("u238")] = 100;
+ m[pyne::nucname::id("pu239")] = 1;
+ return cyclus::Composition::CreateFromMass(m);
+};
+
+cyclus::Composition::Ptr c_spentmox() {
+ cyclus::CompMap m;
+ m[pyne::nucname::id("u235")] = .2;
+ m[pyne::nucname::id("u238")] = 100;
+ m[pyne::nucname::id("pu239")] = .9;
+ return cyclus::Composition::CreateFromMass(m);
+};
+
+cyclus::Composition::Ptr c_water() {
+ cyclus::CompMap m;
+ m[pyne::nucname::id("O16")] = 1;
+ m[pyne::nucname::id("H1")] = 2;
+ return cyclus::Composition::CreateFromAtom(m);
+};
+
+} // namespace gpr_reactor_test
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Functions handling the initiation and ending of the tests.
+void GprReactorTest::SetUp() {
+ cyclus::PyStart();
+
+ facility = new GprReactor(tc.get());
+ InitParameters();
+ SetUpGprReactor();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactorTest::InitParameters() {
+ in_commods = std::vector(1, "fresh_fuel");
+ out_commods = std::vector(1, "spent_fuel");
+ in_recipes = std::vector(1, "fresh_fuel_recipe");
+ fuel_prefs = std::vector(1, 1);
+ n_assem_core = 1;
+ n_assem_batch = 1;
+ assem_size = 110820;
+ n_assem_fresh = 1;
+ n_assem_spent = 1;
+ latitude = 0.;
+ longitude = 0.;
+ decom_transmute_all = false;
+ cycle_time = 3;
+ refuel_time = 1;
+ temperature = 350;
+ power_output = 2400;
+
+ /*
+ // The calculations below need to be done because of the way the SRS kernel
+ // handles the power output (namely as a power density).
+ const int n_assemblies_tot = 600;
+ const int n_assemblies_model = 18;
+ const double feet_to_m = 0.3048; // in m * ft^-1
+ double assembly_length = 12; // in ft
+ assembly_length *= feet_to_m; // in m
+
+ power_output *= n_assemblies_model / n_assemblies_tot / assembly_length;
+ */
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactorTest::SetUpGprReactor() {
+ using cyclus::Material;
+
+ facility->in_commods = in_commods;
+ facility->out_commods = out_commods;
+ facility->in_recipes = in_recipes;
+ facility->fuel_prefs = fuel_prefs;
+ facility->n_assem_core = n_assem_core;
+ facility->n_assem_batch = n_assem_batch;
+ facility->assem_size = assem_size;
+ facility->n_assem_fresh = n_assem_fresh;
+ facility->n_assem_spent = n_assem_spent;
+ facility->latitude = latitude;
+ facility->longitude = longitude;
+ facility->decom_transmute_all = decom_transmute_all;
+ facility->cycle_time = cycle_time;
+ facility->refuel_time = refuel_time;
+ facility->power_output = power_output;
+ facility->temperature = temperature;
+
+ for (int i = 0; i < n_assem_core; ++i) {
+ cyclus::Material::Ptr mat = cyclus::Material::CreateUntracked(assem_size,
+ gpr_reactor_test::valid_composition());
+ facility->core.Push(mat);
+ }
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+void GprReactorTest::TearDown() {
+ delete facility;
+ cyclus::PyStop();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Unit tests
+TEST_F(GprReactorTest, ExportCompositions) {
+ using cyclus::Composition;
+
+ Composition::Ptr invalid_comp = gpr_reactor_test::invalid_composition();
+ EXPECT_THROW(DoCompositionToOutFile(invalid_comp, true),
+ cyclus::ValueError);
+
+ Composition::Ptr valid_comp = gpr_reactor_test::valid_composition();
+ EXPECT_NO_THROW(DoCompositionToOutFile(valid_comp, true));
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+TEST_F(GprReactorTest, ImportCompositions) {
+ using cyclus::Composition;
+
+ Composition::Ptr valid_comp = gpr_reactor_test::valid_composition();
+ cyclus::CompMap valid_cm = valid_comp->mass();
+
+ // Prepare import composition-tests.
+ nlohmann::json json_object;
+ cyclus::compmath::Normalize(&valid_cm);
+ for (const std::pair& x : valid_cm) {
+ // Convert the nuc id to a name (e.g., 922350001 --> U235M) to match the
+ // output of the `spentfuelgpr` program and multiply the mass fraction with
+ // the total mass of the core.
+ std::string nuc_name = pyne::nucname::name(x.first);
+ double mass = n_assem_core * assem_size * x.second;
+ json_object["spent_fuel_composition"][nuc_name] = mass;
+ }
+ std::ofstream file("gpr_reactor_spent_fuel_composition.json",
+ std::ofstream::out | std::ofstream::trunc);
+ file << json_object;
+ file.close();
+
+ // Perform and test the imports
+ ASSERT_NO_THROW(DoImportSpentFuelComposition(n_assem_core * assem_size));
+ cyclus::CompMap return_cm = DoImportSpentFuelComposition(
+ n_assem_core * assem_size)->mass();
+ cyclus::compmath::Normalize(&return_cm);
+ EXPECT_TRUE(cyclus::compmath::AlmostEq(valid_cm, return_cm, kEpsCompMap));
+
+ nlohmann::json json_object2;
+ json_object2["test"] = -1;
+ file.open("gpr_reactor_spent_fuel_composition.json",
+ std::ofstream::out | std::ofstream::trunc);
+ file << json_object2;
+ file.close();
+ EXPECT_THROW(DoImportSpentFuelComposition(n_assem_core * assem_size / 2),
+ cyclus::IOError);
+ std::remove("gpr_reactor_spent_fuel_composition.json");
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+TEST_F(GprReactorTest, SpentFuelWaste) {
+ using cyclus::Composition;
+
+ Composition::Ptr valid_comp = gpr_reactor_test::valid_composition();
+ cyclus::CompMap valid_cm = valid_comp->mass();
+
+ // Prepare import composition-tests.
+ nlohmann::json json_object;
+ cyclus::compmath::Normalize(&valid_cm);
+ for (const std::pair& x : valid_cm) {
+ // Convert the nuc id to a name (e.g., 922350001 --> U235M) to match the
+ // output of the `spentfuelgpr` program and multiply the mass fraction with
+ // the total mass of the core.
+ // In this test, multiply by half a core (as opposed to the test in
+ // `GprReactorTest.ImportCompositions`) because we want to test if half of
+ // the spent fuel gets set to hydrogen.
+ std::string nuc_name = pyne::nucname::name(x.first);
+ double mass = n_assem_core * assem_size * x.second / 2.;
+ json_object["spent_fuel_composition"][nuc_name] = mass;
+ }
+ std::ofstream file("gpr_reactor_spent_fuel_composition.json",
+ std::ofstream::out | std::ofstream::trunc);
+ file << json_object;
+ file.close();
+
+ // Perform and test the import.
+ cyclus::CompMap return_cm = DoImportSpentFuelComposition(
+ n_assem_core * assem_size)->mass();
+ cyclus::compmath::Normalize(&return_cm);
+ valid_cm[10010000] = 1.;
+ cyclus::compmath::Normalize(&valid_cm);
+ EXPECT_TRUE(cyclus::compmath::AlmostEq(valid_cm, return_cm, kEpsCompMap));
+
+ std::remove("gpr_reactor_spent_fuel_composition.json");
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+TEST_F(GprReactorTest, TransmuteFuel) {
+ DoTransmute();
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Below are unit tests taken from CNERG's cycamore module, see
+// https://github.com/cyclus/cycamore/blob/master/src/reactor_tests.cc
+// Copyright under BSD-3 License, University of Wisconsin Computational Nuclear
+// Engineering Research Group, 2010-2016.
+
+// Test that with a zero refuel_time and a zero capacity fresh fuel buffer
+// (the default), fuel can be ordered and the cycle started with no time step
+// delay.
+TEST_F(GprReactorTest, JustInTimeOrdering) {
+ std::string config =
+ " lwr_fresh "
+ " lwr_spent "
+ " enriched_u "
+ " waste "
+ " 1.0 "
+ ""
+ " 1 "
+ " 0 "
+ " 300 "
+ " 1 "
+ " 1 ";
+
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("enriched_u").Finalize();
+ sim.AddRecipe("lwr_fresh", gpr_reactor_test::c_uox());
+ sim.AddRecipe("lwr_spent", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("Transactions", NULL);
+ EXPECT_EQ(simdur, qr.rows.size()) << "failed to order+run on fresh fuel "
+ "inside 1 time step";
+}
+
+// tests that the correct number of assemblies are popped from the core each
+// cycle.
+TEST_F(GprReactorTest, BatchSizes) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 1 "
+ " 0 "
+ " 1 "
+ " 7 "
+ " 3 ";
+
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("Transactions", NULL);
+ // 7 for initial core, 3 per time step for each new batch for remainder
+ EXPECT_EQ(7+3*(simdur-1), qr.rows.size());
+}
+
+// tests that the refueling period between cycle end and start of the next
+// cycle is honored.
+TEST_F(GprReactorTest, RefuelTimes) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 4 "
+ " 3 "
+ " 1 "
+ " 1 "
+ " 1 ";
+
+ int simdur = 49;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("Transactions", NULL);
+ int cyclet = 4;
+ int refuelt = 3;
+ int n_assem_want = simdur/(cyclet+refuelt)+1; // +1 for initial core
+ EXPECT_EQ(n_assem_want, qr.rows.size());
+}
+
+
+// tests that a reactor decommissions on time without producing
+// power at the end of its lifetime.
+TEST_F(GprReactorTest, DecomTimes) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 2 "
+ " 2 "
+ " 1 "
+ " 3 "
+ " 1000 "
+ " 1 ";
+
+ int simdur = 12;
+ int lifetime = 7;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur,
+ lifetime);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ // operating for 2+2 months and shutdown for 2+1
+ int on_time = 4;
+ std::vector conds;
+ conds.push_back(cyclus::Cond("Value", "==", 1000));
+ cyclus::QueryResult qr = sim.db().Query("TimeSeriesPower", &conds);
+ EXPECT_EQ(on_time, qr.rows.size());
+
+ int off_time = 3;
+ conds.clear();
+ conds.push_back(cyclus::Cond("Value", "==", 0));
+ qr = sim.db().Query("TimeSeriesPower", &conds);
+ EXPECT_EQ(off_time, qr.rows.size());
+}
+
+
+// Tests if a reactor produces power at the time of its decommission
+// given a refuel_time of zero.
+TEST_F(GprReactorTest, DecomZeroRefuel) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 2 "
+ " 0 "
+ " 1 "
+ " 3 "
+ " 1000 "
+ " 1 ";
+
+ int simdur = 8;
+ int lifetime = 6;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur,
+ lifetime);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ // operating for 2+2 months and shutdown for 2+1
+ int on_time = 6;
+ std::vector conds;
+ conds.push_back(cyclus::Cond("Value", "==", 1000));
+ cyclus::QueryResult qr = sim.db().Query("TimeSeriesPower", &conds);
+ EXPECT_EQ(on_time, qr.rows.size());
+}
+
+// tests that new fuel is ordered immediately following cycle end - at the
+// start of the refueling period - not before and not after. - thie is subtly
+// different than RefuelTimes test and is not a duplicate of it.
+TEST_F(GprReactorTest, OrderAtRefuelStart) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 4 "
+ " 3 "
+ " 1 "
+ " 1 "
+ " 1 ";
+
+ int simdur = 7;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("Transactions", NULL);
+ int cyclet = 4;
+ int refuelt = 3;
+ int n_assem_want = simdur/(cyclet+refuelt)+1; // +1 for initial core
+ EXPECT_EQ(n_assem_want, qr.rows.size());
+}
+
+// tests that the reactor handles requesting multiple types of fuel correctly
+// - with proper inventory constraint honoring, etc.
+TEST_F(GprReactorTest, MultiFuelMix) {
+ std::string config =
+ " uox mox "
+ " spentuox spentmox "
+ " uox mox "
+ " waste waste "
+ ""
+ " 1 "
+ " 0 "
+ " 1 "
+ " 3 "
+ " 3 "
+ " 3 ";
+
+ // it is important that the sources have cumulative capacity greater than
+ // the reactor can take on a single time step - to test that inventory
+ // capacity constraints are being set properly. It is also important that
+ // each source is smaller capacity thatn the reactor orders on each time
+ // step to make it easy to compute+check the number of transactions.
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").capacity(2).Finalize();
+ sim.AddSource("mox").capacity(2).Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ sim.AddRecipe("mox", gpr_reactor_test::c_spentuox());
+ sim.AddRecipe("spentmox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("Transactions", NULL);
+ // +3 is for fresh fuel inventory
+ EXPECT_EQ(3*simdur+3, qr.rows.size());
+}
+
+// tests that the reactor halts operation when it has no more room in its
+// spent fuel inventory buffer.
+TEST_F(GprReactorTest, FullSpentInventory) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 1 "
+ " 0 "
+ " 1 "
+ " 1 "
+ " 1 "
+ " 3 ";
+
+ int simdur = 10;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("Transactions", NULL);
+ int n_assem_spent = 3;
+
+ // +1 is for the assembly in the core + the three in spent
+ EXPECT_EQ(n_assem_spent+1, qr.rows.size());
+}
+
+// tests that the reactor shuts down, ie., does not generate power, when the
+// spent fuel inventory is full and the core cannot be unloaded.
+TEST_F(GprReactorTest, FullSpentInventoryShutdown) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 1 "
+ " 0 "
+ " 1 "
+ " 1 "
+ " 1 "
+ " 1 "
+ " 100 ";
+
+ int simdur = 3;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("TimeSeriesPower", NULL);
+ EXPECT_EQ(0, qr.GetVal("Value", simdur - 1));
+
+}
+
+// tests that the reactor cycle is delayed as expected when it is unable to
+// acquire fuel in time for the next cycle start. This checks that after a
+// cycle is delayed past an original scheduled start time, as soon as enough
+// fuel is received, a new cycle pattern is established starting from the
+// delayed start time.
+TEST_F(GprReactorTest, FuelShortage) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 7 "
+ " 0 "
+ " 1 "
+ " 3 "
+ " 3 ";
+
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ // provide initial full batch
+ sim.AddSource("uox").lifetime(1).Finalize();
+ // provide partial batch post cycle-end
+ sim.AddSource("uox").start(9).lifetime(1).capacity(2).Finalize();
+ // provide remainder of batch much later
+ sim.AddSource("uox").start(15).Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ // check that we never got a full refueled batch during refuel period
+ std::vector conds;
+ conds.push_back(cyclus::Cond("Time", "<", 15));
+ cyclus::QueryResult qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(5, qr.rows.size());
+
+ // after being delayed past original scheduled start of new cycle, we got
+ // final assembly for core.
+ conds.clear();
+ conds.push_back(cyclus::Cond("Time", "==", 15));
+ qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(1, qr.rows.size());
+
+ // all during the next (delayed) cycle we shouldn't be requesting any new fuel
+ conds.clear();
+ conds.push_back(cyclus::Cond("Time", "<", 21));
+ qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(6, qr.rows.size());
+
+ // as soon as this delayed cycle ends, we should be requesting/getting 3 new
+ // batches
+ conds.clear();
+ conds.push_back(cyclus::Cond("Time", "==", 22));
+ qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(3, qr.rows.size());
+}
+
+// tests that discharged fuel is transmuted properly immediately at cycle end.
+TEST_F(GprReactorTest, DischargedFuelTransmute) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 4 "
+ " 3 "
+ " 1 "
+ " 1 "
+ " 1 ";
+
+ int simdur = 7;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddSink("waste").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ Composition::Ptr spentuox = gpr_reactor_test::c_spentuox();
+ sim.AddRecipe("spentuox", spentuox);
+ int id = sim.Run();
+
+ std::vector conds;
+ conds.push_back(cyclus::Cond("SenderId", "==", id));
+ int resid = sim.db().Query("Transactions", &conds).GetVal("ResourceId");
+ cyclus::Material::Ptr m = sim.GetMaterial(resid);
+ cyclus::toolkit::MatQuery mq(m);
+ EXPECT_EQ(spentuox->id(), m->comp()->id());
+ EXPECT_TRUE(mq.mass(942390000) > 0) << "transmuted spent fuel doesn't have Pu239";
+}
+
+// tests that spent fuel is offerred on correct commods according to the
+// incommod it was received on - esp when dealing with multiple fuel commods
+// simultaneously.
+TEST_F(GprReactorTest, SpentFuelProperCommodTracking) {
+ std::string config =
+ " uox mox "
+ " spentuox spentmox "
+ " uox mox "
+ " waste1 waste2 "
+ ""
+ " 1 "
+ " 0 "
+ " 1 "
+ " 3 "
+ " 3 ";
+
+ int simdur = 7;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").capacity(1).Finalize();
+ sim.AddSource("mox").capacity(2).Finalize();
+ sim.AddSink("waste1").Finalize();
+ sim.AddSink("waste2").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ sim.AddRecipe("mox", gpr_reactor_test::c_mox());
+ sim.AddRecipe("spentmox", gpr_reactor_test::c_spentmox());
+ int id = sim.Run();
+
+ std::vector conds;
+ conds.push_back(cyclus::Cond("SenderId", "==", id));
+ conds.push_back(cyclus::Cond("Commodity", "==", std::string("waste1")));
+ cyclus::QueryResult qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(simdur-1, qr.rows.size());
+
+ conds[1] = cyclus::Cond("Commodity", "==", std::string("waste2"));
+ qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(2*(simdur-1), qr.rows.size());
+}
+
+// The user can optionally omit fuel preferences. In the case where
+// preferences are adjusted, the ommitted preference vector must be populated
+// with default values - if it wasn't then preferences won't be adjusted
+// correctly and the reactor could segfault. Check that this doesn't happen.
+TEST_F(GprReactorTest, PrefChange) {
+ // it is important that the fuel_prefs not be present in the config below.
+ std::string config =
+ " lwr_fresh "
+ " lwr_spent "
+ " enriched_u "
+ " waste "
+ ""
+ " 1 "
+ " 0 "
+ " 300 "
+ " 1 "
+ " 1 "
+ ""
+ " 25 "
+ " enriched_u "
+ " -1 ";
+
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("enriched_u").Finalize();
+ sim.AddRecipe("lwr_fresh", gpr_reactor_test::c_uox());
+ sim.AddRecipe("lwr_spent", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("Transactions", NULL);
+ EXPECT_EQ(25, qr.rows.size()) << "failed to adjust preferences properly";
+}
+
+TEST_F(GprReactorTest, RecipeChange) {
+ // it is important that the fuel_prefs not be present in the config below.
+ std::string config =
+ " lwr_fresh "
+ " lwr_spent "
+ " enriched_u "
+ " waste "
+ ""
+ " 1 "
+ " 0 "
+ " 300 "
+ " 1 "
+ " 1 "
+ ""
+ " 25 35 "
+ " enriched_u enriched_u "
+ " water water "
+ " lwr_spent water ";
+
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("enriched_u").Finalize();
+ sim.AddSink("waste").Finalize();
+ sim.AddRecipe("lwr_fresh", gpr_reactor_test::c_uox());
+ sim.AddRecipe("lwr_spent", gpr_reactor_test::c_spentuox());
+ sim.AddRecipe("water", gpr_reactor_test::c_water());
+ int aid = sim.Run();
+
+ std::vector conds;
+ cyclus::QueryResult qr;
+
+ // check that received recipe is not water
+ conds.clear();
+ conds.push_back(cyclus::Cond("Time", "==", 24));
+ conds.push_back(cyclus::Cond("ReceiverId", "==", aid));
+ qr = sim.db().Query("Transactions", &conds);
+ cyclus::toolkit::MatQuery mq = cyclus::toolkit::MatQuery(
+ sim.GetMaterial(qr.GetVal("ResourceId")));
+
+ EXPECT_TRUE(0 < mq.qty());
+ EXPECT_TRUE(0 == mq.mass(pyne::nucname::id("H1")));
+
+ // check that received recipe changed to water
+ conds.clear();
+ conds.push_back(cyclus::Cond("Time", "==", 26));
+ conds.push_back(cyclus::Cond("ReceiverId", "==", aid));
+ qr = sim.db().Query("Transactions", &conds);
+ mq = cyclus::toolkit::MatQuery(sim.GetMaterial(qr.GetVal("ResourceId")));
+
+ EXPECT_TRUE(0 < mq.qty());
+ EXPECT_TRUE(0 < mq.mass(pyne::nucname::id("H1")));
+
+ // check that sent recipe is not water
+ conds.clear();
+ conds.push_back(cyclus::Cond("Time", "==", 34));
+ conds.push_back(cyclus::Cond("SenderId", "==", aid));
+ qr = sim.db().Query("Transactions", &conds);
+ mq = cyclus::toolkit::MatQuery(sim.GetMaterial(qr.GetVal("ResourceId")));
+
+ EXPECT_TRUE(0 < mq.qty());
+ EXPECT_TRUE(0 == mq.mass(pyne::nucname::id("H1")));
+
+ // check that sent recipe changed to water
+ conds.clear();
+ conds.push_back(cyclus::Cond("Time", "==", 36));
+ conds.push_back(cyclus::Cond("SenderId", "==", aid));
+ qr = sim.db().Query("Transactions", &conds);
+ mq = cyclus::toolkit::MatQuery(sim.GetMaterial(qr.GetVal("ResourceId")));
+
+ EXPECT_TRUE(0 < mq.qty());
+ EXPECT_TRUE(0 < mq.mass(pyne::nucname::id("H1")));
+}
+
+TEST_F(GprReactorTest, Retire) {
+ std::string config =
+ " lwr_fresh "
+ " lwr_spent "
+ " enriched_u "
+ " waste "
+ ""
+ " 7 "
+ " 0 "
+ " 300 "
+ " 1 "
+ " 3 "
+ " 1 "
+ " 1 "
+ "";
+
+ int dur = 50;
+ int life = 36;
+ int cycle_time = 7;
+ int refuel_time = 0;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, dur,
+ life);
+ sim.AddSource("enriched_u").Finalize();
+ sim.AddSink("waste").Finalize();
+ sim.AddRecipe("lwr_fresh", gpr_reactor_test::c_uox());
+ sim.AddRecipe("lwr_spent", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ int ncore = 3;
+ int nbatch = 1;
+
+ // reactor should stop requesting new fresh fuel as it approaches retirement
+ int nassem_recv =
+ static_cast(ceil(static_cast(life) / 7.0)) * nbatch +
+ (ncore - nbatch);
+
+ std::vector conds;
+ conds.push_back(cyclus::Cond("ReceiverId", "==", id));
+ cyclus::QueryResult qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(nassem_recv, qr.rows.size())
+ << "failed to stop ordering near retirement";
+
+ // reactor should discharge all fuel before/by retirement
+ conds.clear();
+ conds.push_back(cyclus::Cond("SenderId", "==", id));
+ qr = sim.db().Query("Transactions", &conds);
+ EXPECT_EQ(nassem_recv, qr.rows.size())
+ << "failed to discharge all material by retirement time";
+
+ // reactor should record power entry on the time step it retires if operating
+ int time_online = life / (cycle_time + refuel_time) * cycle_time
+ + std::min(life % (cycle_time + refuel_time), cycle_time);
+ conds.clear();
+ conds.push_back(cyclus::Cond("AgentId", "==", id));
+ conds.push_back(cyclus::Cond("Value", ">", 0));
+ qr = sim.db().Query("TimeSeriesPower", &conds);
+ EXPECT_EQ(time_online, qr.rows.size())
+ << "failed to generate power for the correct number of time steps";
+}
+
+TEST_F(GprReactorTest, PositionInitialize) {
+ std::string config =
+ " lwr_fresh "
+ " lwr_spent "
+ " enriched_u "
+ " waste "
+ " 1.0 "
+ ""
+ " 1 "
+ " 0 "
+ " 300 "
+ " 1 "
+ " 1 ";
+
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("enriched_u").Finalize();
+ sim.AddRecipe("lwr_fresh", gpr_reactor_test::c_uox());
+ sim.AddRecipe("lwr_spent", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("AgentPosition", NULL);
+ EXPECT_EQ(qr.GetVal("Latitude"), 0.0);
+ EXPECT_EQ(qr.GetVal("Longitude"), 0.0);
+}
+
+TEST_F(GprReactorTest, PositionInitialize2) {
+ std::string config =
+ " lwr_fresh "
+ " lwr_spent "
+ " enriched_u "
+ " waste "
+ " 1.0 "
+ ""
+ " 1 "
+ " 0 "
+ " 300 "
+ " 1 "
+ " 1 "
+ " 30.0 "
+ " 30.0 ";
+
+ int simdur = 50;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("enriched_u").Finalize();
+ sim.AddRecipe("lwr_fresh", gpr_reactor_test::c_uox());
+ sim.AddRecipe("lwr_spent", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ cyclus::QueryResult qr = sim.db().Query("AgentPosition", NULL);
+ EXPECT_EQ(qr.GetVal("Latitude"), 30.0);
+ EXPECT_EQ(qr.GetVal("Longitude"), 30.0);
+}
+
+TEST_F(GprReactorTest, ByProduct) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 1 "
+ " 1 "
+ " 1 "
+ " 7 "
+ " 3 "
+ ""
+ " process_heat "
+ " 10 ";
+
+ int simdur = 10;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+ std::vector conds;
+ // test if it produces side products only when reactor is running
+ int quantity = 10;
+ conds.push_back(cyclus::Cond("Value", "==", quantity));
+ cyclus::QueryResult qr = sim.db().Query("ReactorSideProducts", &conds);
+ EXPECT_EQ(5, qr.rows.size());
+
+ // test if it doesn't produce side products when reactor is refueling
+ conds.clear();
+ conds.push_back(cyclus::Cond("Value", "==", 0));
+ qr = sim.db().Query("ReactorSideProducts", &conds);
+ EXPECT_EQ(5, qr.rows.size());
+}
+
+TEST_F(GprReactorTest, MultipleByProduct) {
+ std::string config =
+ " uox "
+ " spentuox "
+ " uox "
+ " waste "
+ ""
+ " 1 "
+ " 1 "
+ " 1 "
+ " 7 "
+ " 3 "
+ ""
+ " process_heat water "
+ " 10 100 ";
+
+ int simdur = 10;
+ cyclus::MockSim sim(cyclus::AgentSpec(":cycamore:Reactor"), config, simdur);
+ sim.AddSource("uox").Finalize();
+ sim.AddRecipe("uox", gpr_reactor_test::c_uox());
+ sim.AddRecipe("spentuox", gpr_reactor_test::c_spentuox());
+ int id = sim.Run();
+
+
+ std::vector conds;
+ // test if it produces heat when reactor is running
+ int quantity = 10;
+ conds.push_back(cyclus::Cond("Product", "==", std::string("process_heat")));
+ conds.push_back(cyclus::Cond("Value", "==", quantity));
+ cyclus::QueryResult qr = sim.db().Query("ReactorSideProducts", &conds);
+ EXPECT_EQ(5, qr.rows.size());
+
+ // test if it produces water when reactor is running
+ conds.clear();
+ quantity = 100;
+ conds.push_back(cyclus::Cond("Product", "==", std::string("water")));
+ conds.push_back(cyclus::Cond("Value", "==", quantity));
+ qr = sim.db().Query("ReactorSideProducts", &conds);
+ EXPECT_EQ(5, qr.rows.size());
+
+ conds.clear();
+ conds.push_back(cyclus::Cond("Value", "==", 0));
+ qr = sim.db().Query("ReactorSideProducts", &conds);
+ EXPECT_EQ(10, qr.rows.size());
+}
+
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Below are standard tests introduced by cycstub. Probably not so useful but
+// they don't do harm, either.
+TEST_F(GprReactorTest, Print) {
+ EXPECT_NO_THROW(std::string s = facility->str());
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+TEST_F(GprReactorTest, Tick) {
+ EXPECT_NO_THROW(facility->Tick());
+}
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+TEST_F(GprReactorTest, Tock) {
+ EXPECT_NO_THROW(facility->Tock());
+}
+
+} // namespace misoenrichment
+
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Do Not Touch! Below section required for connection with Cyclus
+cyclus::Agent* GprReactorConstructor(cyclus::Context* ctx) {
+ return new misoenrichment::GprReactor(ctx);
+}
+// Required to get functionality in cyclus agent unit tests library
+#ifndef CYCLUS_AGENT_TESTS_CONNECTED
+int ConnectAgentTests();
+static int cyclus_agent_tests_connected = ConnectAgentTests();
+#define CYCLUS_AGENT_TESTS_CONNECTED cyclus_agent_tests_connected
+#endif // CYCLUS_AGENT_TESTS_CONNECTED
+INSTANTIATE_TEST_CASE_P(GprReactor, FacilityTests,
+ ::testing::Values(&GprReactorConstructor));
+INSTANTIATE_TEST_CASE_P(GprReactor, AgentTests,
+ ::testing::Values(&GprReactorConstructor));
diff --git a/src/gpr_reactor_tests.h b/src/gpr_reactor_tests.h
new file mode 100644
index 0000000..c104a84
--- /dev/null
+++ b/src/gpr_reactor_tests.h
@@ -0,0 +1,64 @@
+#ifndef MISOENRICHMENT_SRC_GPR_REACTOR_TESTS_H_
+#define MISOENRICHMENT_SRC_GPR_REACTOR_TESTS_H_
+
+#include
+
+#include "agent_tests.h"
+#include "gpr_reactor.h"
+
+namespace misoenrichment {
+
+class GprReactorTest : public ::testing::Test {
+ protected:
+ cyclus::TestContext tc;
+ GprReactor* facility;
+
+ virtual void SetUp();
+ virtual void TearDown();
+
+ void InitParameters();
+ void SetUpGprReactor();
+
+ bool decom_transmute_all;
+
+ double assem_size;
+ double latitude;
+ double longitude;
+ double power_output;
+ double temperature;
+
+ int n_assem_core;
+ int n_assem_batch;
+ int n_assem_fresh;
+ int n_assem_spent;
+ int cycle_time;
+ int refuel_time;
+
+ std::vector in_commods;
+ std::vector out_commods;
+ std::vector in_recipes;
+ std::vector fuel_prefs;
+
+ // The Do* functions are a hack: GprReactorTest is a friend class to
+ // GprReactor and it can therefore access private functions. However,
+ // all of the tests are sub-classes of the fixture and they cannot access
+ // the private functions, hence we need to use the Do* functions as an
+ // intermediary. See, e.g., 4th bullet point here:
+ // https://github.com/google/googletest/blob/master/googletest/docs/advanced.md#testing-private-code
+ //
+ // TODO implement functions below
+ inline cyclus::Composition::Ptr DoImportSpentFuelComposition(double qty) {
+ return facility->ImportSpentFuelComposition_(qty);
+ }
+ inline void DoCompositionToOutFile(cyclus::Composition::Ptr comp,
+ bool delete_txt) {
+ facility->CompositionToOutFile_(comp, delete_txt);
+ }
+ inline void DoTransmute() {
+ facility->Transmute_();
+ }
+};
+
+} // namespace misoenrichment
+
+#endif // MISOENRICHMENT_SRC_GPR_REACTOR_TESTS_H_