Merge pull request #6 from gwenchee/framework

First major ROLLO PR with source code, tests, and one example problem

README.md

@@ -1,2 +1,35 @@
-# generative-reactor-designs
-Repository to hold the framework and code to produce generative reactor designs. 
+# ROLLO (Reactor evOLutionary aLgorithm Optimizer)
+
+ROLLO is a Python package that applies evolutionary algorithm techniques to optimize nuclear reactor design.
+
+Documentation on the usage of ROLLO is hosted at: ## INSERT WEBSITE ##
+
+## Installation 
+`python -m pip install rollo`
+
+## Running ROLLO 
+`python -m rollo -i <input file> -c <checkpoint file>`
+
+Command line flags: 
+| Flag | Description | Mandatory ? |
+| ----------- | ----------- | ----------- |
+| -i | name of input file | Yes |
+| -c| name of checkpoint file | No |
+
+The checkpoint file holds the results from the ROLLO simulation and also acts 
+as a restart file. Thus, if a ROLLO simulation ends prematurely, the checkpoint 
+file can be used to restart the code from the most recent population and 
+continue the simulation.
+
+## To build documentation 
+`cd docs` 
+
+`make html`
+
+## To upload to PyPI
+`python3 -m build`
+
+`python3 -m twine upload dist/*`
+
+## License 
+ROLLO is distributed under the 3-Clause BSD License.

examples/fhr-slab/README.md

@@ -0,0 +1,8 @@
+**ROLLO input file and supporting OpenMC templating scripts for maximizing keff
+in a FHR fuel slab by varying the slab's TRISO fuel particle distribution.**
+
+`rollo_input.json`: ROLLO input file 
+
+`openmc_input.py`: OpenMC templating script 
+
+`constants.py`: Supporting functions for `openmcinp.py`

examples/fhr-slab/constants.py

