From 07ef67fee019b3c2eb580fc651e6bbbf0bf2e0cb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 26 May 2016 07:46:07 -0500 Subject: [PATCH] Initial TRISO modeling capabilities --- docs/source/pythonapi/index.rst | 27 +++++ openmc/model/__init__.py | 1 + openmc/model/triso.py | 172 ++++++++++++++++++++++++++++++ setup.py | 3 +- tests/test_triso/inputs_true.dat | 1 + tests/test_triso/results_true.dat | 2 + tests/test_triso/test_triso.py | 107 +++++++++++++++++++ 7 files changed, 312 insertions(+), 1 deletion(-) create mode 100644 openmc/model/__init__.py create mode 100644 openmc/model/triso.py create mode 100644 tests/test_triso/inputs_true.dat create mode 100644 tests/test_triso/results_true.dat create mode 100644 tests/test_triso/test_triso.py diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index bf35e7587..54cb93594 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -296,6 +296,33 @@ Multi-group Cross Section Libraries openmc.mgxs.Library +------------------------------------- +:mod:`openmc.model` -- Model Building +------------------------------------- + +TRISO Fuel Modeling +------------------- + +Classes ++++++++ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.model.TRISO + +Functions ++++++++++ + +.. autosummary:: + :toctree: generated + :nosignatures: + + openmc.model.create_triso_lattice + + .. _Jupyter: https://jupyter.org/ .. _NumPy: http://www.numpy.org/ .. _Codecademy: https://www.codecademy.com/tracks/python diff --git a/openmc/model/__init__.py b/openmc/model/__init__.py new file mode 100644 index 000000000..ffb1f4282 --- /dev/null +++ b/openmc/model/__init__.py @@ -0,0 +1 @@ +from .triso import * diff --git a/openmc/model/triso.py b/openmc/model/triso.py new file mode 100644 index 000000000..dc53e6ec0 --- /dev/null +++ b/openmc/model/triso.py @@ -0,0 +1,172 @@ +import copy +from collections import Iterable +from numbers import Real + +import numpy as np + +import openmc +import openmc.checkvalue as cv + +class TRISO(object): + """Tristructural-isotopic (TRISO) micro fuel particle + + Parameters + ---------- + materials : Iterable of openmc.Material + Material to be assigned to each layer of the TRISO particle starting + with the innermost and proceeding outwards + radii : Iterable of float + Outer radii in cm of each layer of the TRISO particle in ascending order + center : Iterable of float + Cartesian coordinates of the center of the TRISO particle in cm + + Attributes + ---------- + cells : list of opemc.Cell + Each layer of the TRISO particle + center : numpy.ndarray + Cartesian coordinates of the center of the TRISO particle in cm + outside : openmc.Region + Region of space outside of the TRISO particle + bounding_box : tuple of numpy.ndarray + Lower-left and upper-right coordinates of an axis-aligned bounding box + for the TRISO particle + + """ + + def __init__(self, materials, radii, center=(0., 0., 0.)): + surfaces = [openmc.Sphere(R=r) for r in radii] + cells = [] + for i, m in enumerate(materials): + c = openmc.Cell(fill=m) + if i == 0: + c.region = -surfaces[i] + else: + c.region = +surfaces[i-1] & -surfaces[i] + cells.append(c) + self._cells = cells + self._surfaces = surfaces + self.center = np.asarray(center) + + @property + def bounding_box(self): + return self.cells[-1].region.bounding_box + + @property + def cells(self): + return self._cells + + @property + def center(self): + return self._center + + @property + def outside(self): + return ~self.cells[-1].region.nodes[-1] + + @center.setter + def center(self, center): + cv.check_type('TRISO center', center, Iterable, Real) + for s in self._surfaces: + s.x0, s.y0, s.z0 = center + self._center = center + + def classify(self, lattice): + """Determine lattice element indices which might contain the TRISO particle. + + Parameters + ---------- + lattice : openmc.RectLattice + Lattice to check + + Returns + ------- + list of tuple + (z,y,x) lattice element indices which might contain the TRISO + particle. + + """ + + ll, ur = self.bounding_box + if lattice.ndim == 2: + (i_min, j_min), p = lattice.find_element(ll) + (i_max, j_max), p = lattice.find_element(ur) + return list(np.broadcast(*np.ogrid[ + j_min:j_max+1, i_min:i_max+1])) + else: + (i_min, j_min, k_min), p = lattice.find_element(ll) + (i_max, j_max, k_max), p = lattice.find_element(ur) + return list(np.broadcast(*np.ogrid[ + k_min:k_max+1, j_min:j_max+1, i_min:i_max+1])) + + +def create_triso_lattice(trisos, lower_left, pitch, shape, background): + """Create a lattice containing TRISO particles for optimized tracking. + + Parameters + ---------- + trisos : list of openmc.model.TRISO + List of TRISO particles to put in lattice + lower_left : Iterable of float + Lower-left Cartesian coordinates of the lattice + pitch : Iterable of float + Pitch of the lattice elements in the x-, y-, and z-directions + shape : Iterable of float + Number of lattice elements in the x-, y-, and z-directions + background : openmc.Material + A background material that is used anywhere within the lattice but + outside a TRISO particle + + Returns + ------- + lattice : openmc.RectLattice + A lattice containing the TRISO particles + + """ + + lattice = openmc.RectLattice() + lattice.lower_left = lower_left + lattice.pitch = pitch + + indices = list(np.broadcast(*np.ogrid[:shape[2], :shape[1], :shape[0]])) + triso_locations = {idx: [] for idx in indices} + for t in trisos: + for idx in t.classify(lattice): + if idx in sorted(triso_locations): + # Create copy of TRISO particle with materials preserved and + # different cell/surface IDs + t_copy = copy.deepcopy(t) + for c, c_copy in zip(t.cells, t_copy.cells): + c_copy.id = None + c_copy.fill = c.fill + for s in t_copy._surfaces: + s.id = None + triso_locations[idx].append(t_copy) + + # Create universes + universes = np.empty(shape[::-1], dtype=openmc.Universe) + for idx, triso_list in sorted(triso_locations.items()): + if len(triso_list) > 0: + outside_trisos = openmc.Intersection(*[t.outside for t in triso_list]) + background_cell = openmc.Cell(fill=background, region=outside_trisos) + else: + background_cell = openmc.Cell(fill=background) + + u = openmc.Universe() + u.add_cell(background_cell) + for t in triso_list: + u.add_cells(t.cells) + iz, iy, ix = idx + t.center = lattice.get_local_coordinates(t.center, (ix, iy, iz)) + + if len(shape) == 2: + universes[-1 - idx[0], idx[1]] = u + else: + universes[idx[0], -1 - idx[1], idx[2]] = u + lattice.universes = universes + + # Set outer universe + background_cell = openmc.Cell(fill=background) + lattice.outer = openmc.Universe(cells=[background_cell]) + + return lattice diff --git a/setup.py b/setup.py index 770f280ad..99e465c9d 100644 --- a/setup.py +++ b/setup.py @@ -11,7 +11,8 @@ except ImportError: kwargs = {'name': 'openmc', 'version': '0.7.1', - 'packages': ['openmc', 'openmc.data', 'openmc.mgxs', 'openmc.stats'], + 'packages': ['openmc', 'openmc.data', 'openmc.mgxs', 'openmc.model', + 'openmc.stats'], 'scripts': glob.glob('scripts/openmc-*'), # Metadata diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat new file mode 100644 index 000000000..40f312850 --- /dev/null +++ b/tests/test_triso/inputs_true.dat @@ -0,0 +1 @@ +6c6fecedf3db0b91b69b7b7b57b3a52c42947253bea4ee3889d6bbd6e74935cc4e599c1c743eb589a1045868412f3f88c5fef7cf10e180bdc4bd182970c9ee15 \ No newline at end of file diff --git a/tests/test_triso/results_true.dat b/tests/test_triso/results_true.dat new file mode 100644 index 000000000..8f968f65a --- /dev/null +++ b/tests/test_triso/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.685303E+00 1.121936E-01 diff --git a/tests/test_triso/test_triso.py b/tests/test_triso/test_triso.py new file mode 100644 index 000000000..f7f4daa7b --- /dev/null +++ b/tests/test_triso/test_triso.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python + +import os +import sys +import glob +import random +from math import sqrt + +import numpy as np + +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc +import openmc.model + + +class TRISOTestHarness(PyAPITestHarness): + def _build_inputs(self): + # Define TRISO matrials + fuel = openmc.Material() + fuel.set_density('g/cm3', 10.5) + fuel.add_nuclide('U-235', 0.14154) + fuel.add_nuclide('U-238', 0.85846) + fuel.add_nuclide('C-Nat', 0.5) + fuel.add_nuclide('O-16', 1.5) + + porous_carbon = openmc.Material() + porous_carbon.set_density('g/cm3', 1.0) + porous_carbon.add_nuclide('C-Nat', 1.0) + porous_carbon.add_s_alpha_beta('Graph', '71t') + + ipyc = openmc.Material() + ipyc.set_density('g/cm3', 1.90) + ipyc.add_nuclide('C-Nat', 1.0) + ipyc.add_s_alpha_beta('Graph', '71t') + + sic = openmc.Material() + sic.set_density('g/cm3', 3.20) + sic.add_element('Si', 1.0) + sic.add_nuclide('C-Nat', 1.0) + + opyc = openmc.Material() + opyc.set_density('g/cm3', 1.87) + opyc.add_nuclide('C-Nat', 1.0) + opyc.add_s_alpha_beta('Graph', '71t') + + graphite = openmc.Material() + graphite.set_density('g/cm3', 1.1995) + graphite.add_nuclide('C-Nat', 1.0) + graphite.add_s_alpha_beta('Graph', '71t') + + # Create TRISO particles + materials = [fuel, porous_carbon, ipyc, sic, opyc] + radii = np.array([212.5, 312.5, 347.5, 382.5, 422.5])*1e-4 + trisos = [] + random.seed(1) + for i in range(100): + # Randomly sample location + x = random.uniform(-0.5, 0.5) + y = random.uniform(-0.5, 0.5) + z = random.uniform(-0.5, 0.5) + t = openmc.model.TRISO(materials, radii, (x, y, z)) + + # Make sure TRISO doesn't overlap with another + for tp in trisos: + xp, yp, zp = tp.center + distance = sqrt((x - xp)**2 + (y - yp)**2 + (z - zp)**2) + if distance <= 2*radii[-1]: + break + else: + trisos.append(t) + + # Define box to contain lattice + min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective') + max_x = openmc.XPlane(x0=0.5, boundary_type='reflective') + min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective') + max_y = openmc.YPlane(y0=0.5, boundary_type='reflective') + min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective') + max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective') + box = openmc.Cell(region=+min_x & -max_x & +min_y & -max_y & +min_z & -max_z) + + # Create lattice + ll, ur = box.region.bounding_box + shape = (3, 3, 3) + lattice = openmc.model.create_triso_lattice( + trisos, ll, (ur - ll)/shape, shape, graphite) + box.fill = lattice + + root = openmc.Universe(0, cells=[box]) + geom = openmc.Geometry(root) + geom.export_to_xml() + + settings = openmc.Settings() + settings.batches = 5 + settings.inactive = 0 + settings.particles = 50 + settings.source = openmc.Source(space=openmc.stats.Point()) + settings.export_to_xml() + + mats = openmc.Materials([fuel, porous_carbon, ipyc, sic, opyc, graphite]) + mats.default_xs = '71c' + mats.export_to_xml() + + +if __name__ == '__main__': + harness = TRISOTestHarness('statepoint.5.h5') + harness.main()