@@ -0,0 +1,123 @@
+import openmc
+from numpy import pi
+
+# Constants
+T_r1 = 2135e-5
+T_r2 = 3135e-5
+T_r3 = 3485e-5
+T_r4 = 3835e-5
+T_r5 = 4235e-5
+T_pitch = 0.09266
+
+uoc_9 = openmc.Material()
+uoc_9.set_density("g/cc", 11)
+uoc_9.add_nuclide("U235", 2.27325e-3)
+uoc_9.add_nuclide("U238", 2.269476e-2)
+uoc_9.add_nuclide("O16", 3.561871e-2)
+uoc_9.add_nuclide("C0", 9.79714e-3)
+uoc_9.temperature = 1110
+uoc_9.volume = 4 / 3 * pi * (T_r1 ** 3) * 101 * 210 * 4 * 36
+
+por_c = openmc.Material()
+por_c.set_density("g/cc", 1)
+por_c.add_nuclide("C0", 5.013980e-2)
+por_c.temperature = 948
+
+si_c = openmc.Material()
+si_c.set_density("g/cc", 3.2)
+si_c.add_nuclide("Si28", 4.431240e-2)
+si_c.add_nuclide("Si29", 2.25887e-3)
+si_c.add_nuclide("Si30", 1.48990e-3)
+si_c.add_nuclide("C0", 4.806117e-2)
+si_c.temperature = 948
+
+graphite = openmc.Material()
+graphite.set_density("g/cc", 1.8)
+graphite.add_nuclide("C0", 9.025164e-2)
+graphite.temperature = 948
+
+triso_4_layers = openmc.Material()
+triso_4_layers.add_nuclide("C0", 0.06851594519357823)
+triso_4_layers.add_nuclide("Si28", 0.009418744960032735)
+triso_4_layers.add_nuclide("Si29", 0.00048013017638108395)
+triso_4_layers.add_nuclide("Si30", 0.0003166830980933728)
+triso_4_layers.set_density("sum")
+triso_4_layers.temperature = 948
+
+lm_graphite = openmc.Material()
+lm_graphite.set_density("g/cc", 1.8)
+lm_graphite.add_nuclide("C0", 9.025164e-2)
+lm_graphite.temperature = 948
+
+flibe = openmc.Material()
+flibe.set_density("g/cc", 1.95)
+flibe.add_nuclide("Li6", 1.383014e-6)
+flibe.add_nuclide("Li7", 2.37132e-2)
+flibe.add_nuclide("Be9", 1.18573e-2)
+flibe.add_nuclide("F19", 4.74291e-2)
+flibe.temperature = 948
+
+mats = openmc.Materials(
+    (uoc_9, por_c, si_c, graphite, lm_graphite, flibe, triso_4_layers)
+)
+
+# 4 layer triso
+two_spheres = [openmc.Sphere(r=r) for r in [T_r1, T_r5]]
+two_triso_cells = [
+    openmc.Cell(fill=uoc_9, region=-two_spheres[0]),
+    openmc.Cell(fill=triso_4_layers, region=+two_spheres[0] & -two_spheres[1]),
+    openmc.Cell(fill=lm_graphite, region=+two_spheres[1]),
+]
+two_triso_univ = openmc.Universe(cells=two_triso_cells)
+
+
+def create_prism(left, right, left_refl, right_refl):
+    if left_refl:
+        xplane_left = +openmc.XPlane(x0=left, boundary_type="reflective")
+    else:
+        xplane_left = +openmc.XPlane(x0=left)
+    if right_refl:
+        xplane_right = -openmc.XPlane(x0=right, boundary_type="reflective")
+    else:
+        xplane_right = -openmc.XPlane(x0=right)
+    prism = (
+        xplane_left &
+        xplane_right &
+        +openmc.YPlane(y0=0.35) &
+        -openmc.YPlane(y0=0.35 + 2.55) &
+        +openmc.ZPlane(z0=0, boundary_type="reflective") &
+        -openmc.ZPlane(z0=T_pitch * 20, boundary_type="reflective")
+    )
+    return prism
+
+
+def create_prism_vertical(bot, top):
+    yplane_bot = +openmc.YPlane(y0=bot)
+    yplane_top = -openmc.YPlane(y0=top)
+    prism = (
+        +openmc.XPlane(x0=2) &
+        -openmc.XPlane(x0=2 + 23.1) &
+        yplane_bot &
+        yplane_top &
+        +openmc.ZPlane(z0=0, boundary_type="reflective") &
+        -openmc.ZPlane(z0=T_pitch * 20, boundary_type="reflective")
+    )
+    return prism
+
+
+def create_lattice(region, pf):
+    try:
+        centers_1 = openmc.model.pack_spheres(radius=T_r5, region=region, pf=pf)
+        trisos_1 = [openmc.model.TRISO(T_r5, two_triso_univ, c) for c in centers_1]
+        prism = openmc.Cell(region=region)
+        lower_left_1, upper_right_1 = prism.region.bounding_box
+        shape = tuple(((upper_right_1 - lower_left_1) / 0.4).astype(int))
+        pitch_1 = (upper_right_1 - lower_left_1) / shape
+        lattice_1 = openmc.model.create_triso_lattice(
+            trisos_1, lower_left_1, pitch_1, shape, lm_graphite
+        )
+        prism.fill = lattice_1
+    except:
+        prism = openmc.Cell(region=region)
+        prism.fill = lm_graphite
+    return prism

examples/fhr-slab/openmc_input.py

@@ -0,0 +1,133 @@
+import openmc
+import numpy as np
+from numpy import sin, cos, tan, pi
+import sys
+
+sys.path.insert(1, "../")
+from constants import *
+
+# Templating
+total_pf = 0.0979
+sine_a = {{sine_a}}
+sine_b = {{sine_b}}
+sine_c = {{sine_c}}
+vol_total = 23.1 * 2.55 * T_pitch * 20
+vol_slice = 2.31 * 2.55 * T_pitch * 20
+
+x_left = +openmc.XPlane(x0=0, boundary_type="periodic")
+x_right = -openmc.XPlane(x0=27.1, boundary_type="periodic")
+y_top = -openmc.YPlane(y0=3.25, boundary_type="periodic")
+y_bot = +openmc.YPlane(y0=0, boundary_type="periodic")
+y_top.periodic_surface = y_bot
+x_left.periodic_surface = x_right
+z_top = -openmc.ZPlane(z0=T_pitch * 20, boundary_type="reflective")
+z_bot = +openmc.ZPlane(z0=0, boundary_type="reflective")
+bounds = openmc.Cell(fill=flibe)
+bounds.region = x_left & x_right & y_top & y_bot & z_top & z_bot
+
+plank_x_left = +openmc.XPlane(x0=2)
+plank_x_right = -openmc.XPlane(x0=2 + 23.1)
+plank_y_top = -openmc.YPlane(y0=0.35 + 2.55)
+plank_y_bot = +openmc.YPlane(y0=0.35)
+plank_region = plank_x_left & plank_x_right & plank_y_top & plank_y_bot & z_top & z_bot
+bounds.region &= ~plank_region
+
+graphite1_x_right = -openmc.XPlane(x0=2)
+graphite1 = openmc.Cell(fill=graphite)
+graphite1.region = x_left & graphite1_x_right & y_top & y_bot & z_top & z_bot
+bounds.region &= ~graphite1.region
+
+graphite2_x_left = +openmc.XPlane(x0=25.1)
+graphite2 = openmc.Cell(fill=graphite)
+graphite2.region = graphite2_x_left & x_right & y_top & y_bot & z_top & z_bot
+bounds.region &= ~graphite2.region
+
+boundaries = np.arange(2, 27.1, 2.31)
+prism_1 = create_prism(boundaries[0], boundaries[1], False, False)
+prism_2 = create_prism(boundaries[1], boundaries[2], False, False)
+prism_3 = create_prism(boundaries[2], boundaries[3], False, False)
+prism_4 = create_prism(boundaries[3], boundaries[4], False, False)
+prism_5 = create_prism(boundaries[4], boundaries[5], False, False)
+prism_6 = create_prism(boundaries[5], boundaries[6], False, False)
+prism_7 = create_prism(boundaries[6], boundaries[7], False, False)
+prism_8 = create_prism(boundaries[7], boundaries[8], False, False)
+prism_9 = create_prism(boundaries[8], boundaries[9], False, False)
+prism_10 = create_prism(boundaries[9], boundaries[10], False, False)
+
+# triso PF distribution
+vol_triso = 4 / 3 * pi * T_r5 ** 3
+no_trisos = total_pf * vol_total / vol_triso
+midpoints = []
+for x in range(len(boundaries) - 1):
+    midpoints.append((boundaries[x] + boundaries[x + 1]) / 2)
+midpoints = np.array(midpoints)
+sine_val = sine_a * sin(sine_b * midpoints + sine_c) + 2
+sine_val = np.where(sine_val < 0, 0, sine_val)
+triso_z = sine_val / sum(sine_val) * no_trisos
+pf_z = triso_z * vol_triso / vol_slice
+
+prism_cell_1 = create_lattice(prism_1, pf_z[0])
+prism_cell_2 = create_lattice(prism_2, pf_z[1])
+prism_cell_3 = create_lattice(prism_3, pf_z[2])
+prism_cell_4 = create_lattice(prism_4, pf_z[3])
+prism_cell_5 = create_lattice(prism_5, pf_z[4])
+prism_cell_6 = create_lattice(prism_6, pf_z[5])
+prism_cell_7 = create_lattice(prism_7, pf_z[6])
+prism_cell_8 = create_lattice(prism_8, pf_z[7])
+prism_cell_9 = create_lattice(prism_9, pf_z[8])
+prism_cell_10 = create_lattice(prism_10, pf_z[9])
+
+univ = openmc.Universe(
+    cells=[
+        bounds,
+        graphite1,
+        graphite2,
+        prism_cell_1,
+        prism_cell_2,
+        prism_cell_3,
+        prism_cell_4,
+        prism_cell_5,
+        prism_cell_6,
+        prism_cell_7,
+        prism_cell_8,
+        prism_cell_9,
+        prism_cell_10,
+    ]
+)
+geom = openmc.Geometry(univ)
+
+# settings
+point = openmc.stats.Point((13.5, 1.7, T_pitch * 9.5))
+src = openmc.Source(space=point)
+settings = openmc.Settings()
+settings.source = src
+settings.batches = 10
+settings.inactive = 2
+settings.particles = 100
+settings.temperature = {"multipole": True, "method": "interpolation"}
+
+plot = openmc.Plot()
+plot.basis = "xy"
+plot.origin = (13.5, 1.7, T_pitch * 9.5)
+plot.width = (30, 4)
+plot.pixels = (1000, 200)
+colors = {
+    uoc_9: "yellow",
+    por_c: "black",
+    si_c: "orange",
+    graphite: "grey",
+    flibe: "blue",
+    lm_graphite: "red",
+}
+plot.color_by = "material"
+plot.colors = colors
+plots = openmc.Plots()
+plots.append(plot)
+
+# export
+mats.export_to_xml()
+geom.export_to_xml()
+settings.export_to_xml()
+plots.export_to_xml()
+openmc.run()
+# openmc.run(openmc_exec="openmc-ccm-nompi",threads=16)

examples/fhr-slab/rollo_input.json

@@ -0,0 +1,32 @@
+{
+    "control_variables": {
+        "sine_a": {"min": 0.0, "max": 2.0},
+        "sine_b": {"min": 0.0, "max": 1.57},
+        "sine_c": {"min": 0.0, "max": 6.28}
+    },
+    "evaluators": {
+        "openmc": {
+            "input_script": "openmcinp.py",
+            "inputs": ["sine_a", "sine_b", "sine_c"],
+            "outputs": ["keff"],
+            "keep_files": false
+        }
+    },
+    "constraints": {"keff": {"operator": [">="], "constrained_val": [1.0]}},
+    "algorithm": {
+        "parallel": "multiprocessing",
+        "objective": ["max"],
+        "optimized_variable": ["keff"],
+        "pop_size": 4,
+        "generations": 10,
+        "mutation_probability": 0.2374127402121101,
+        "mating_probability": 0.4699131568275016,
+        "selection_operator": {"operator": "selTournament", "inds": 1, "tournsize": 5},
+        "mutation_operator": {
+            "operator": "mutPolynomialBounded",
+            "eta": 0.2374127402121101,
+            "indpb": 0.2374127402121101
+        },
+        "mating_operator": {"operator": "cxBlend", "alpha": 0.4699131568275016}
+    }
+}

pyproject.toml

@@ -0,0 +1,6 @@
+[build-system]
+requires = [
+    "setuptools>=42",
+    "wheel"
+]
+build-backend = "setuptools.build_meta"

docs/lit-review/Makefile → reports/lit-review/Makefile

@@ -24,7 +24,7 @@ epub:
 	ebook-convert $(manuscript).html $(manuscript).epub
 
 clean:
-	rm -f *.bcf *.run.xml *.fdb_latexmk *.fls *.dvi *.toc *.aux *.gz *.out *.log *.bbl *.blg *.log *.spl *~ *.spl *.zip *.acn *.glo *.ist *.epub
+	rm -f *.bcf *.run.xml *.fdb_latexmk *.fls *.dvi *.toc *.aux *.gz *.out *.log *.bbl *.blg *.log *.spl *~ *.spl *.zip *.acn *.glo *.ist *.epub, *.glsdefs
 
 realclean: clean
 	rm -rf $(manuscript).dvi

rollo/__init__.py

@@ -0,0 +1,12 @@
+from rollo.algorithm import *
+from rollo.backend import *
+from rollo.constraints import *
+from rollo.evaluation import *
+from rollo.executor import *
+from rollo.input_validation import *
+from rollo.moltres_evaluation import *
+from rollo.openmc_evaluation import *
+from rollo.special_variables import *
+from rollo.toolbox_generator import *
+
+__version__ = "0.1.1"

rollo/__main__.py

@@ -0,0 +1,29 @@
+from rollo import executor
+import sys
+import getopt
+
+
+def main():
+    argv = sys.argv[1:]
+    msg = "python rollo -i <inputfile> -c <checkpoint file> "
+    try:
+        opts, args = getopt.getopt(argv, "i:c:")
+        opts_dict = {}
+        for opt, arg in opts:
+            opts_dict[opt] = arg
+        if "-i" in opts_dict:
+            if "-c" in opts_dict:
+                new_run = executor.Executor(
+                    input_file=opts_dict["-i"], checkpoint_file=opts_dict["-c"]
+                )
+            else:
+                new_run = executor.Executor(input_file=opts_dict["-i"])
+            new_run.execute()
+        if len(opts) == 0:
+            raise Exception("To run rollo: " + msg)
+    except getopt.GetoptError:
+        raise Exception("To run rollo: " + msg)
+
+
+if __name__ == "__main__":
+    main()