From d0b1287b2ef273ee142a6ebb11030ffdbf580bf1 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 13 Sep 2018 20:10:45 -0500 Subject: [PATCH 1/9] Added spherical shell domain --- openmc/model/triso.py | 147 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 139 insertions(+), 8 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 1e4b4aa9b..50ad70680 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -387,7 +387,7 @@ class _CylindricalDomain(_Domain): @property def volume(self): - return self.length * pi * self.radius**2 + return self.length*pi*self.radius**2 @length.setter def length(self, length): @@ -493,7 +493,7 @@ class _SphericalDomain(_Domain): @property def volume(self): - return 4/3 * pi * self.radius**3 + return 4/3*pi*self.radius**3 @radius.setter def radius(self, radius): @@ -535,6 +535,125 @@ class _SphericalDomain(_Domain): q *= r_max/r +class _SphericalShellDomain(_Domain): + """Spherical shell container in which to pack particles. + + Parameters + ---------- + radius : float + Outer radius of the spherical shell container. + inner_radius : float + Inner radius of the spherical shell container. + center : Iterable of float + Cartesian coordinates of the center of the container. Default is + [0., 0., 0.] + + Attributes + ---------- + radius : float + Outer radius of the spherical shell container. + inner_radius : float + Inner radius of the spherical shell container. + particle_radius : float + Radius of particles to be packed in container. + center : list of float + Cartesian coordinates of the center of the container. Default is + [0., 0., 0.] + cell_length : list of float + Length in x-, y-, and z- directions of each cell in mesh overlaid on + domain. + limits : list of float + Maximum radial distance and minimum radial distance where particle + center can be placed. + volume : float + Volume of the container. + + """ + + def __init__(self, radius, inner_radius, particle_radius, + center=[0., 0., 0.]): + super().__init__(particle_radius, center) + self.radius = radius + self.inner_radius = inner_radius + + @property + def radius(self): + return self._radius + + @property + def inner_radius(self): + return self._inner_radius + + @property + def limits(self): + if self._limits is None: + self._limits = [self.radius - self.particle_radius, + self.inner_radius + self.particle_radius] + return self._limits + + @property + def cell_length(self): + if self._cell_length is None: + mesh_length = 3*[2*self.radius] + self._cell_length = [x/int(x/(4*self.particle_radius)) + for x in mesh_length] + return self._cell_length + + @property + def volume(self): + return 4/3*pi*(self.radius**3 - self.inner_radius**3) + + @radius.setter + def radius(self, radius): + self._radius = float(radius) + self._limits = None + self._cell_length = None + + @inner_radius.setter + def inner_radius(self, inner_radius): + self._inner_radius = float(inner_radius) + self._limits = None + + @limits.setter + def limits(self, limits): + self._limits = limits + + def random_point(self): + r_max = self.limits[0] + r_min = self.limits[1] + x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) + r = (uniform(r_min**3, r_max**3)**(1/3)/sqrt(x[0]**2 + x[1]**2 + x[2]**2)) + return [r*s for s in x] + + def repel_particles(self, p, q, d, d_new): + # Moving each particle distance 's' away from the other along the line + # joining the particle centers will ensure their final distance is + # equal to the outer diameter + s = (d_new - d)/2 + + v = (p - q)/d + p += s*v + q -= s*v + + # Enforce the rigid boundary by moving each particle back along the + # surface normal until it is completely within the container if it + # overlaps the surface + r_max = self.limits[0] + r_min = self.limits[1] + + r = sqrt(p[0]**2 + p[1]**2 + p[2]**2) + if r > r_max: + p *= r_max/r + elif r < r_min: + p *= r_min/r + + r = sqrt(q[0]**2 + q[1]**2 + q[2]**2) + if r > r_max: + q *= r_max/r + elif r < r_min: + q *= r_min/r + + def create_triso_lattice(trisos, lower_left, pitch, shape, background): """Create a lattice containing TRISO particles for optimized tracking. @@ -937,9 +1056,10 @@ def _close_random_pack(domain, particles, contraction_rate): def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, - domain_radius=None, domain_center=[0., 0., 0.], - n_particles=None, packing_fraction=None, - initial_packing_fraction=0.3, contraction_rate=1.e-3, seed=1): + domain_radius=None, domain_inner_radius=None, + domain_center=[0., 0., 0.], n_particles=None, + packing_fraction=None, initial_packing_fraction=0.3, + contraction_rate=1.e-3, seed=1): """Generate a random, non-overlapping configuration of TRISO particles within a container. @@ -955,6 +1075,8 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, Length of the container (if cube or cylinder). domain_radius : float Radius of the container (if cylinder or sphere). + domain_inner_radius : float + Inner radius of the container (if spherical shell). domain_center : Iterable of float Cartesian coordinates of the center of the container. n_particles : int @@ -1018,16 +1140,20 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, """ # Check for valid container geometry and dimensions - if domain_shape not in ['cube', 'cylinder', 'sphere']: + if domain_shape not in ['cube', 'cylinder', 'sphere', 'spherical shell']: raise ValueError('Unable to set domain_shape to "{}". Only "cube", ' - '"cylinder", and "sphere" are ' + '"cylinder", "sphere", and "spherical shell" are ' 'supported."'.format(domain_shape)) if not domain_length and domain_shape in ['cube', 'cylinder']: raise ValueError('"domain_length" must be specified for {} domain ' 'geometry '.format(domain_shape)) - if not domain_radius and domain_shape in ['cylinder', 'sphere']: + if not domain_radius and domain_shape in ['cylinder', 'sphere', + 'spherical shell']: raise ValueError('"domain_radius" must be specified for {} domain ' 'geometry '.format(domain_shape)) + if not domain_inner_radius and domain_shape in ['spherical shell']: + raise ValueError('"domain_inner_radius" must be specified for {} domain ' + 'geometry '.format(domain_shape)) if domain_shape == 'cube': domain = _CubicDomain(length=domain_length, particle_radius=radius, @@ -1038,6 +1164,11 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, elif domain_shape == 'sphere': domain = _SphericalDomain(radius=domain_radius, particle_radius=radius, center=domain_center) + elif domain_shape == 'spherical shell': + domain = _SphericalShellDomain(radius=domain_radius, + inner_radius=domain_inner_radius, + particle_radius=radius, + center=domain_center) # Calculate the packing fraction if the number of particles is specified; # otherwise, calculate the number of particles from the packing fraction. From 8d6541951bdeec065531b72bc55080ecc63eb441 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 17 Sep 2018 17:48:05 -0500 Subject: [PATCH 2/9] SphericalShellDomain extends SphericalDomain --- openmc/model/triso.py | 25 +++---------------------- 1 file changed, 3 insertions(+), 22 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 50ad70680..0d77ad670 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -486,7 +486,7 @@ class _SphericalDomain(_Domain): @property def cell_length(self): if self._cell_length is None: - mesh_length = [2*self.radius, 2*self.radius, 2*self.radius] + mesh_length = 3*[2*self.radius] self._cell_length = [x/int(x/(4*self.particle_radius)) for x in mesh_length] return self._cell_length @@ -535,7 +535,7 @@ class _SphericalDomain(_Domain): q *= r_max/r -class _SphericalShellDomain(_Domain): +class _SphericalShellDomain(_SphericalDomain): """Spherical shell container in which to pack particles. Parameters @@ -572,14 +572,9 @@ class _SphericalShellDomain(_Domain): def __init__(self, radius, inner_radius, particle_radius, center=[0., 0., 0.]): - super().__init__(particle_radius, center) - self.radius = radius + super().__init__(radius, particle_radius, center) self.inner_radius = inner_radius - @property - def radius(self): - return self._radius - @property def inner_radius(self): return self._inner_radius @@ -591,24 +586,10 @@ class _SphericalShellDomain(_Domain): self.inner_radius + self.particle_radius] return self._limits - @property - def cell_length(self): - if self._cell_length is None: - mesh_length = 3*[2*self.radius] - self._cell_length = [x/int(x/(4*self.particle_radius)) - for x in mesh_length] - return self._cell_length - @property def volume(self): return 4/3*pi*(self.radius**3 - self.inner_radius**3) - @radius.setter - def radius(self, radius): - self._radius = float(radius) - self._limits = None - self._cell_length = None - @inner_radius.setter def inner_radius(self, inner_radius): self._inner_radius = float(inner_radius) From 1107457a90176886f28abe903e9c368f4cf8d99f Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 25 Oct 2018 23:45:05 -0500 Subject: [PATCH 3/9] Separate functions for packing spheres and creating TRISOs --- openmc/model/triso.py | 47 ++++++++++++++++++++++++++++--------------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 0d77ad670..adbfe15a7 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -1036,20 +1036,18 @@ def _close_random_pack(domain, particles, contraction_rate): break -def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, - domain_radius=None, domain_inner_radius=None, - domain_center=[0., 0., 0.], n_particles=None, - packing_fraction=None, initial_packing_fraction=0.3, - contraction_rate=1.e-3, seed=1): - """Generate a random, non-overlapping configuration of TRISO particles - within a container. +def pack_spheres(radius, domain_shape='cylinder', domain_length=None, + domain_radius=None, domain_inner_radius=None, + domain_center=[0., 0., 0.], n_particles=None, + packing_fraction=None, initial_packing_fraction=0.3, + contraction_rate=1.e-3, seed=1): + """Generate a random, non-overlapping configuration of spheres within a + container. Parameters ---------- radius : float Outer radius of TRISO particles. - fill : openmc.Universe - Universe which contains all layers of the TRISO particle. domain_shape : {'cube', 'cylinder', or 'sphere'} Geometry of the container in which the TRISO particles are packed. domain_length : float @@ -1080,9 +1078,9 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, RNG seed. Returns - ------- - trisos : list of openmc.model.TRISO - List of TRISO particles in the domain. + ------ + numpy.ndarray + Cartesian coordinates of sphere centers. Notes ----- @@ -1202,7 +1200,24 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None, domain.particle_radius = radius _close_random_pack(domain, particles, contraction_rate) - trisos = [] - for p in particles: - trisos.append(TRISO(radius, fill, [x + c for x, c in zip(p, domain.center)])) - return trisos + return particles + domain.center + +def create_trisos(outer_radius, fill, centers): + """Create TRISO particles at the given coordinates. + + Parameters + ---------- + outer_radius : float + Outer radius of TRISO particle + fill : openmc.Universe + Universe which contains all layers of the TRISO particle + centers : Iterable of float + Cartesian coordinates of the centers of the TRISO particles in cm + + Returns + ------- + list of openmc.model.TRISO + List of TRISO particles in the domain. + + """ + return [TRISO(outer_radius, fill, c) for c in centers] From c443208e21a1b245243455e84cc89f5af6946b16 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 8 Nov 2018 13:07:42 -0800 Subject: [PATCH 4/9] Generate sphere packing container from openmc.Region --- openmc/model/triso.py | 851 ++++++++++++++++++++++++++---------------- openmc/region.py | 2 +- openmc/surface.py | 4 +- 3 files changed, 540 insertions(+), 317 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index adbfe15a7..03ca6f5e0 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -2,19 +2,24 @@ import copy import warnings import itertools import random -from collections import defaultdict +from abc import ABCMeta, abstractproperty, abstractmethod +from collections import Counter, defaultdict from collections.abc import Iterable -from numbers import Real -from random import uniform, gauss from heapq import heappush, heappop from math import pi, sin, cos, floor, log10, sqrt -from abc import ABCMeta, abstractproperty, abstractmethod +from numbers import Real +from operator import attrgetter +from random import uniform, gauss import numpy as np import scipy.spatial import openmc -import openmc.checkvalue as cv +from openmc.checkvalue import check_type + + +MAX_PF_RSP = 0.38 +MAX_PF_CRP = 0.64 class TRISO(openmc.Cell): @@ -55,7 +60,7 @@ class TRISO(openmc.Cell): @center.setter def center(self, center): - cv.check_type('TRISO center', center, Iterable, Real) + check_type('TRISO center', center, Iterable, Real) self._surface.x0 = center[0] self._surface.y0 = center[1] self._surface.z0 = center[2] @@ -91,43 +96,43 @@ class TRISO(openmc.Cell): k_min:k_max+1, j_min:j_max+1, i_min:i_max+1])) -class _Domain(metaclass=ABCMeta): - """Container in which to pack particles. +class _Container(metaclass=ABCMeta): + """Container in which to pack spheres. Parameters ---------- - particle_radius : float - Radius of particles to be packed in container. + sphere_radius : float + Radius of spheres to be packed in container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- - particle_radius : float - Radius of particles to be packed in container. + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Constraint on where particle center can be placed. + Constraint on where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, particle_radius, center=[0., 0., 0.]): + def __init__(self, sphere_radius, center=(0., 0., 0.)): self._cell_length = None self._limits = None - self.particle_radius = particle_radius + self.sphere_radius = sphere_radius self.center = center @property - def particle_radius(self): - return self._particle_radius + def sphere_radius(self): + return self._sphere_radius @property def center(self): @@ -145,9 +150,9 @@ class _Domain(metaclass=ABCMeta): def volume(self): pass - @particle_radius.setter - def particle_radius(self, particle_radius): - self._particle_radius = float(particle_radius) + @sphere_radius.setter + def sphere_radius(self, sphere_radius): + self._sphere_radius = float(sphere_radius) self._limits = None self._cell_length = None @@ -160,12 +165,12 @@ class _Domain(metaclass=ABCMeta): def mesh_cell(self, p): """Calculate the index of the cell in a mesh overlaid on the domain in - which the given particle center falls. + which the given sphere center falls. Parameters ---------- p : Iterable of float - Cartesian coordinates of particle center. + Cartesian coordinates of sphere center. Returns ------- @@ -177,12 +182,12 @@ class _Domain(metaclass=ABCMeta): def nearby_mesh_cells(self, p): """Calculates the indices of all cells in a mesh overlaid on the domain - within one diameter of the given particle. + within one diameter of the given sphere. Parameters ---------- p : Iterable of float - Cartesian coordinates of particle center. + Cartesian coordinates of sphere center. Returns ------- @@ -190,27 +195,27 @@ class _Domain(metaclass=ABCMeta): Indices of mesh cells. """ - d = 2*self.particle_radius + d = 2*self.sphere_radius r = [[a/self.cell_length[i] for a in [p[i]-d, p[i], p[i]+d]] for i in range(3)] return list(itertools.product(*({int(x) for x in y} for y in r))) @abstractmethod def random_point(self): - """Generate Cartesian coordinates of center of a particle that is + """Generate Cartesian coordinates of center of a sphere that is contained entirely within the domain with uniform probability. Returns ------- list of float - Cartesian coordinates of particle center. + Cartesian coordinates of sphere center. """ pass @abstractmethod - def repel_particles(self, p, q, d, d_new): - """Move particles p and q apart according to the following + def repel_spheres(self, p, q, d, d_new): + """Move spheres p and q apart according to the following transformation (accounting for boundary conditions on domain): r_i^(n+1) = r_i^(n) + 1/2(d_out^(n+1) - d^(n)) @@ -219,78 +224,110 @@ class _Domain(metaclass=ABCMeta): Parameters ---------- p, q : numpy.ndarray - Cartesian coordinates of particle center. + Cartesian coordinates of sphere center. d : float - distance between centers of particles i and j. + distance between centers of spheres i and j. d_new : float - final distance between centers of particles i and j. + final distance between centers of spheres i and j. """ pass -class _CubicDomain(_Domain): - """Cubic container in which to pack particles. +class _RectangularPrism(_Container): + """Rectangular prism container in which to pack spheres. Parameters ---------- - length : float - Length of each side of the cubic container. - particle_radius : float - Radius of particles to be packed in container. + width : float + Prism length along the x-axis + depth : float + Prism length along the y-axis + height : float + Prism length along the z-axis + sphere_radius : float + Radius of spheres to be packed in container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- - length : float - Length of each side of the cubic container. - particle_radius : float - Radius of particles to be packed in container. + width : float + Prism length along the x-axis + depth : float + Prism length along the y-axis + height : float + Prism length along the z-axis + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum distance from center in x-, y-, or z-direction where particle + Maximum distance from center in x-, y-, or z-direction where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, length, particle_radius, center=[0., 0., 0.]): - super().__init__(particle_radius, center) - self.length = length + def __init__(self, width, depth, height, sphere_radius, center=(0., 0., 0.)): + super().__init__(sphere_radius, center) + self.width = width + self.depth = depth + self.height = height @property - def length(self): - return self._length + def width(self): + return self._width + + @property + def depth(self): + return self._depth + + @property + def height(self): + return self._height @property def limits(self): if self._limits is None: - self._limits = [self.length/2 - self.particle_radius] + self._limits = [self.width/2 - self.sphere_radius, + self.depth/2 - self.sphere_radius, + self.height/2 - self.sphere_radius] return self._limits @property def cell_length(self): if self._cell_length is None: - mesh_length = [self.length, self.length, self.length] - self._cell_length = [x/int(x/(4*self.particle_radius)) + mesh_length = [self.width, self.depth, self.height] + self._cell_length = [x/int(x/(4*self.sphere_radius)) for x in mesh_length] return self._cell_length @property def volume(self): - return self.length**3 + return self.width*self.depth*self.height - @length.setter - def length(self, length): - self._length = float(length) + @width.setter + def width(self, width): + self._width = float(width) + self._limits = None + self._cell_length = None + + @depth.setter + def depth(self, depth): + self._depth = float(depth) + self._limits = None + self._cell_length = None + + @height.setter + def height(self, height): + self._height = float(height) self._limits = None self._cell_length = None @@ -299,14 +336,14 @@ class _CubicDomain(_Domain): self._limits = limits def random_point(self): - x_max = self.limits[0] + x_max, y_max, z_max = self.limits return [uniform(-x_max, x_max), - uniform(-x_max, x_max), - uniform(-x_max, x_max)] + uniform(-y_max, y_max), + uniform(-z_max, z_max)] - def repel_particles(self, p, q, d, d_new): - # Moving each particle distance 's' away from the other along the line - # joining the particle centers will ensure their final distance is + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is # equal to the outer diameter s = (d_new - d)/2 @@ -314,53 +351,64 @@ class _CubicDomain(_Domain): p += s*v q -= s*v - # Enforce the rigid boundary by moving each particle back along the + # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - x_max = self.limits[0] - p[:] = np.clip(p, -x_max, x_max) - q[:] = np.clip(q, -x_max, x_max) + p[:] = np.clip(p, [-x for x in self.limits], self.limits) + q[:] = np.clip(q, [-x for x in self.limits], self.limits) -class _CylindricalDomain(_Domain): - """Cylindrical container in which to pack particles. +class _Cylinder(_Container): + """Cylindrical container in which to pack spheres. Parameters ---------- length : float - Length along z-axis of the cylindrical container. + Length of the cylindrical container. radius : float Radius of the cylindrical container. + axis : string + Axis along which the length of the cylinder is aligned. + sphere_radius : float + Radius of spheres to be packed in container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- length : float - Length along z-axis of the cylindrical container. + Length of the cylindrical container. radius : float Radius of the cylindrical container. - particle_radius : float - Radius of particles to be packed in container. + axis : string + Axis along which the length of the cylinder is aligned. + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) + shift : list of int + Rolled indices of the x-, y-, and z- coordinates of a sphere so the + configuration is aligned with the correct axis. No shift corresponds to + a cylinder along the z-axis. cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float Maximum radial distance and maximum distance from center in z-direction - where particle center can be placed. + where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, length, radius, particle_radius, center=[0., 0., 0.]): - super().__init__(particle_radius, center) + def __init__(self, length, radius, axis, sphere_radius, center=(0., 0., 0.)): + super().__init__(sphere_radius, center) + self._shift = None self.length = length self.radius = radius + self.axis = axis @property def length(self): @@ -370,19 +418,37 @@ class _CylindricalDomain(_Domain): def radius(self): return self._radius + @property + def axis(self): + return self._axis + + @property + def shift(self): + if self._shift is None: + if self.axis == 'x': + self._shift = [1, 2, 0] + elif self.axis == 'y': + self._shift = [2, 0, 1] + else: + self._shift = [0, 1, 2] + return self._shift + @property def limits(self): if self._limits is None: - self._limits = [self.radius - self.particle_radius, - self.length/2 - self.particle_radius] + self._limits = [self.radius - self.sphere_radius, + self.length/2 - self.sphere_radius] return self._limits @property def cell_length(self): if self._cell_length is None: - mesh_length = [2*self.radius, 2*self.radius, self.length] - self._cell_length = [x/int(x/(4*self.particle_radius)) - for x in mesh_length] + h = 4*self.sphere_radius + i, j, k = self.shift + self._cell_length = [None]*3 + self._cell_length[i] = 2*self.radius/int(2*self.radius/h) + self._cell_length[j] = 2*self.radius/int(2*self.radius/h) + self._cell_length[k] = self.length/int(self.length/h) return self._cell_length @property @@ -401,20 +467,29 @@ class _CylindricalDomain(_Domain): self._limits = None self._cell_length = None + @axis.setter + def axis(self, axis): + self._axis = axis + self._shift = None + @limits.setter def limits(self, limits): self._limits = limits def random_point(self): - r_max = self.limits[0] - z_max = self.limits[1] + r_max, z_max = self.limits r = sqrt(uniform(0, r_max**2)) t = uniform(0, 2*pi) - return [r*cos(t), r*sin(t), uniform(-z_max, z_max)] + i, j, k = self.shift + p = [None]*3 + p[i] = r*cos(t) + p[j] = r*sin(t) + p[k] = uniform(-z_max, z_max) + return p - def repel_particles(self, p, q, d, d_new): - # Moving each particle distance 's' away from the other along the line - # joining the particle centers will ensure their final distance is + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is # equal to the outer diameter s = (d_new - d)/2 @@ -422,25 +497,27 @@ class _CylindricalDomain(_Domain): p += s*v q -= s*v - # Enforce the rigid boundary by moving each particle back along the + # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - r_max = self.limits[0] - z_max = self.limits[1] + r_max, z_max = self.limits + i, j, k = self.shift - r = sqrt(p[0]**2 + p[1]**2) + r = sqrt(p[i]**2 + p[j]**2) if r > r_max: - p[0:2] *= r_max/r - p[2] = np.clip(p[2], -z_max, z_max) + p[i] *= r_max/r + p[j] *= r_max/r + p[k] = np.clip(p[k], -z_max, z_max) - r = sqrt(q[0]**2 + q[1]**2) + r = sqrt(q[i]**2 + q[j]**2) if r > r_max: - q[0:2] *= r_max/r - q[2] = np.clip(q[2], -z_max, z_max) + q[i] *= r_max/r + q[j] *= r_max/r + q[k] = np.clip(q[k], -z_max, z_max) -class _SphericalDomain(_Domain): - """Spherical container in which to pack particles. +class _Sphere(_Container): + """Spherical container in which to pack spheres. Parameters ---------- @@ -448,29 +525,29 @@ class _SphericalDomain(_Domain): Radius of the spherical container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- radius : float Radius of the spherical container. - particle_radius : float - Radius of particles to be packed in container. + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum radial distance where particle center can be placed. + Maximum radial distance where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, radius, particle_radius, center=[0., 0., 0.]): - super().__init__(particle_radius, center) + def __init__(self, radius, sphere_radius, center=(0., 0., 0.)): + super().__init__(sphere_radius, center) self.radius = radius @property @@ -480,14 +557,14 @@ class _SphericalDomain(_Domain): @property def limits(self): if self._limits is None: - self._limits = [self.radius - self.particle_radius] + self._limits = [self.radius - self.sphere_radius] return self._limits @property def cell_length(self): if self._cell_length is None: mesh_length = 3*[2*self.radius] - self._cell_length = [x/int(x/(4*self.particle_radius)) + self._cell_length = [x/int(x/(4*self.sphere_radius)) for x in mesh_length] return self._cell_length @@ -511,9 +588,9 @@ class _SphericalDomain(_Domain): r = (uniform(0, r_max**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2)) return [r*s for s in x] - def repel_particles(self, p, q, d, d_new): - # Moving each particle distance 's' away from the other along the line - # joining the particle centers will ensure their final distance is + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is # equal to the outer diameter s = (d_new - d)/2 @@ -521,7 +598,7 @@ class _SphericalDomain(_Domain): p += s*v q -= s*v - # Enforce the rigid boundary by moving each particle back along the + # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface r_max = self.limits[0] @@ -535,8 +612,8 @@ class _SphericalDomain(_Domain): q *= r_max/r -class _SphericalShellDomain(_SphericalDomain): - """Spherical shell container in which to pack particles. +class _SphericalShell(_Sphere): + """Spherical shell container in which to pack spheres. Parameters ---------- @@ -546,7 +623,7 @@ class _SphericalShellDomain(_SphericalDomain): Inner radius of the spherical shell container. center : Iterable of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) Attributes ---------- @@ -554,25 +631,25 @@ class _SphericalShellDomain(_SphericalDomain): Outer radius of the spherical shell container. inner_radius : float Inner radius of the spherical shell container. - particle_radius : float - Radius of particles to be packed in container. + sphere_radius : float + Radius of spheres to be packed in container. center : list of float Cartesian coordinates of the center of the container. Default is - [0., 0., 0.] + (0., 0., 0.) cell_length : list of float Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum radial distance and minimum radial distance where particle + Maximum radial distance and minimum radial distance where sphere center can be placed. volume : float Volume of the container. """ - def __init__(self, radius, inner_radius, particle_radius, - center=[0., 0., 0.]): - super().__init__(radius, particle_radius, center) + def __init__(self, radius, inner_radius, sphere_radius, + center=(0., 0., 0.)): + super().__init__(radius, sphere_radius, center) self.inner_radius = inner_radius @property @@ -582,8 +659,8 @@ class _SphericalShellDomain(_SphericalDomain): @property def limits(self): if self._limits is None: - self._limits = [self.radius - self.particle_radius, - self.inner_radius + self.particle_radius] + self._limits = [self.radius - self.sphere_radius, + self.inner_radius + self.sphere_radius] return self._limits @property @@ -600,15 +677,14 @@ class _SphericalShellDomain(_SphericalDomain): self._limits = limits def random_point(self): - r_max = self.limits[0] - r_min = self.limits[1] + r_max, r_min = self.limits x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) r = (uniform(r_min**3, r_max**3)**(1/3)/sqrt(x[0]**2 + x[1]**2 + x[2]**2)) return [r*s for s in x] - def repel_particles(self, p, q, d, d_new): - # Moving each particle distance 's' away from the other along the line - # joining the particle centers will ensure their final distance is + def repel_spheres(self, p, q, d, d_new): + # Moving each sphere distance 's' away from the other along the line + # joining the sphere centers will ensure their final distance is # equal to the outer diameter s = (d_new - d)/2 @@ -616,11 +692,10 @@ class _SphericalShellDomain(_SphericalDomain): p += s*v q -= s*v - # Enforce the rigid boundary by moving each particle back along the + # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - r_max = self.limits[0] - r_min = self.limits[1] + r_max, r_min = self.limits r = sqrt(p[0]**2 + p[1]**2 + p[2]**2) if r > r_max: @@ -708,28 +783,215 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): return lattice -def _random_sequential_pack(domain, n_particles): - """Random sequential packing of particles within a container. +def _create_container(region, sphere_radius): + """Create a container to pack spheres in using the region to determine the + container shape. Parameters ---------- - domain : openmc.model._Domain - Container in which to pack particles. - n_particles : int - Number of particles to pack. + region : openmc.Region + Container in which the spheres are packed. Supported shapes are + rectangular_prism, cylinder, sphere, and spherical shell. + sphere_radius : float + Outer radius of spheres. + + Returns + ------- + domain : openmc.model._Container + Container in which to pack spheres. + + """ + + def rectangular_prism(): + # Assume the simplest case where the prism volume is the intersection + # of the half-spaces of six planes + if not isinstance(region, openmc.Intersection): + return None + + if any(not isinstance(node, openmc.Halfspace) for node in region): + return None + + if len(region) != 6: + return None + + # Sort half-spaces by surface type + px1, px2, py1, py2, pz1, pz2 = sorted(region, key=attrgetter('surface.type')) + + # Make sure the region consists of the correct surfaces + if (not isinstance(px1.surface, openmc.XPlane) or + not isinstance(px2.surface, openmc.XPlane) or + not isinstance(py1.surface, openmc.YPlane) or + not isinstance(py2.surface, openmc.YPlane) or + not isinstance(pz1.surface, openmc.ZPlane) or + not isinstance(pz2.surface, openmc.ZPlane)): + return None + + # Secondary sorting by location of the plane + if px1.surface.x0 > px2.surface.x0: + px1, px2 = px2, px1 + + if py1.surface.y0 > py2.surface.y0: + py1, py2 = py2, py1 + + if pz1.surface.z0 > pz2.surface.z0: + pz1, pz2 = pz2, pz1 + + # Make sure the half-spaces are on the correct side of the surfaces + if (px1.side != '+' or px2.side != '-' or + py1.side != '+' or py2.side != '-' or + pz1.side != '+' or pz2.side != '-'): + return None + + # Calculate the parameters for the container + width = px2.surface.x0 - px1.surface.x0 + depth = py2.surface.y0 - py1.surface.y0 + height = pz2.surface.z0 - pz1.surface.z0 + center = (px1.surface.x0 + width/2, + py1.surface.y0 + depth/2, + pz1.surface.z0 + height/2) + + # The region is the volume of a rectangular prism, so create container + return _RectangularPrism(width, depth, height, sphere_radius, center) + + def cylinder(): + # Assume the simplest case where the cylinder volume is the + # intersection of the half-spaces of a cylinder and two planes + if not isinstance(region, openmc.Intersection): + return None + + if any(not isinstance(node, openmc.Halfspace) for node in region): + return None + + if len(region) != 3: + return None + + # Identify the axis that the cylinder lies along + axis = region[0].surface.type[0] + + # Make sure the region is composed of a cylinder and two planes on the + # same axis + count = Counter((node.surface.type for node in region)) + if count[axis + '-cylinder'] != 1 or count[axis + '-plane'] != 2: + return None + + # Sort the half-spaces by surface type + cyl, p1, p2 = sorted(region, key=attrgetter('surface.type')) + + # Calculate the parameters for a cylinder along the x-axis + if axis == 'x': + if p1.surface.x0 > p2.surface.x0: + p1, p2 = p2, p1 + length = p2.surface.x0 - p1.surface.x0 + center = (length/2, cyl.surface.y0, cyl.surface.z0) + + # Calculate the parameters for a cylinder along the y-axis + elif axis == 'y': + if p1.surface.y0 > p2.surface.y0: + p1, p2 = p2, p1 + length = p2.surface.y0 - p1.surface.y0 + center = (cyl.surface.x0, length/2, cyl.surface.z0) + + # Calculate the parameters for a cylinder along the z-axis + else: + if p1.surface.z0 > p2.surface.z0: + p1, p2 = p2, p1 + length = p2.surface.z0 - p1.surface.z0 + center = (cyl.surface.x0, cyl.surface.y0, length/2) + + # Make sure the half-spaces are on the correct side of the surfaces + if cyl.side != '-' or p1.side != '+' or p2.side != '-': + return None + + radius = cyl.surface.r + + # The region is the volume of a cylinder, so create container + return _Cylinder(length, radius, axis, sphere_radius, center) + + def sphere(): + # Assume the simplest case where the sphere volume is the negative + # half-space of a sphere + if not isinstance(region, openmc.Halfspace): + return None + + if not isinstance(region.surface, openmc.Sphere): + return None + + if region.side != '-': + return None + + radius = region.surface.r + center = (region.surface.x0, region.surface.y0, region.surface.z0) + + # The region is the volume of a sphere, so create container + return _Sphere(radius, sphere_radius, center) + + def spherical_shell(): + # Assume the simplest case where the spherical shell volume is the + # intersection of the half-spaces of two spheres + if not isinstance(region, openmc.Intersection): + return None + + if any(not isinstance(node, openmc.Halfspace) for node in region): + return None + + if len(region) != 2: + return None + + if any(not isinstance(node.surface, openmc.Sphere) for node in region): + return None + + s1, s2 = sorted(region, key=attrgetter('surface.r')) + radius = s2.surface.r + inner_radius = s1.surface.r + center = (s1.surface.x0, s1.surface.y0, s1.surface.z0) + + if center != (s2.surface.x0, s2.surface.y0, s2.surface.z0): + return None + + if s1.side != '+' or s2.side != '-': + return None + + # The region is the volume of a spherical shell, so create container + return _SphericalShell(radius, inner_radius, sphere_radius, center) + + check_type('region', region, openmc.Region) + + # Check whether the region matches any of the supported container shapes + # and create the container if it does + shapes = [rectangular_prism, cylinder, sphere, spherical_shell] + for shape in shapes: + container = shape() + if container: + return container + + msg = ('Could not translate region {} into a container: supported ' + 'container shapes are rectangular prism, cylinder, sphere, ' + 'and spherical shell.'.format(region)) + raise ValueError(msg) + + +def _random_sequential_pack(domain, num_spheres): + """Random sequential packing of spheres within a container. + + Parameters + ---------- + domain : openmc.model._Container + Container in which to pack spheres. + num_spheres : int + Number of spheres to pack. Returns ------ numpy.ndarray - Cartesian coordinates of centers of particles. + Cartesian coordinates of centers of spheres. """ - sqd = (2*domain.particle_radius)**2 - particles = [] + sqd = (2*domain.sphere_radius)**2 + spheres = [] mesh = defaultdict(list) - for i in range(n_particles): + for i in range(num_spheres): # Randomly sample new center coordinates while there are any overlaps while True: p = domain.random_point() @@ -739,23 +1001,23 @@ def _random_sequential_pack(domain, n_particles): continue else: break - particles.append(p) + spheres.append(p) for idx in domain.nearby_mesh_cells(p): mesh[idx].append(p) - return np.array(particles) + return np.array(spheres) -def _close_random_pack(domain, particles, contraction_rate): - """Close random packing of particles using the Jodrey-Tory algorithm. +def _close_random_pack(domain, spheres, contraction_rate): + """Close random packing of spheres using the Jodrey-Tory algorithm. Parameters ---------- - domain : openmc.model._Domain - Container in which to pack particles. - particles : numpy.ndarray - Initial Cartesian coordinates of centers of particles. + domain : openmc.model._Container + Container in which to pack spheres. + spheres : numpy.ndarray + Initial Cartesian coordinates of centers of spheres. contraction_rate : float Contraction rate of outer diameter. @@ -767,9 +1029,9 @@ def _close_random_pack(domain, particles, contraction_rate): Parameters ---------- d : float - distance between centers of particles i and j. + distance between centers of spheres i and j. i, j : int - Index of particles in particles array. + Index of spheres in spheres array. """ @@ -779,12 +1041,12 @@ def _close_random_pack(domain, particles, contraction_rate): heappush(rods, rod) def remove_rod(i): - """Mark the rod containing particle i as removed. + """Mark the rod containing sphere i as removed. Parameters ---------- i : int - Index of particle in particles array. + Index of sphere in spheres array. """ @@ -800,9 +1062,9 @@ def _close_random_pack(domain, particles, contraction_rate): Returns ------- d : float - distance between centers of particles i and j. + distance between centers of spheres i and j. i, j : int - Index of particles in particles array. + Index of spheres in spheres array. """ @@ -815,15 +1077,15 @@ def _close_random_pack(domain, particles, contraction_rate): return None, None, None def create_rod_list(): - """Generate sorted list of rods (distances between particle centers). + """Generate sorted list of rods (distances between sphere centers). Rods are arranged in a heap where each element contains the rod length - and the particle indices. A rod between particles p and q is only + and the sphere indices. A rod between spheres p and q is only included if the distance between p and q could not be changed by the elimination of a greater overlap, i.e. q has no nearer neighbors than p. - A mapping of particle ids to rods is maintained in 'rods_map'. Each key - in the dict is the id of a particle that is in the rod list, and the + A mapping of sphere ids to rods is maintained in 'rods_map'. Each key + in the dict is the id of a sphere that is in the rod list, and the value is the id of its nearest neighbor and the rod that contains them. The dict is used to find rods in the priority queue and to mark removed rods so rods can be "removed" without breaking the heap structure @@ -832,25 +1094,25 @@ def _close_random_pack(domain, particles, contraction_rate): """ # Create KD tree for quick nearest neighbor search - tree = scipy.spatial.cKDTree(particles) + tree = scipy.spatial.cKDTree(spheres) # Find distance to nearest neighbor and index of nearest neighbor for - # all particles - d, n = tree.query(particles, k=2) + # all spheres + d, n = tree.query(spheres, k=2) d = d[:, 1] n = n[:, 1] - # Array of particle indices, indices of nearest neighbors, and + # Array of sphere indices, indices of nearest neighbors, and # distances to nearest neighbors a = np.vstack((list(range(n.size)), n, d)).T # Sort along second column and swap first and second columns to create - # array of nearest neighbor indices, indices of particles they are + # array of nearest neighbor indices, indices of spheres they are # nearest neighbors of, and distances between them b = a[a[:, 1].argsort()] b[:, [0, 1]] = b[:, [1, 0]] - # Find the intersection between 'a' and 'b': a list of particles who + # Find the intersection between 'a' and 'b': a list of spheres who # are each other's nearest neighbors and the distance between them r = list({tuple(x) for x in a} & {tuple(x) for x in b}) @@ -866,30 +1128,30 @@ def _close_random_pack(domain, particles, contraction_rate): add_rod(d, i, j) def update_mesh(i): - """Update which mesh cells the particle is in based on new particle + """Update which mesh cells the sphere is in based on new sphere center coordinates. 'mesh'/'mesh_map' is a two way dictionary used to look up which - particles are located within one diameter of a given mesh cell and - which mesh cells a given particle center is within one diameter of. + spheres are located within one diameter of a given mesh cell and + which mesh cells a given sphere center is within one diameter of. This is used to speed up the nearest neighbor search. Parameters ---------- i : int - Index of particle in particles array. + Index of sphere in spheres array. """ - # Determine which mesh cells the particle is in and remove the - # particle id from those cells + # Determine which mesh cells the sphere is in and remove the + # sphere id from those cells for idx in mesh_map[i]: mesh[idx].remove(i) del mesh_map[i] - # Determine which mesh cells are within one diameter of particle's - # center and add this particle to the list of particles in those cells - for idx in domain.nearby_mesh_cells(particles[i]): + # Determine which mesh cells are within one diameter of sphere's + # center and add this sphere to the list of spheres in those cells + for idx in domain.nearby_mesh_cells(spheres[i]): mesh[idx].add(i) mesh_map[i].add(idx) @@ -898,7 +1160,7 @@ def _close_random_pack(domain, particles, contraction_rate): d_out^(i+1) = d_out^(i) - (1/2)^(j) * d_out0 * k / n, - where k is the contraction rate, n is the number of particles, and + where k is the contraction rate, n is the number of spheres, and j = floor(-log10(pf_out - pf_in)). @@ -909,26 +1171,26 @@ def _close_random_pack(domain, particles, contraction_rate): """ - inner_pf = 4/3*pi*(inner_diameter/2)**3*n_particles/domain.volume - outer_pf = 4/3*pi*(outer_diameter/2)**3*n_particles/domain.volume + inner_pf = 4/3*pi*(inner_diameter/2)**3*num_spheres/domain.volume + outer_pf = 4/3*pi*(outer_diameter/2)**3*num_spheres/domain.volume j = floor(-log10(outer_pf - inner_pf)) return (outer_diameter - 0.5**j * contraction_rate * - initial_outer_diameter / n_particles) + initial_outer_diameter / num_spheres) def nearest(i): - """Find index of nearest neighbor of particle i. + """Find index of nearest neighbor of sphere i. Parameters ---------- i : int - Index in particles array of particle for which to find nearest + Index in spheres array of sphere for which to find nearest neighbor. Returns ------- int - Index in particles array of nearest neighbor of i + Index in spheres array of nearest neighbor of i float distance between i and nearest neighbor. @@ -937,8 +1199,8 @@ def _close_random_pack(domain, particles, contraction_rate): # Need the second nearest neighbor of i since the nearest neighbor # will be itself. Using argpartition, the k-th nearest neighbor is # placed at index k. - idx = list(mesh[domain.mesh_cell(particles[i])]) - dists = scipy.spatial.distance.cdist([particles[i]], particles[idx])[0] + idx = list(mesh[domain.mesh_cell(spheres[i])]) + dists = scipy.spatial.distance.cdist([spheres[i]], spheres[idx])[0] if dists.size > 1: j = dists.argpartition(1)[1] return idx[j], dists[j] @@ -946,17 +1208,17 @@ def _close_random_pack(domain, particles, contraction_rate): return None, None def update_rod_list(i): - """Update the rod list with the new nearest neighbors of particle since + """Update the rod list with the new nearest neighbors of sphere since its overlap was eliminated. Parameters ---------- i : int - Index of particle in particles array. + Index of sphere in spheres array. """ - # If the nearest neighbor k of particle i has no nearer neighbors, + # If the nearest neighbor k of sphere i has no nearer neighbors, # remove the rod currently containing k from the rod list and add rod # k-i, keeping the rod list sorted k, d_ik = nearest(i) @@ -965,57 +1227,57 @@ def _close_random_pack(domain, particles, contraction_rate): remove_rod(k) add_rod(d_ik, i, k) - n_particles = len(particles) - diameter = 2*domain.particle_radius + num_spheres = len(spheres) + diameter = 2*domain.sphere_radius # Flag for marking rods that have been removed from priority queue removed = -1 # Outer diameter initially set to arbitrary value that yields pf of 1 - initial_outer_diameter = 2*(domain.volume/(n_particles*4/3*pi))**(1/3) + initial_outer_diameter = 2*(domain.volume/(num_spheres*4/3*pi))**(1/3) - # Inner and outer diameter of particles will change during packing + # Inner and outer diameter of spheres will change during packing outer_diameter = initial_outer_diameter inner_diameter = 0. - # List of rods arranged in a heap and mapping of particle ids to rods + # List of rods arranged in a heap and mapping of sphere ids to rods rods = [] rods_map = {} - # Initialize two-way dictionary that identifies which particles are near a - # given mesh cell and which mesh cells a particle is near + # Initialize two-way dictionary that identifies which spheres are near a + # given mesh cell and which mesh cells a sphere is near mesh = defaultdict(set) mesh_map = defaultdict(set) - for i in range(n_particles): - for idx in domain.nearby_mesh_cells(particles[i]): + for i in range(num_spheres): + for idx in domain.nearby_mesh_cells(spheres[i]): mesh[idx].add(i) mesh_map[i].add(idx) while True: - # Rebuild the sorted list of rods according to the current particle + # Rebuild the sorted list of rods according to the current sphere # configuration create_rod_list() # Set the inner diameter to the shortest center-to-center distance - # between any two particles + # between any two spheres if rods: inner_diameter = rods[0][0] - # Reached the desired particle radius + # Reached the desired sphere radius if inner_diameter >= diameter: break - # The algorithm converged before reaching the desired particle radius. + # The algorithm converged before reaching the desired sphere radius. # This can happen when the desired packing fraction is close to the # packing fraction limit. The packing fraction is a random variable - # that is determined by the particle locations and the contraction + # that is determined by the sphere locations and the contraction # rate. A higher packing fraction can be achieved with a smaller # contraction rate, though at the cost of a longer simulation time -- # the number of iterations needed to remove all overlaps is inversely # proportional to the contraction rate. if inner_diameter >= outer_diameter or not rods: warnings.warn('Close random pack converged before reaching true ' - 'particle radius; some particles may overlap. Try ' + 'sphere radius; some spheres may overlap. Try ' 'reducing contraction rate or packing fraction.') break @@ -1024,7 +1286,7 @@ def _close_random_pack(domain, particles, contraction_rate): if not d: break outer_diameter = reduce_outer_diameter() - domain.repel_particles(particles[i], particles[j], d, outer_diameter) + domain.repel_spheres(spheres[i], spheres[j], d, outer_diameter) update_mesh(i) update_mesh(j) update_rod_list(i) @@ -1036,10 +1298,7 @@ def _close_random_pack(domain, particles, contraction_rate): break -def pack_spheres(radius, domain_shape='cylinder', domain_length=None, - domain_radius=None, domain_inner_radius=None, - domain_center=[0., 0., 0.], n_particles=None, - packing_fraction=None, initial_packing_fraction=0.3, +def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, contraction_rate=1.e-3, seed=1): """Generate a random, non-overlapping configuration of spheres within a container. @@ -1047,33 +1306,27 @@ def pack_spheres(radius, domain_shape='cylinder', domain_length=None, Parameters ---------- radius : float - Outer radius of TRISO particles. - domain_shape : {'cube', 'cylinder', or 'sphere'} - Geometry of the container in which the TRISO particles are packed. - domain_length : float - Length of the container (if cube or cylinder). - domain_radius : float - Radius of the container (if cylinder or sphere). - domain_inner_radius : float - Inner radius of the container (if spherical shell). - domain_center : Iterable of float - Cartesian coordinates of the center of the container. - n_particles : int - Number of TRISO particles to pack in the domain. Exactly one of - 'n_particles' and 'packing_fraction' should be specified -- the other - will be calculated. - packing_fraction : float - Packing fraction of particles. Exactly one of 'n_particles' and - 'packing_fraction' should be specified -- the other will be calculated. - initial_packing_fraction : float, optional - Packing fraction used to initialize the configuration of particles in + Outer radius of spheres. + region : openmc.Region + Container in which the spheres are packed. Supported shapes are + rectangular prism, cylinder, sphere, and spherical shell. + pf : float + Packing fraction of the spheres. One of 'pf' and 'num_spheres' must + be specified; the other will be calculated. If both are specified, 'pf' + takes precedence over 'num_spheres'. + num_spheres : int + Number of spheres to pack in the domain. One of 'num_spheres' and 'pf' + must be specified; the other will be calculated. + initial_pf : float, optional + Packing fraction used to initialize the configuration of spheres in the domain. Default value is 0.3. It is not recommended to set the initial packing fraction much higher than 0.3 as the random sequential packing algorithm becomes prohibitively slow as it approaches its limit (~0.38). contraction_rate : float, optional - Contraction rate of outer diameter. This can affect the speed of the - close random packing algorithm. Default value is 1/400. + Contraction rate of the outer diameter. Higher packing fractions can be + reached using a smaller contraction rate, but the algorithm will take + longer to converge. seed : int, optional RNG seed. @@ -1084,7 +1337,7 @@ def pack_spheres(radius, domain_shape='cylinder', domain_length=None, Notes ----- - The particle configuration is generated using a combination of random + The sphere configuration is generated using a combination of random sequential packing (RSP) and close random packing (CRP). RSP performs better than CRP for lower packing fractions (pf), but it becomes prohibitively slow as it approaches its packing limit (~0.38). CRP can @@ -1092,23 +1345,23 @@ def pack_spheres(radius, domain_shape='cylinder', domain_length=None, If the desired pf is below some threshold for which RSP will be faster than CRP ('initial_packing_fraction'), only RSP is used. If a higher pf is - required, particles with a radius smaller than the desired final radius + required, spheres with a radius smaller than the desired final radius (and therefore with a smaller pf) are initialized within the domain using - RSP. This initial configuration of particles is then used as a starting + RSP. This initial configuration of spheres is then used as a starting point for CRP using Jodrey and Tory's algorithm [1]_. - In RSP, particle centers are placed one by one at random, and placement - attempts for a particle are made until the particle is not overlapping any + In RSP, sphere centers are placed one by one at random, and placement + attempts for a sphere are made until the sphere is not overlapping any others. This implementation of the algorithm uses a mesh over the domain - to speed up the nearest neighbor search by only searching for a particle's + to speed up the nearest neighbor search by only searching for a sphere's neighbors within that mesh cell. - In CRP, each particle is assigned two diameters, an inner and an outer, + In CRP, each sphere is assigned two diameters, an inner and an outer, which approach each other during the simulation. The inner diameter, defined as the minimum center-to-center distance, is the true diameter of - the particles and defines the pf. At each iteration the worst overlap - between particles based on outer diameter is eliminated by moving the - particles apart along the line joining their centers. Iterations continue + the spheres and defines the pf. At each iteration the worst overlap + between spheres based on outer diameter is eliminated by moving the + spheres apart along the line joining their centers. Iterations continue until the two diameters converge or until the desired pf is reached. References @@ -1117,90 +1370,60 @@ def pack_spheres(radius, domain_shape='cylinder', domain_length=None, packing of equal spheres", Phys. Rev. A 32 (1985) 2347-2351. """ - - # Check for valid container geometry and dimensions - if domain_shape not in ['cube', 'cylinder', 'sphere', 'spherical shell']: - raise ValueError('Unable to set domain_shape to "{}". Only "cube", ' - '"cylinder", "sphere", and "spherical shell" are ' - 'supported."'.format(domain_shape)) - if not domain_length and domain_shape in ['cube', 'cylinder']: - raise ValueError('"domain_length" must be specified for {} domain ' - 'geometry '.format(domain_shape)) - if not domain_radius and domain_shape in ['cylinder', 'sphere', - 'spherical shell']: - raise ValueError('"domain_radius" must be specified for {} domain ' - 'geometry '.format(domain_shape)) - if not domain_inner_radius and domain_shape in ['spherical shell']: - raise ValueError('"domain_inner_radius" must be specified for {} domain ' - 'geometry '.format(domain_shape)) - - if domain_shape == 'cube': - domain = _CubicDomain(length=domain_length, particle_radius=radius, - center=domain_center) - elif domain_shape == 'cylinder': - domain = _CylindricalDomain(length=domain_length, radius=domain_radius, - particle_radius=radius, center=domain_center) - elif domain_shape == 'sphere': - domain = _SphericalDomain(radius=domain_radius, particle_radius=radius, - center=domain_center) - elif domain_shape == 'spherical shell': - domain = _SphericalShellDomain(radius=domain_radius, - inner_radius=domain_inner_radius, - particle_radius=radius, - center=domain_center) - - # Calculate the packing fraction if the number of particles is specified; - # otherwise, calculate the number of particles from the packing fraction. - if ((n_particles is None and packing_fraction is None) or - (n_particles is not None and packing_fraction is not None)): - raise ValueError('Exactly one of "n_particles" and "packing_fraction" ' - 'must be specified.') - elif packing_fraction is None: - n_particles = int(n_particles) - packing_fraction = 4/3*pi*radius**3*n_particles / domain.volume - elif n_particles is None: - packing_fraction = float(packing_fraction) - n_particles = int(packing_fraction*domain.volume // (4/3*pi*radius**3)) - - # Check for valid packing fractions for each algorithm - if packing_fraction >= 0.64: - raise ValueError('Packing fraction of {} is greater than the ' - 'packing fraction limit for close random ' - 'packing (0.64)'.format(packing_fraction)) - if initial_packing_fraction >= 0.38: - raise ValueError('Initial packing fraction of {} is greater than the ' - 'packing fraction limit for random sequential' - 'packing (0.38)'.format(initial_packing_fraction)) - if initial_packing_fraction > packing_fraction: - initial_packing_fraction = packing_fraction - if packing_fraction > 0.3: - initial_packing_fraction = 0.3 - + # Seed RNG random.seed(seed) - # Calculate the particle radius used in the initial random sequential + # Create container with the correct shape based on the supplied region + domain = _create_container(region, radius) + + # Determine the packing fraction/number of spheres + volume = 4/3*pi*radius**3 + if pf is None and num_spheres is None: + raise ValueError('`pf` or `num_spheres` must be specified.') + elif pf is None: + num_spheres = int(num_spheres) + pf = volume*num_spheres/domain.volume + else: + pf = float(pf) + num_spheres = int(pf*domain.volume//volume) + + # Make sure initial packing fraction is less than packing fraction + if initial_pf > pf: + initial_pf = pf + + # Check packing fraction for close random packing + if pf > MAX_PF_CRP: + raise ValueError('Packing fraction {0} is greater than the limit for ' + 'close random packing, {1}'.format(pf, MAX_PF_CRP)) + + # Check packing fraction for random sequential packing + if initial_pf > MAX_PF_RSP: + raise ValueError('Initial packing fraction {0} is greater than the ' + 'limit for random sequential packing, ' + '{1}'.format(initial_pf, MAX_PF_RSP)) + + # Calculate the sphere radius used in the initial random sequential # packing from the initial packing fraction - initial_radius = (3/4 * initial_packing_fraction * domain.volume / - (pi * n_particles))**(1/3) - domain.particle_radius = initial_radius + initial_radius = (3/4*initial_pf*domain.volume/(pi*num_spheres))**(1/3) + domain.sphere_radius = initial_radius # Recalculate the limits for the initial random sequential packing using - # the desired final particle radius to ensure particles are fully contained + # the desired final sphere radius to ensure spheres are fully contained # within the domain during the close random pack domain.limits = [x + initial_radius - radius for x in domain.limits] - # Generate non-overlapping particles for an initial inner radius using + # Generate non-overlapping spheres for an initial inner radius using # random sequential packing algorithm - particles = _random_sequential_pack(domain, n_particles) + spheres = _random_sequential_pack(domain, num_spheres) - # Use the particle configuration produced in random sequential packing as a - # starting point for close random pack with the desired final particle + # Use the sphere configuration produced in random sequential packing as a + # starting point for close random pack with the desired final sphere # radius - if initial_packing_fraction != packing_fraction: - domain.particle_radius = radius - _close_random_pack(domain, particles, contraction_rate) + if initial_pf != pf: + domain.sphere_radius = radius + _close_random_pack(domain, spheres, contraction_rate) - return particles + domain.center + return spheres + domain.center def create_trisos(outer_radius, fill, centers): """Create TRISO particles at the given coordinates. diff --git a/openmc/region.py b/openmc/region.py index f82db8553..ceffadbe3 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -475,7 +475,7 @@ class Complement(Region): >>> xr = openmc.XPlane(x0=10.0) >>> yl = openmc.YPlane(y0=-10.0) >>> yr = openmc.YPlane(y0=10.0) - >>> inside_box = +xl & -xr & +yl & -yl + >>> inside_box = +xl & -xr & +yl & -yr >>> outside_box = ~inside_box >>> type(outside_box) diff --git a/openmc/surface.py b/openmc/surface.py index 826ac5cc8..c861880c5 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -62,7 +62,7 @@ class Surface(IDManagerMixin): self.boundary_type = boundary_type # A dictionary of the quadratic surface coefficients - # Key - coefficeint name + # Key - coefficient name # Value - coefficient value self._coefficients = {} @@ -1661,7 +1661,7 @@ class Quadric(Surface): a, b, c, d, e, f, g, h, j, k : float, optional coefficients for the surface. All default to 0. name : str, optional - Name of the sphere. If not specified, the name will be the empty string. + Name of the surface. If not specified, the name will be the empty string. Attributes ---------- From 66fe2185f37a0a98ba66a3ddd31caed3ecb0a7a5 Mon Sep 17 00:00:00 2001 From: amandalund Date: Wed, 14 Nov 2018 10:54:24 -0600 Subject: [PATCH 5/9] Update TRISO tests --- openmc/model/triso.py | 26 +-- tests/regression_tests/triso/inputs_true.dat | 14 +- tests/regression_tests/triso/results_true.dat | 2 +- tests/regression_tests/triso/test.py | 16 +- tests/unit_tests/test_model_triso.py | 193 +++++++++--------- 5 files changed, 121 insertions(+), 130 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index 03ca6f5e0..a99f136bb 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -882,21 +882,21 @@ def _create_container(region, sphere_radius): if p1.surface.x0 > p2.surface.x0: p1, p2 = p2, p1 length = p2.surface.x0 - p1.surface.x0 - center = (length/2, cyl.surface.y0, cyl.surface.z0) + center = (p1.surface.x0 + length/2, cyl.surface.y0, cyl.surface.z0) # Calculate the parameters for a cylinder along the y-axis elif axis == 'y': if p1.surface.y0 > p2.surface.y0: p1, p2 = p2, p1 length = p2.surface.y0 - p1.surface.y0 - center = (cyl.surface.x0, length/2, cyl.surface.z0) + center = (cyl.surface.x0, p1.surface.y0 + length/2, cyl.surface.z0) # Calculate the parameters for a cylinder along the z-axis else: if p1.surface.z0 > p2.surface.z0: p1, p2 = p2, p1 length = p2.surface.z0 - p1.surface.z0 - center = (cyl.surface.x0, cyl.surface.y0, length/2) + center = (cyl.surface.x0, cyl.surface.y0, p1.surface.z0 + length/2) # Make sure the half-spaces are on the correct side of the surfaces if cyl.side != '-' or p1.side != '+' or p2.side != '-': @@ -1424,23 +1424,3 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, _close_random_pack(domain, spheres, contraction_rate) return spheres + domain.center - -def create_trisos(outer_radius, fill, centers): - """Create TRISO particles at the given coordinates. - - Parameters - ---------- - outer_radius : float - Outer radius of TRISO particle - fill : openmc.Universe - Universe which contains all layers of the TRISO particle - centers : Iterable of float - Cartesian coordinates of the centers of the TRISO particles in cm - - Returns - ------- - list of openmc.model.TRISO - List of TRISO particles in the domain. - - """ - return [TRISO(outer_radius, fill, c) for c in centers] diff --git a/tests/regression_tests/triso/inputs_true.dat b/tests/regression_tests/triso/inputs_true.dat index a13187628..7e5007eca 100644 --- a/tests/regression_tests/triso/inputs_true.dat +++ b/tests/regression_tests/triso/inputs_true.dat @@ -5,7 +5,7 @@ - + @@ -220,12 +220,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/triso/results_true.dat b/tests/regression_tests/triso/results_true.dat index dfc700c05..2bbecb4a9 100644 --- a/tests/regression_tests/triso/results_true.dat +++ b/tests/regression_tests/triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.683226E+00 7.383559E-02 +1.683227E+00 7.383566E-02 diff --git a/tests/regression_tests/triso/test.py b/tests/regression_tests/triso/test.py index 174bba94f..9ede6dbb0 100644 --- a/tests/regression_tests/triso/test.py +++ b/tests/regression_tests/triso/test.py @@ -53,19 +53,21 @@ class TRISOTestHarness(PyAPITestHarness): c5 = openmc.Cell(fill=opyc, region=+spheres[3]) inner_univ = openmc.Universe(cells=[c1, c2, c3, c4, c5]) - outer_radius = 422.5*1e-4 - trisos = openmc.model.pack_trisos( - radius=outer_radius, fill=inner_univ, domain_shape='cube', - domain_length=1., domain_center=(0., 0., 0.), n_particles=100) - - # Define box to contain lattice + # Define box to contain lattice and to pack TRISO particles in 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) + box_region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + box = openmc.Cell(region=box_region) + + outer_radius = 422.5*1e-4 + centers = openmc.model.pack_spheres(radius=outer_radius, + region=box_region, num_spheres=100) + trisos = [openmc.model.TRISO(outer_radius, inner_univ, c) + for c in centers] # Create lattice ll, ur = box.region.bounding_box diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index 19c4f9081..3c261ed6f 100644 --- a/tests/unit_tests/test_model_triso.py +++ b/tests/unit_tests/test_model_triso.py @@ -10,19 +10,54 @@ import pytest import scipy.spatial +_RADIUS = 0.1 _PACKING_FRACTION = 0.35 -_RADIUS = 4.25e-2 -domain_params = [ - {'shape': 'cube', 'length': 0.75, 'radius': 0., 'volume': 0.75**3}, - {'shape': 'cylinder', 'length': 0.5, 'radius': 0.5, 'volume': 0.5*pi*0.5**2}, - {'shape': 'sphere', 'length': 0., 'radius': 0.5, 'volume': 4/3*pi*0.5**3} -] +_SHAPES = ['rectangular_prism', 'cylinder', 'sphere', 'spherical_shell'] +_VOLUMES = [1.**3, 1.*pi*1.**2, 4/3*pi*1.**3, 4/3*pi*(1.**3 - 0.5**3)] -@pytest.fixture(scope='module', params=domain_params, - ids=['cube', 'cylinder', 'sphere']) -def domain(request): - return request.param +@pytest.fixture(scope='module') +def container(request): + return request.getfixturevalue(request.param) + + +@pytest.fixture(scope='module') +def centers(request): + container = request.getfixturevalue(request.param) + return openmc.model.pack_spheres(radius=_RADIUS, region=container, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def rectangular_prism(): + min_x = openmc.XPlane(x0=-0.5) + max_x = openmc.XPlane(x0=0.5) + min_y = openmc.YPlane(y0=-0.5) + max_y = openmc.YPlane(y0=0.5) + min_z = openmc.ZPlane(z0=-0.5) + max_z = openmc.ZPlane(z0=0.5) + return +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + + +@pytest.fixture(scope='module') +def cylinder(): + cylinder = openmc.ZCylinder(R=1.) + min_z = openmc.ZPlane(z0=-0.5) + max_z = openmc.ZPlane(z0=0.5) + return +min_z & -max_z & -cylinder + + +@pytest.fixture(scope='module') +def sphere(): + sphere = openmc.Sphere(R=1.) + return -sphere + + +@pytest.fixture(scope='module') +def spherical_shell(): + sphere = openmc.Sphere(R=1.) + inner_sphere = openmc.Sphere(R=0.5) + return -sphere & +inner_sphere @pytest.fixture(scope='module') @@ -33,73 +68,68 @@ def triso_universe(): return univ -@pytest.fixture(scope='module') -def trisos(domain, triso_universe): - trisos = openmc.model.pack_trisos( - radius=_RADIUS, - fill=triso_universe, - domain_shape=domain['shape'], - domain_length=domain['length'], - domain_radius=domain['radius'], - domain_center=(0., 0., 0.), - initial_packing_fraction=0.2, - packing_fraction=_PACKING_FRACTION - ) - return trisos - - -def test_overlap(trisos): - """Check that no TRISO particles overlap.""" - centers = [t.center for t in trisos] - +@pytest.mark.parametrize('centers', _SHAPES, indirect=True) +def test_overlap(centers): + """Check that none of the spheres in the packed configuration overlap.""" # Create KD tree for quick nearest neighbor search tree = scipy.spatial.cKDTree(centers) - # Find distance to nearest neighbor for all particles + # Find distance to nearest neighbor for all spheres d = tree.query(centers, k=2)[0] - # Get the smallest distance between any two particles + # Get the smallest distance between any two spheres d_min = min(d[:, 1]) assert d_min > 2*_RADIUS or d_min == pytest.approx(2*_RADIUS) -def test_contained(trisos, domain): - """Make sure all particles are entirely contained within the domain.""" - if domain['shape'] == 'cube': - x = max(np.hstack([abs(t.center) for t in trisos])) + _RADIUS - assert x < 0.5*domain['length'] or x == pytest.approx(0.5*domain['length']) +@pytest.mark.parametrize('centers,shape', zip(_SHAPES, _SHAPES), + indirect=['centers']) +def test_contained(centers, shape): + """Make sure all spheres are entirely contained within the domain.""" + if shape == 'rectangular_prism': + x = np.amax(abs(centers)) + _RADIUS + assert x < 0.5 or x == pytest.approx(0.5) - elif domain['shape'] == 'cylinder': - r = max([norm(t.center[0:2]) for t in trisos]) + _RADIUS - z = max([abs(t.center[2]) for t in trisos]) + _RADIUS - assert r < domain['radius'] or r == pytest.approx(domain['radius']) - assert z < 0.5*domain['length'] or z == pytest.approx(0.5*domain['length']) + elif shape == 'cylinder': + r = max(np.linalg.norm(centers[:,0:2], axis=1)) + _RADIUS + z = max(abs(centers[:,2])) + _RADIUS + assert r < 1. or r == pytest.approx(1.) + assert z < 0.5 or z == pytest.approx(0.5) - elif domain['shape'] == 'sphere': - r = max([norm(t.center) for t in trisos]) + _RADIUS - assert r < domain['radius'] or r == pytest.approx(domain['radius']) + elif shape == 'sphere': + r = max(np.linalg.norm(centers, axis=1)) + _RADIUS + assert r < 1. or r == pytest.approx(1.) + + elif shape == 'spherical_shell': + d = np.linalg.norm(centers, axis=1) + r_max = max(d) + _RADIUS + r_min = min(d) - _RADIUS + assert r_max < 1. or r_max == pytest.approx(1.) + assert r_min > 0.5 or r_min == pytest.approx(0.5) -def test_packing_fraction(trisos, domain): +@pytest.mark.parametrize('centers,volume', zip(_SHAPES, _VOLUMES), + indirect=['centers']) +def test_packing_fraction(centers, volume): """Check that the actual PF is close to the requested PF.""" - pf = len(trisos)*4/3*pi*_RADIUS**3/domain['volume'] + pf = len(centers) * 4/3 * pi *_RADIUS**3 / volume assert pf == pytest.approx(_PACKING_FRACTION, rel=1e-2) -def test_n_particles(triso_universe): - """Check that the function returns the correct number of particles""" - trisos = openmc.model.pack_trisos( - radius=_RADIUS, fill=triso_universe, domain_shape='cube', - domain_length=1.0, n_particles=800 +@pytest.mark.parametrize('container', _SHAPES, indirect=True) +def test_num_spheres(container): + """Check that the function returns the correct number of spheres""" + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=container, num_spheres=50 ) - assert len(trisos) == 800 + assert len(centers) == 50 -def test_triso_lattice(triso_universe): - trisos = openmc.model.pack_trisos( - radius=_RADIUS, fill=triso_universe, domain_shape='cube', - domain_length=1.0, domain_center=(0., 0., 0.), packing_fraction=0.2 +def test_triso_lattice(triso_universe, rectangular_prism): + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=rectangular_prism, pf=0.2 ) + trisos = [openmc.model.TRISO(_RADIUS, triso_universe, c) for c in centers] lower_left = np.array((-.5, -.5, -.5)) upper_right = np.array((.5, .5, .5)) @@ -112,50 +142,29 @@ def test_triso_lattice(triso_universe): ) -def test_domain_input(triso_universe): - # Invalid domain shape +def test_container_input(triso_universe): + # Invalid container shape with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, n_particles=100, - domain_shape='circle' + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=+openmc.Sphere(R=1.), num_spheres=100 ) - # Don't specify domain length on a cube - with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, n_particles=100, - domain_shape='cube' - ) - # Don't specify domain radius on a sphere - with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, n_particles=100, - domain_shape='sphere' - ) - -def test_packing_fraction_input(triso_universe): - # Provide neither packing fraction nor number of particles + +def test_packing_fraction_input(sphere): + # Provide neither packing fraction nor number of spheres with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10 - ) - # Provide both packing fraction and number of particles - with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10, n_particles=100, packing_fraction=0.2 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=sphere ) + # Specify a packing fraction that is too high for CRP with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10, packing_fraction=1 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=sphere, pf=1. ) + # Specify a packing fraction that is too high for RSP with pytest.raises(ValueError): - trisos = openmc.model.pack_trisos( - radius=1, fill=triso_universe, domain_shape='cube', - domain_length=10, packing_fraction=0.5, - initial_packing_fraction=0.4 + centers = openmc.model.pack_spheres( + radius=_RADIUS, region=sphere, pf=0.5, initial_pf=0.4 ) From bbaad4ed93f0c6351d31aebe69afeee601815681 Mon Sep 17 00:00:00 2001 From: amandalund Date: Wed, 14 Nov 2018 16:03:44 -0600 Subject: [PATCH 6/9] Updated TRISO example notebook --- docs/source/pythonapi/index.rst | 2 +- docs/source/pythonapi/model.rst | 2 +- examples/jupyter/triso.ipynb | 115 +++++++++++++++++--------------- 3 files changed, 63 insertions(+), 56 deletions(-) diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 3c28a2185..8abb52528 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -22,7 +22,7 @@ there are many substantial benefits to using the Python API, including: - Ability to plot individual universes as geometry is being created - A :math:`k_\text{eff}` search function (:func:`openmc.search_for_keff`) - Random sphere packing for generating TRISO particle locations - (:func:`openmc.model.pack_trisos`) + (:func:`openmc.model.pack_spheres`) - Ability to create materials based on natural elements or uranium enrichment For those new to Python, there are many good tutorials available online. We diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 9ed77f366..90e8b779d 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -37,7 +37,7 @@ Functions :template: myfunction.rst openmc.model.create_triso_lattice - openmc.model.pack_trisos + openmc.model.pack_spheres Model Container --------------- diff --git a/examples/jupyter/triso.ipynb b/examples/jupyter/triso.ipynb index bc7801e2d..36e0c1f14 100644 --- a/examples/jupyter/triso.ipynb +++ b/examples/jupyter/triso.ipynb @@ -101,7 +101,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now that we have a universe that can be used for each TRISO particle, we need to randomly select locations. In this example, we will select locations at random within a 1 cm x 1 cm x 1 cm box centered at the origin with a packing fraction of 30%. Note that `pack_trisos` can handle up to the theoretical maximum of 60% (it will just be slow)." + "Next, we need a region to pack the TRISO particles in. We will use a 1 cm x 1 cm x 1 cm box centered at the origin." ] }, { @@ -111,16 +111,49 @@ "collapsed": false }, "outputs": [], + "source": [ + "min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective')\n", + "max_x = openmc.XPlane(x0=0.5, boundary_type='reflective')\n", + "min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective')\n", + "max_y = openmc.YPlane(y0=0.5, boundary_type='reflective')\n", + "min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective')\n", + "max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective')\n", + "region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we need to randomly select locations for the TRISO particles. In this example, we will select locations at random within the box with a packing fraction of 30%. Note that `pack_spheres` can handle up to the theoretical maximum of 60% (it will just be slow)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], "source": [ "outer_radius = 425.*1e-4\n", - "\n", - "trisos = openmc.model.pack_trisos(\n", - " radius=outer_radius,\n", - " fill=triso_univ,\n", - " domain_shape='cube',\n", - " domain_length=1,\n", - " packing_fraction=0.3\n", - ")" + "centers = openmc.model.pack_spheres(radius=outer_radius, region=region, pf=0.3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have the locations of the TRISO particles determined and a universe that can be used for each particle, we can create the TRISO particles." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "trisos = [openmc.model.TRISO(outer_radius, triso_univ, c) for c in centers]" ] }, { @@ -132,7 +165,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -142,10 +175,10 @@ "output_type": "stream", "text": [ "Cell\n", - "\tID =\t10005\n", + "\tID =\t6\n", "\tName =\t\n", - "\tFill =\t10000\n", - "\tRegion =\t-10004\n", + "\tFill =\t1\n", + "\tRegion =\t-11\n", "\tRotation =\tNone\n", "\tTranslation =\t[-0.33455672 0.31790187 0.24135378]\n", "\n" @@ -165,7 +198,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -175,7 +208,7 @@ "output_type": "stream", "text": [ "[-0.45718713 -0.45730405 -0.45725048]\n", - "[ 0.45705454 0.45743843 0.45741142]\n" + "[0.45705454 0.45743843 0.45741142]\n" ] } ], @@ -194,7 +227,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -205,7 +238,7 @@ "0.2996893513959326" ] }, - "execution_count": 7, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } @@ -218,42 +251,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Now that we have our TRISO particles created, we need to place them in a lattice to provide optimal tracking performance in OpenMC. We'll start by creating a box that the lattice will be placed within." + "Now that we have our TRISO particles created, we need to place them in a lattice to provide optimal tracking performance in OpenMC. We can use the box we created above to place the lattice in. Actually creating a lattice containing TRISO particles can be done with the `model.create_triso_lattice()` function. This function requires that we give it a list of TRISO particles, the lower-left coordinates of the lattice, the pitch of each lattice cell, the overall shape of the lattice (number of cells in each direction), and a background material." ] }, { "cell_type": "code", - "execution_count": 8, - "metadata": { - "collapsed": false, - "scrolled": true - }, - "outputs": [], - "source": [ - "min_x = openmc.XPlane(x0=-0.5, boundary_type='reflective')\n", - "max_x = openmc.XPlane(x0=0.5, boundary_type='reflective')\n", - "min_y = openmc.YPlane(y0=-0.5, boundary_type='reflective')\n", - "max_y = openmc.YPlane(y0=0.5, boundary_type='reflective')\n", - "min_z = openmc.ZPlane(z0=-0.5, boundary_type='reflective')\n", - "max_z = openmc.ZPlane(z0=0.5, boundary_type='reflective')\n", - "box = openmc.Cell(region=+min_x & -max_x & +min_y & -max_y & +min_z & -max_z)" - ] - }, - { - "cell_type": "markdown", + "execution_count": 10, "metadata": {}, - "source": [ - "Our last step is to actually create a lattice containing TRISO particles which can be done with `model.create_triso_lattice()` function. This function requires that we give it a list of TRISO particles, the lower-left coordinates of the lattice, the pitch of each lattice cell, the overall shape of the lattice (number of cells in each direction), and a background material." - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "collapsed": false - }, "outputs": [], "source": [ + "box = openmc.Cell(region=region)\n", "lower_left, upper_right = box.region.bounding_box\n", "shape = (3, 3, 3)\n", "pitch = (upper_right - lower_left)/shape\n", @@ -270,7 +277,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 11, "metadata": { "collapsed": true }, @@ -288,14 +295,14 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAALVBMVEX///803tEsYIdnEy2T\nUVA4vPLpgJFyEhJNv8QsP88Otf3sYrFLC4epQXVfq1V9dXMyAAAAAWJLR0QAiAUdSAAAAAd0SU1F\nB+EEAxMJO4hQJO0AACJwSURBVHja7Z07dhpLEEC9BZ/DCmx2IDZguQmEhCJGOxDaAWgHsAAS5URS\nSqaUUCGpXkSsNTympz9V1VX9mRkk+53XmZFsz6W6/tU93779V9b37utH5/XzovP6H+R/kP9B/oMg\nI6XUfwFE6fX7rwdRqjqtGSH5+0BGJ475vCb5u0FOHHfzelWz3381iBrfzR+Wy+X8vlJ/N0g1X+o1\nv0ci+RSQ4cnI9AQyGt89LA0JEkl/IMPLGEecJB/ECWS5XNzPzgCiIg/b/OyyHxAnkOUS7a2eQLRt\nn/Ekw/qHs6hIskGanfW4Xq/13uodRI0rvdiv3Tqwy15A6p114livTnsLKklbEKUAyNDY9nuWxFD2\nA1KrSC2QWiSL+3wQBR8Yf/7Lg9QctfadLGK4gYYWRPUAoqqHRiANyCwXZIAeGH/uCYfju7mziJch\nSLPtZj2DrGptzwVRMRD3A28ST9/RmUHu2oOoFMjA+6haJB1BmFTjjCAKg9QCedQWkRFJGQiTaXya\nRJSzJLVIGJASZdf/43l0JA1yBy1JuLdi5nfYfOpARioukhG2WtnmV+VYrfGdNe3NPx6ARByiCV3y\nQcYZfoR55EGGHxkAkBUH0uwtdmcNKYhKgtx5zz5nPTvvMgasQL5Dz651fQ32Lfe0Mz6itJFmCciD\nJllpffzNgEjqkFwQZM2BRKLfACS1tUgYz4AMlKQPxSCX7PNyHN1A+MSq+dbaiISAzBmQH1I2MuRB\nZI56b82dQM4KwklEXIGyp/yIFomOIxakjAJAqp5A8jn03vqBQEYJjpNIqvnpP5zPsUAAyLjO4jqC\nrEpBfjS27Cd4ojhHXQ+qTO7zmwdRLUG+p/xIxvp5kV5WViMX8lwIIG3tbypE6QfEm4GRqQFc9A5S\nmXpA420LdL0EBLhK1rRF/EiuZxnYgsBqWWi0CkBGtRLHDIEMomO9HBKYWLVTkQwQZQK2JEgQohiV\nmmWQYG+b9eSkiJoDwjVFOBBSMBmMmwrP6e8mQYZOJIvMnUWLqEkQY6qyQOi3rEOBmiQpkiHwtlkC\nGdIwMhekKgYZ2JCmfrYUyA/tox60/LIEomiGcj4Qv+1PoUBKJKcnq7thJ2c7y9MQu08uW4KMspuh\nqMKTATIMHi2+sypSjigDYU2xANIIRHuGRXJvcbs+BdK0M1spO7/NeBDnrFfaEv1KgjQ5UhaG3lm1\nRQSFlSQINL/+r/9Ogoxd+LTK2FuZzw9B5uUg1iGOlP7b2jOkQAaowjOvzgNyXwLi9cI4Bm1PfydB\nfIpRK8kfAOKCxjqxihUfRBCdYvwBIDbkldq6PIgvTJ4RpEhHXGIl1YO+CKTYarnVCORxTdu6XwFS\n+5G7Mj/i5aJrpkxb9xNAhtRTtvDsAERo6wrK3idI6POD+nwBiNTWzTK/IkiTjGVwYK9fHv0CEKFB\nkuEQZT9iChYpkHEw91Ccj7iFWlawP90lRKlrYXVWn9CQcRgXF2eIPMgyCaLsxMdaGzkeZNBEoYlY\n0Wg2+a28nH0Uln1KQSrYqBVUZGCNz2UMRLEgZP2UOUgGJY4+ZCRWkq7XAsFpRb8gI1WFBYeRqwbm\nKDsZZpBUpKpMNB0FsWWlUpCR//fhp+MS81t/20YkcvFhYEaCEnurNUg9hTk//f+4DyJ2Q7+JT6kr\nPAu59jCIzjZ13lojWGNCIHw3VKxrNRFq/ZjCbwz8b5wDRMGqHxCJ1vbc6Pd7Y1vvTg9ZiSVTCzKP\ng7B+JAdE6HpKn4uTD85LfO8K0hTDC0FGuDLO/WCRB2Lq83K8mLm1uFgrB6RCvQq0txrdIRPMsVmU\neOMnU9lNOHJZCEK6R8huaZIFmSlvP2UKzG/0IX+kBcKBkH4eYrxjuqEdxmWbjDWj4Jsu3YUgtMOK\nDbNxTRe9gVR398moMWcFINrvgZ43ema2e9UBZGBjoctzgJDmPd53TDe0EAQZssbX9SCQOMg6AKlR\nfl90ATGTV04kTQf1C0DCHnURiAs3oHxSqe55QEYmg24FMmgCaxC2NB3tM4OsGJBwjqMExB7Tuied\n3i8ACXskBSCDsQ3tSaf3DCAxP3JheyQLOOtUAOKCnAURyTlAIp79okm6yPRZPsgAdHqrfkBc7MKA\niLGWEQidBywA8YE1yeMznviS+RA0UEMQF/2uafR7wfdI8kFARrMsA+HDX9jSZkBQPkJ0neuRFIOs\nw72V5BhzmZU1PHXMCUBGxg6Jc7EX7uhb8yjlICbrXwe99xQHG5INx8bw1CQ/PYYtyTnTQnJ2Jy3X\nI1GFIO5vr2j1McoxFMY7YJHk8qfjaOLa3xe+6bkIBIJKQtY0F4Ksi0EUn9rj0aifkMNkGk1daxkc\nkxaKdNkgA/c1rEtA3AgUORGHVeCndQ+VHwXgK40MyOxzQOC8FARBRum3eUA9m1wLoZk5YWq/zjQ/\nkvClBKTN1kI6ewn5QE9zZg3VHaqPmFlkytGArNuCtFR2uIUgiPtSfU9zhEIHZ8MCjo4S8c2fEvM7\nFI72KaYVOEKhQ/j8PelIM/p0+vsrySEOGQfebCFXMEQq8kh2Bwwd5gUgpeYX/z+/ApBhs6MvCQg6\no+g+R1G6iaWIx42BdPIjiaCRnwYcgson1HYAsjbf/4hYkyhIF8/uBk/5ML7xAdSBK1zCRSBrCgI2\nDFdu8KtbrAUSKwZk6Bw4mphtB7IKI3cM0in6dSJhU13V+DJ9XUlLEOQdEtreJR8xVetF2MWCDnyJ\nRVKiI2CSJwXSKUO0Veua4xcFERx4idUqAOmWs38HJxi+MyCc3yvxIyUgTRVl2bKKomuk3JkS7MDR\nEZICz16gIx3rWt9ttZd+6r545my+umNFxcRaBVara6WxQQk/s6qwXodn84X7IJjol/qRKEh4Yqyf\nC178F0xB/N7CKeI4yEdGJOG5uEiRtK/GiyDywTeQIVaJDJGE1xdl69wgP8BICfwY5ewGBEalpLb4\nB4DY63jSVZTsfORrQGxdq+LqWubw8A8LQjLErwCZy4dDhSur4GC5BrlA5f6vARlWkvltSLiRAXdf\nyqUtmZqZBltF+RIQG1KRINc9s3AFx7hC1XhU1/oiEOQWLn9kLZ/kGxBdyZpHz02eGQTW23JPucJl\nQEzk0YajNxD5uO4wY6LAgqRv4cgFGbQ51P6D1qTp9kkOefz0D5S68yEPxM4ylIMIx3XtIGDxmFM3\nkIEt6XCPG+C5D+pHMVeM0mvo9JAz4w3PDNKkKyxJuOcG7s/NV8/N0zl7irLGcLf1DAKmFH9xHMFN\nCr88CBuJDE1tZQEjxqHx9JfnA4G3FNJtFE5qEpAfzNDTEM5zw91GifsFGcCuXQAypsk6BWEiES49\nBN77XCB+9jmYIx9Y4xMDqTcN2vlc5YHdbX2D+O9vzoLgaWas7Nwa4sPkRiDMbusZJHJmxPXKRfPL\nLXKY/FLabT2DkOtN0iBuiSBuwscF90Kdq2eQOaxN/eoMMvR10ZUFIWeezgjyyN4q0BJkjg4BXloh\nPbra1eVZQHAHoCeQtQdRvJDOB8Ic2YufwcgAMRdvcUL6AhDpHqECEFRx93dxfdrW4jx7W5BH/FFv\nIE92bbCyXz3BZWKtJ3Y982tLdeT5WVEQ86svr50XBLnzfvh+hp91ozmuykGAQkxOIHcQZHk2EOgQ\nnxiSp0IQfCnJZ4E84UpI8LiSPESQZ3IGOwA509Z6QmGQ+NQFIMT7nUC2n6EjYG/NewLB8chzA7L8\nBBCQWBVwiCBbHCEyhmxyFpCnOtW1lZC2IFulJl4kd2CzTgJDptXmHCCbsauElOysJ4RRL7C3HqBA\nOP0/B8iTPWqvigTyBDmakMztLTBiivW/8bqe/rpXkI2NDYs4PIg+0Ff/AxP3wZ0tMDUfbfFcOpBi\nVxIEYhx42cbyIFtfpJu4T+wUV6D/Rm2erRR7BYk68CQIMBYTz4Y22zOaJJ3436i6iuRb6VPLIFs4\nq+weHKu//yWrNraVOOsokh5B0KjTxO83yGHFVvfdG4H4ybtuIukXBLs/3ksqM/RlBeJ91x8CAv34\nQgZ5dvk/EmO9Hf8YEDQsN4mQqLHdblt3NV7XvdUfiMI5vzPAWEXIZ1s/C3n/54A8gOkHA2KMliie\nLew9/ikgKIyqgLNz7jAB0klJzgFiSw2hYw/1xZvsjtqeB6KiDr95prFPYy2I9xmS7n8yiDk7flUI\nwrr6rwRx476FWwsVAXiRfKqyK39heR7Iqsliw0z3a0E2yp9Fu8oFqc3vdgw6b/ciyKf5Ed/1uxdq\nEtSPLBsQ9IGwt7YwRDk3CLyw/EoG2YJXLGidEPJzqu137rjG6+uufa6YBNmgC8ujIKjkSiqmUhip\n7IxmnZDsOmS9SRBUZeP3ln0mrNtb2lXglcSnx9c7O9PThiQNAuqeQZEeg2Bri6uKy3vJk7irlF79\nfYlnANmQrt+VDLKFL6qcoDqvvqPlWRKJydl3QDjnAMFFzgjIsz9nVcdWmSDPuhl2ihtegbrkgOyw\nMiVBcNmZVRKkuK4cBFsIHkSF+Umj4a+v8LxZDkbzt3JBgmYZ8ztb+3R+vz8LIKp5XVZIcno0kCwm\n99YuqIYlQWiPiRGaObBU7/exTz+4raVomcusl/rZYFcjBRLahRKQNQuycTGlK39OBJCtvtX7HtRU\nAQiqJiVAdn4TX7cECe2vv85+gtPx0PyCoTMWBLihhEBCdXIgQuqUBNk0tkbfNojtauAQ68BKJ1r0\nV0tBdkyo+c3tc756ndxa8P7HCQahIcoWHIOcdAHhfvUb3OezkCTowwYCgadXsDGiQSM6H9IBZIej\nOgjizydSkk3C/G7QQWb0PQdhvD9hS5lLQdCMxrUH0fv8oTkmxIDEPDu4tZaG6iaxegThsCS8MhB8\nfgiCwOyGPOmGtC+DH+PD/kRJSDgcNHlagpDqgAPZgMQ5iEHgMcfl/YzjBKcgsZL4XqgOh3NBFvN4\n/IvrNfceBB3YpyC4fRnqOj6XSgxwZZs6pjpEQbbG67wgFU45RHKG3IPgV3rQ7xwdc4yo0CpMzHUv\nVL9xTDvLMbULLkh7EZyDrOvA084siN8e+llZLbD/Pavra69CGASXTEMQF9BokNxWCQRZIxByp7zk\n8biZIehmmAoD6oWGW8s1GF+0Dmc2rwKQawNCX94eqIHv+l8VgqBRiABka45yGxAmEiyTiEIXSXDb\nx3X9n0pB9Lc+CczxvJkVMvebTcwIB0jfW4HQqz0oyEYH35zbzwHB5hg5RGVmVy3ILq/lviOXGluQ\nDQZchPsnMtkRV/YQxMZaBqRCIK8ufY8vej0oB7Ji6wsbccY0bn4pCBgWmHAgrzQR55fiHWIGiM5W\nnriFHOJCrJR4B2nL4RNma70GpREZBNT9i0CktSE3REVBnq3ZsFk9UfbctcM3TFz3AfIUCRoZkSib\nLUyMbYbmNx8EX17CgrDKHluRMJ7VEtjl3WKHmL9wlfnVgCTNb2JviYkVKxJlNpP1MT5EKQGBVeZr\nB5JyiCklkbKl+kEnIQmYH0BBY9neCmp56RAlubeE4oP/vumnmGtSDvIKj36+OpBE0JgUiVAO8l4i\ntpqfloIwhftkGJ8WCSzQYY7o0ANcpSCv7pys5YCJ1YpNrNIicRcrMZaWHXrYKq5kWkhi7ih5hSDR\nVDeHxLoDZGj9u1bJUzO6Uw4SRDO0+FA2TN6src853JcO37VKbUAFTXBbkBPKNfxTqhyUs5hNf0pr\nH3irvHXv2ewKgleqQNcSBL9taoI4jL2BnxqQXYcx5lTJtB3I1oZgKxq5mMR8ia31i9/218Ej5nXe\nU0XstiD4ijkoKs5/vpgHrm0GeehdXo6SbCu0A6GX/hFRWX1EIDubhxIO7lMZRGz0tAQhUQ8UyJIx\nAjWIO+pwjTgqezFILkiHxYOwNQl8p/cEgOz8peRQPzJLRDJIiXxCEHpVqX1iyQi8CFUtV31M994F\nkE2JwvAgMFfzIOA9xWBvvTQB7QOde3D3lqdb1gJIkeYXgJD3FCMQP//hQZhUsBwke2/lg0jW7EXv\nrKCEvcNXuv4lIGC01u8sdDL9+mtAOGWnL1xmQOAkCukLE5CsoZqOIJL5lfzLC5xE8SC0Ux9wKJxY\nCSBXvYCgSmpUIg7EfvU7OjuBOXCd+DzmV9KFqI4EEuF7bNbDjHGuKzvEbI6CoFE2v0hH7MPS+SIo\nkLEtCSRAShYHgiqpkxTgix87W2SCVBWOA3oA2TAndoTEaosuXvaAL2wwQkCg2dopd+pP9QWyYctw\nwlE+3AkHIFzNLQLiG6dVJsgmaYhVWHqAo7M401X8DM5L82hNT949b7C1AEjlizRZINbIATCWlBEJ\nWw7aopcJPAMQUDx8lUCgitCzAXGQjW1PXnmyECSsaqECHfyRUCZ6eeV6urvxg2R+YcacAbLxtuHK\ngVzh37Hvhg3UfcyWTHlJaRB7Qu7ag0gOcee2aB4InEWQJCLWqvkfbIG5eUYgYRV0h0OU9iAbL8DK\nloQ32SBSW2GM2m8AJKyXkKON161Bckqpke7BNrRmbsi5YkqmO1r3IcOErUGwMygHATykuZNZxHbD\ni813iRD9XkmCWIGslrFLeJIgKnA0DQeUlwgipbrQZmSAoKvHRZBxUFwPdT7s9yA8AYTtekLPvszy\n7Bvyag3JIrB+BFgpeKEAVhX31yQQNgFufuCvz0+C0ORBBIkd9CZXPAgqL4HwRaJXHwf4z2MgJM0T\ntT0C4o/r4l4pcZciyM6/mYEQOiueBqEvBJJ2IBc0eoGAQ4well4FI4Igf78DWbqb7cpIrMaZgwQb\nmQPNMEMQ8ObTOAhoNuyM678GhLOc8yO0piOBiPdr+avCtAUHe4u2eyIg3k3iljQ9fJkAyZpRiYDg\nQ4wTAOii8MlzXg8RnKLxgFlnrHK3lgyyJYcYAxCrORkgNhdeeEOF4pnuyh4FCW5qdCqCspIMEFeY\nF+abu5vfGAg9xEiNcj6ILxcJo/MJkLRDTIE8epB7AALC6kwQNKKfAbJxxQaFpoaqFiDSRZktQJr6\n0tpelZgGUc3FPrViZ08NZYHA47poa+WB+FP+QochuPGscg13E8av+eMYvYFQHQnSK6MiyP6lQDa+\ntfqUPzVUDgITPAJCHR1QEfcPLaoUCHgD9xU6OhK9S7MY5BkmeBjE3TJAQOjdoXEQcgnoZgy+uKuO\nIFDZg6LQi+dwR2cCXX/MB0Hnfq8aAdnC8lMLEMn82psoF0GsJQ0IBF9JAoTc+Lvx0fZVKxBQKsR9\nEjIoaEFATR6LJNikURB8mbRW93oU5D5ybVAKBB1i9G0Ed6khyUcU10ksByGXr1092RpJoi8qgjz7\nQ4wrkpCYvIhkiHKxoQwEF/b0bhJPjmSBcHewmp+gKtKL8xZS+adIR2g32OAl29QyyBYfYnz2l1fg\nMp0DQdcmtbVaG1r8TgFkgJCDbmDIH42mvvjHBXsb7K0iP7Kh3eDuIPR8mDt2MdFCcb9mQcRJhyLP\nvgHGcp0LUseXERB8B2v0zodXGBmucPW9LNba0E5XBkcd9MdAnsEhxsnWXznGXJXwGh3ZUAXRbwBy\nlSOQk0hiIO64rqpvGADZKgOyi4FY15CRj7SQyCYJAu502Y7ppeVxEPS4cHQrlSG20JEMEN9W2MJs\ntRhE5efsLaxWtFwKUSbPdhpiba8cK9tabBVlp9jpoBZ+JEcibsFslZy9TCt7M7RB61ogBcPvVsCS\nFYxUkbIDtZeqdc+M+WWcBa00kuovH2v5m0M2MIAP585VyvxCEP69CgAkMfo3ruCTN50J11jA0a+v\n1JpaL3ryjTlbi0muckGEah0Ekbq4ViO0+bu2W814qEZCOB+hb+kxNX3IwSQneSDSS1QAyA695yMM\n1aFOgLNPmvjbm197HEW/vTkfoH/sjrfNbt7Qes9aAQj42T8feh3xA3xE1hHd33P6AIK8obzm5m3v\n7vWqn3yvZbnUXroriL4XKgT5QJe83EZBYFytbgnIHg4W3rwBOwF+WMvyJgQ5aLl3B7Gp7n0cBL6V\nXv8qAnmDA+r1n1zA1/zswZU5AxA72NUN5OgGy+/jO+uDHD+nIHtfjbnRf3hwItij1IKCHFxvsouO\nnL5oXSHUzcXEzsJ3FlGQN32sxmjFHuYS0BJoBUIgB183mnYCOfp2dFQgx2qORroDkDc7NIX2Ui2C\nPTq4TkDgrKQEcsBpd8WCfLh2dFQgH2ike1kbBgry1nhvbMNqEKxdGOQAi9KSSA7Es09ZkA/zerw4\nxwcdvmdATgA31BhTkDkFQWfZYiAg1hJAPvQXmeAgIIsKgzhZNIrvM4DT93yj0Nt1MAjuP0h7i/Rq\nJJCsBW9FWFEQczu8fcY9Suhummzi0b4mCIIccEdI2luke/beBQTOAFAQp2YtQFCPLqokbrNOzwWy\nt4NVViT81loyW+tA7loQQaAB7wYS2VrYjRcpOzqnE1ESZM/fO4JIyo7dOAOyr0Q/QicLRE8C/g8n\nkEMdDHQEgeYX5fYOBDlE0bPTMaKpKBIfhQKOE0k5iOgQ0dd/g7Wdhig01qJjRCLIQaendUjmfkWb\n/GkxyJFcq+hCFO77JrsNZTLtQBqSkzymSCBKtQDBt0gBEEaVcSZlO7HahrUEMQ8+JX9OOnJeSer/\ndUXCeOy3rZKQxKoR0NztvRYg9WZCf2gmRlqAINtzC0CA17UgONV9AydC3gKQLGVnqGo33AbkiFT2\nw4PgG0uxszcB2F7ZQhIWSGB+i0DSsS7/YwVm10SQmQ9bxiY5MVsNcAkOMeJHWknkqJWIFYk1Rfqv\nx0GagNgKYI+5/C+Rke1CkFkM5Fj/nywJ2OkfaZDTw9/AP6BA3y6faayWJPVLraTVsrv7lkG0V958\niCDhs/pvn/nsYI9JaltYoiJJP0KeljCOQVIMQZYZINx6x3drxEEO+McJEF/dvOVIgDSRH3Ej01Y/\nbjJB0NGfhE7gMt4hGqIcQaWa33iODztEEErtFa8QHMgBjvZNoxxB8Us7ehEEOGFZjyAIql/XEcje\nTrLlgLzDQbIERzNsRH4ggcCwqBikCRr9Oc0ckIM/XRgRiFzGk0FAmTYLhGRRLswKCtY8yLvnjgkE\nvLMqG8RH3JkgOPmAxwlyQN5tQy8mEHgN+TQH5Ihus8wD2YM3utUNBTAqnRQJsEdRTZdtmwzyUAqC\nY3b8f2aBaEP6HlsRb9MnCIzZ9/i9GHkgqRXz/32CwJh9j4cbU3srE6SSI7IedeTN1uF1VovnpvoB\nicXI/YJoZ64fmo5A9AMSqX4JIDB7KgFxQPQ+815A6JGtaR7IQ5lnJyDkqvxMkKj9PcRKFCII6CS3\nAnnAIDdZIHFPcsDfTh7IEVzQ/tEChF6skgWiaPEtAnLasHkgJJ/9DJCDGxLvEeTDnC2+T0bxHUB8\nAt8IxB1nnqZBlrkgR/f1fJwLZG9uFjAgB7AJ+tMRW3xIdas/ZGVPgezdV2VAwB3fsvktBzGJeVIg\nrc2vOVWm85VmZ3lDORVBpHpkBOQjjyMJUv+fivuVO1/R1iDAdfEgyLMvYp69VX2eA8kIUWAOakF8\nKZYHyYy16hrpLEsIaRDygo9QRXD76h1mf2L9wTexV9H2tK0wJ2RCVYcHwS/4uOF/AeZIB0rGgchj\nBf8gDmPGb+MczfBPHGSPX/ARgthnWtsHzwCJDXp4kKOtx8TLJu7ahNsoyBs+xyeA+g5eHkhk9OYf\n8IQ5YeLRF5Zu4yDwZCWjIv66RZNaZIH4YahKqqIcaf9GEIgrLLlf4kH28Kwrx/mAutxToOwEBAYs\n8N6xdwEkK0s/giFSKzce5M38n0KBjs4dTLVjZ7bNgbRw3U1wUwEk6HEKAvFTtRZXANn78+CsCuGR\nlin2I/4pg2a0UVHCAUDAPY7i5O8RF3miIDaWmrGdBTpkpEHuQs9+oF2EQ9MyDJToH/SIj667LIH4\n0fellZsEohuGUoeEgNRmS4He5BQJhGmHvMsg+NJWuJ38H/xZAO1aEyBNvsH/iAE5KPueCaAiptUZ\nPngmyAzYW+DG6ctPEyCRxYC8+3AYqojuoheAkFfHAXkAElrk6Q1Ex00H804eqMkmjsgHUQTkFgrE\nhsRHehasPQhjtXSuS/pqh0wQ/4wSiJ28D0DW7rdagjyEwycHcxkdesIqvbUOYCRAADnaV14YEJLA\ndgAhr3OceouEnjpHRw5OkkkQk7kfq95AwljLkdA9M050f95R2C4ou62l9A8SRL/irkmcjgFmuw7/\nBD+iOoDsVWxcIMhHYiBxgVQwbD/CLFIEKdIRNBbEKAnJECP2KC4QlCLjEMWHvxSkwGr5po+0t5Lp\nhyWJqzp431UtAOWvTwVB41EESfmRvfEAMZG4KkrJWFO4s3DYzofxGKTEs++t4Y6IBNe1Wi68lW6F\nxAqb35JYi55c4ERyhyqN7dYhGN+FeY0ViFUSG2wVRL+uSCHuLVr7bQtCzS18ax8EASFKST6ivdR9\nDMRV4286gQQOkG/toKCxIEPcu9dmyyC0P9JqhbP6H8pF0V4gtK6Yl7OrG3nCl18dQILTE1JHBNVH\ns6oopy/6a0HwqCK/cupaJ5CbrK11NhCtEVGOrEpjDaKtVkLZzwiSsdwb7tzfkEHirr0/EKLsuSTp\narxSOQ6xJxAhtc0iyeiPGCcxy5uW7eRHhPJPzjqiP4nR79h4O3ntbaTfAsBExQe5IFe6vkUeMi4Q\nc2KiTYjiMscDn3/0CPKW5LDBlioPGkGZiM8/+gTRjjHOYS15cRjvTEmYj5wBJLHkpk2GQGzcMCWd\nrg4cbUH2Y7GNlhYIPB9aMCMXrrT5zRCIazLOS0XSHG9cNI0re353nhxaDNYRhWRtQVBGkKPerghx\nUOjosWvHlQoE5ygtQWAjPuO8Hmr4HOAUPhxYK95XMP1tCwKHkeU7BZwM4JCgO05g2qb2lTLFAsEj\nA+1A6JXySe0GRxnQgXO774pnaD78QY+uILCyFBeInUtsrAK46MjW9g7lhjfISNqB0CEyrNaEC1zC\nIICUH8xvBmkbI94BZC/OJ6pw+JdekNAPCHTIZwA5KKgNVkOQleoLBBxz6gLywM7wHvy8BNxZ+O6N\nQNlbgphbwOddQR4ZEHB//hSCoMI9Nb+tQI7wEFjvIAc4U+RVhFzrQhziF4OwOoKmvKYiiApvR+kd\n5KYLyAFfTiiB4KCxueGlvY6wIPG+IQK5w80vAwLvGJ/JIOial0P09HQMRLRa+1h6ixNG1iGSCzyn\nEghMrIznuW0B4l/DS0FUpAFKGNlYixwxE0FgqmvrEOUg3kDe8iCsSPaEEfWM7FE8Or8tgsAZFZW6\n9EHSdn8RIwXZm1kSTiRmhAFru3/no1cRb5IjIHCQwPSiykGs7w37I6AFxQgEV4Lx0VEOxEeSyCG6\nEMDEY4eghpu/t2qSignjbZmUBaFNXuYA8YFeci6EKGgd0NxPqUhgQygLJKxp7+Hx6ikHModxC/Xk\nPYCQBmkmSNhlgA7pPQrC3nPmQVRbkCO2d21B9uAKgjgITawYkDY6QgY2qdXilZ1ru4NLIVgQH20p\nnOrilXfjZ8ZCIGOhJcLOD9hN4Q7m0/OkRA/44cbasLfxIxGQiEPkO3Hu1y0Ib36ZC9vw3mrl2WMg\n1H3HrJb9C/q3ORBUWgknssnPugskcilrCCJFlOaR6AsF6PfOcrSNfqMgb1Lwm2gp2u8W38M/fc9b\n7fKROIi44i1FZ2ZNALbKKgjD9YkgsVzFfrXjshL9V4DEZzWd3mbfTvV1INEs2IH4N7WV7axPBIku\nr7bu3XllAvnTQFxZYV5lHLX4k0Hybqf6G0BcelHG8eeB5Iz1/x0g6SHyvwak1fof5D8I8i8rFUXT\nyFgUuwAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxNy0wNC0wM1QxNDowOTo1OS0wNTowMEJcIC0AAAAl\ndEVYdGRhdGU6bW9kaWZ5ADIwMTctMDQtMDNUMTQ6MDk6NTktMDU6MDAzAZiRAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAALVBMVEX///803tEsYIdnEy2TUVA4vPLpgJFyEhJNv8QsP88Otf3sYrFLC4epQXVfq1V9dXMyAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+ILDgkoDSxqv5MAACJwSURBVHja7Z07dhpLEEC9BZ/DCmx2IDZguQmEhCJGOxDaAWgHsAAS5URSSqaUUCGpXkSsNTympz9V1VX9mRkk+53XmZFsz6W6/tU93779V9b37utH5/XzovP6H+R/kP9B/oMgI6XUfwFE6fX7rwdRqjqtGSH5+0BGJ475vCb5u0FOHHfzelWz3381iBrfzR+Wy+X8vlJ/N0g1X+o1v0ci+RSQ4cnI9AQyGt89LA0JEkl/IMPLGEecJB/ECWS5XNzPzgCiIg/b/OyyHxAnkOUS7a2eQLRtn/Ekw/qHs6hIskGanfW4Xq/13uodRI0rvdiv3Tqwy15A6p114livTnsLKklbEKUAyNDY9nuWxFD2A1KrSC2QWiSL+3wQBR8Yf/7Lg9QctfadLGK4gYYWRPUAoqqHRiANyCwXZIAeGH/uCYfju7mziJchSLPtZj2DrGptzwVRMRD3A28ST9/RmUHu2oOoFMjA+6haJB1BmFTjjCAKg9QCedQWkRFJGQiTaXyaRJSzJLVIGJASZdf/43l0JA1yBy1JuLdi5nfYfOpARioukhG2WtnmV+VYrfGdNe3NPx6ARByiCV3yQcYZfoR55EGGHxkAkBUH0uwtdmcNKYhKgtx5zz5nPTvvMgasQL5Dz651fQ32Lfe0Mz6itJFmCciDJllpffzNgEjqkFwQZM2BRKLfACS1tUgYz4AMlKQPxSCX7PNyHN1A+MSq+dbaiISAzBmQH1I2MuRBZI56b82dQM4KwklEXIGyp/yIFomOIxakjAJAqp5A8jn03vqBQEYJjpNIqvnpP5zPsUAAyLjO4jqCrEpBfjS27Cd4ojhHXQ+qTO7zmwdRLUG+p/xIxvp5kV5WViMX8lwIIG3tbypE6QfEm4GRqQFc9A5SmXpA420LdL0EBLhK1rRF/EiuZxnYgsBqWWi0CkBGtRLHDIEMomO9HBKYWLVTkQwQZQK2JEgQohiVmmWQYG+b9eSkiJoDwjVFOBBSMBmMmwrP6e8mQYZOJIvMnUWLqEkQY6qyQOi3rEOBmiQpkiHwtlkCGdIwMhekKgYZ2JCmfrYUyA/tox60/LIEomiGcj4Qv+1PoUBKJKcnq7thJ2c7y9MQu08uW4KMspuhqMKTATIMHi2+sypSjigDYU2xANIIRHuGRXJvcbs+BdK0M1spO7/NeBDnrFfaEv1KgjQ5UhaG3lm1RQSFlSQINL/+r/9Ogoxd+LTK2FuZzw9B5uUg1iGOlP7b2jOkQAaowjOvzgNyXwLi9cI4Bm1PfydBfIpRK8kfAOKCxjqxihUfRBCdYvwBIDbkldq6PIgvTJ4RpEhHXGIl1YO+CKTYarnVCORxTdu6XwFS+5G7Mj/i5aJrpkxb9xNAhtRTtvDsAERo6wrK3idI6POD+nwBiNTWzTK/IkiTjGVwYK9fHv0CEKFBkuEQZT9iChYpkHEw91Ccj7iFWlawP90lRKlrYXVWn9CQcRgXF2eIPMgyCaLsxMdaGzkeZNBEoYlY0Wg2+a28nH0Uln1KQSrYqBVUZGCNz2UMRLEgZP2UOUgGJY4+ZCRWkq7XAsFpRb8gI1WFBYeRqwbmKDsZZpBUpKpMNB0FsWWlUpCR//fhp+MS81t/20YkcvFhYEaCEnurNUg9hTk//f+4DyJ2Q7+JT6krPAu59jCIzjZ13lojWGNCIHw3VKxrNRFq/ZjCbwz8b5wDRMGqHxCJ1vbc6Pd7Y1vvTg9ZiSVTCzKPg7B+JAdE6HpKn4uTD85LfO8K0hTDC0FGuDLO/WCRB2Lq83K8mLm1uFgrB6RCvQq0txrdIRPMsVmUeOMnU9lNOHJZCEK6R8huaZIFmSlvP2UKzG/0IX+kBcKBkH4eYrxjuqEdxmWbjDWj4Jsu3YUgtMOKDbNxTRe9gVR398moMWcFINrvgZ43ema2e9UBZGBjoctzgJDmPd53TDe0EAQZssbX9SCQOMg6AKlRfl90ATGTV04kTQf1C0DCHnURiAs3oHxSqe55QEYmg24FMmgCaxC2NB3tM4OsGJBwjqMExB7Tuied3i8ACXskBSCDsQ3tSaf3DCAxP3JheyQLOOtUAOKCnAURyTlAIp79okm6yPRZPsgAdHqrfkBc7MKAiLGWEQidBywA8YE1yeMznviS+RA0UEMQF/2uafR7wfdI8kFARrMsA+HDX9jSZkBQPkJ0neuRFIOsw72V5BhzmZU1PHXMCUBGxg6Jc7EX7uhb8yjlICbrXwe99xQHG5INx8bw1CQ/PYYtyTnTQnJ2Jy3XI1GFIO5vr2j1McoxFMY7YJHk8qfjaOLa3xe+6bkIBIJKQtY0F4Ksi0EUn9rj0aifkMNkGk1daxkckxaKdNkgA/c1rEtA3AgUORGHVeCndQ+VHwXgK40MyOxzQOC8FARBRum3eUA9m1wLoZk5YWq/zjQ/kvClBKTN1kI6ewn5QE9zZg3VHaqPmFlkytGArNuCtFR2uIUgiPtSfU9zhEIHZ8MCjo4S8c2fEvM7FI72KaYVOEKhQ/j8PelIM/p0+vsrySEOGQfebCFXMEQq8kh2Bwwd5gUgpeYX/z+/ApBhs6MvCQg6o+g+R1G6iaWIx42BdPIjiaCRnwYcgson1HYAsjbf/4hYkyhIF8/uBk/5ML7xAdSBK1zCRSBrCgI2DFdu8KtbrAUSKwZk6Bw4mphtB7IKI3cM0in6dSJhU13V+DJ9XUlLEOQdEtreJR8xVetF2MWCDnyJRVKiI2CSJwXSKUO0Veua4xcFERx4idUqAOmWs38HJxi+MyCc3yvxIyUgTRVl2bKKomuk3JkS7MDREZICz16gIx3rWt9ttZd+6r545my+umNFxcRaBVara6WxQQk/s6qwXodn84X7IJjol/qRKEh4YqyfC178F0xB/N7CKeI4yEdGJOG5uEiRtK/GiyDywTeQIVaJDJGE1xdl69wgP8BICfwY5ewGBEalpLb4B4DY63jSVZTsfORrQGxdq+LqWubw8A8LQjLErwCZy4dDhSur4GC5BrlA5f6vARlWkvltSLiRAXdfyqUtmZqZBltF+RIQG1KRINc9s3AFx7hC1XhU1/oiEOQWLn9kLZ/kGxBdyZpHz02eGQTW23JPucJlQEzk0YajNxD5uO4wY6LAgqRv4cgFGbQ51P6D1qTp9kkOefz0D5S68yEPxM4ylIMIx3XtIGDxmFM3kIEt6XCPG+C5D+pHMVeM0mvo9JAz4w3PDNKkKyxJuOcG7s/NV8/N0zl7irLGcLf1DAKmFH9xHMFNCr88CBuJDE1tZQEjxqHx9JfnA4G3FNJtFE5qEpAfzNDTEM5zw91GifsFGcCuXQAypsk6BWEiES49BN77XCB+9jmYIx9Y4xMDqTcN2vlc5YHdbX2D+O9vzoLgaWas7Nwa4sPkRiDMbusZJHJmxPXKRfPLLXKY/FLabT2DkOtN0iBuiSBuwscF90Kdq2eQOaxN/eoMMvR10ZUFIWeezgjyyN4q0BJkjg4BXlohPbra1eVZQHAHoCeQtQdRvJDOB8Ic2YufwcgAMRdvcUL6AhDpHqECEFRx93dxfdrW4jx7W5BH/FFvIE92bbCyXz3BZWKtJ3Y982tLdeT5WVEQ86svr50XBLnzfvh+hp91ozmuykGAQkxOIHcQZHk2EOgQnxiSp0IQfCnJZ4E84UpI8LiSPESQZ3IGOwA509Z6QmGQ+NQFIMT7nUC2n6EjYG/NewLB8chzA7L8BBCQWBVwiCBbHCEyhmxyFpCnOtW1lZC2IFulJl4kd2CzTgJDptXmHCCbsauElOysJ4RRL7C3HqBAOP0/B8iTPWqvigTyBDmakMztLTBiivW/8bqe/rpXkI2NDYs4PIg+0Ff/AxP3wZ0tMDUfbfFcOpBiVxIEYhx42cbyIFtfpJu4T+wUV6D/Rm2erRR7BYk68CQIMBYTz4Y22zOaJJ3436i6iuRb6VPLIFs4q+weHKu//yWrNraVOOsokh5B0KjTxO83yGHFVvfdG4H4ybtuIukXBLs/3ksqM/RlBeJ91x8CAv34QgZ5dvk/EmO9Hf8YEDQsN4mQqLHdblt3NV7XvdUfiMI5vzPAWEXIZ1s/C3n/54A8gOkHA2KMliieLew9/ikgKIyqgLNz7jAB0klJzgFiSw2hYw/1xZvsjtqeB6KiDr95prFPYy2I9xmS7n8yiDk7flUIwrr6rwRx476FWwsVAXiRfKqyK39heR7Iqsliw0z3a0E2yp9Fu8oFqc3vdgw6b/ciyKf5Ed/1uxdqEtSPLBsQ9IGwt7YwRDk3CLyw/EoG2YJXLGidEPJzqu137rjG6+uufa6YBNmgC8ujIKjkSiqmUhip7IxmnZDsOmS9SRBUZeP3ln0mrNtb2lXglcSnx9c7O9PThiQNAuqeQZEeg2Bri6uKy3vJk7irlF79fYlnANmQrt+VDLKFL6qcoDqvvqPlWRKJydl3QDjnAMFFzgjIsz9nVcdWmSDPuhl2ihtegbrkgOywMiVBcNmZVRKkuK4cBFsIHkSF+Umj4a+v8LxZDkbzt3JBgmYZ8ztb+3R+vz8LIKp5XVZIcno0kCwm99YuqIYlQWiPiRGaObBU7/exTz+4raVomcusl/rZYFcjBRLahRKQNQuycTGlK39OBJCtvtX7HtRUAQiqJiVAdn4TX7cECe2vv85+gtPx0PyCoTMWBLihhEBCdXIgQuqUBNk0tkbfNojtauAQ68BKJ1r0V0tBdkyo+c3tc756ndxa8P7HCQahIcoWHIOcdAHhfvUb3OezkCTowwYCgadXsDGiQSM6H9IBZIejOgjizydSkk3C/G7QQWb0PQdhvD9hS5lLQdCMxrUH0fv8oTkmxIDEPDu4tZaG6iaxegThsCS8MhB8fgiCwOyGPOmGtC+DH+PD/kRJSDgcNHlagpDqgAPZgMQ5iEHgMcfl/YzjBKcgsZL4XqgOh3NBFvN4/IvrNfceBB3YpyC4fRnqOj6XSgxwZZs6pjpEQbbG67wgFU45RHKG3IPgV3rQ7xwdc4yo0CpMzHUvVL9xTDvLMbULLkh7EZyDrOvA084siN8e+llZLbD/Pavra69CGASXTEMQF9BokNxWCQRZIxByp7zk8biZIehmmAoD6oWGW8s1GF+0Dmc2rwKQawNCX94eqIHv+l8VgqBRiABka45yGxAmEiyTiEIXSXDbx3X9n0pB9Lc+CczxvJkVMvebTcwIB0jfW4HQqz0oyEYH35zbzwHB5hg5RGVmVy3ILq/lviOXGluQDQZchPsnMtkRV/YQxMZaBqRCIK8ufY8vej0oB7Ji6wsbccY0bn4pCBgWmHAgrzQR55fiHWIGiM5WnriFHOJCrJR4B2nL4RNma70GpREZBNT9i0CktSE3REVBnq3ZsFk9UfbctcM3TFz3AfIUCRoZkSibLUyMbYbmNx8EX17CgrDKHluRMJ7VEtjl3WKHmL9wlfnVgCTNb2JviYkVKxJlNpP1MT5EKQGBVeZrB5JyiCklkbKl+kEnIQmYH0BBY9neCmp56RAlubeE4oP/vumnmGtSDvIKj36+OpBE0JgUiVAO8l4itpqfloIwhftkGJ8WCSzQYY7o0ANcpSCv7pys5YCJ1YpNrNIicRcrMZaWHXrYKq5kWkhi7ih5hSDRVDeHxLoDZGj9u1bJUzO6Uw4SRDO0+FA2TN6src853JcO37VKbUAFTXBbkBPKNfxTqhyUs5hNf0prH3irvHXv2ewKgleqQNcSBL9taoI4jL2BnxqQXYcx5lTJtB3I1oZgKxq5mMR8ia31i9/218Ej5nXeU0XstiD4ijkoKs5/vpgHrm0GeehdXo6SbCu0A6GX/hFRWX1EIDubhxIO7lMZRGz0tAQhUQ8UyJIxAjWIO+pwjTgqezFILkiHxYOwNQl8p/cEgOz8peRQPzJLRDJIiXxCEHpVqX1iyQi8CFUtV31M994FkE2JwvAgMFfzIOA9xWBvvTQB7QOde3D3lqdb1gJIkeYXgJD3FCMQP//hQZhUsBwke2/lg0jW7EXvrKCEvcNXuv4lIGC01u8sdDL9+mtAOGWnL1xmQOAkCukLE5CsoZqOIJL5lfzLC5xE8SC0Ux9wKJxYCSBXvYCgSmpUIg7EfvU7OjuBOXCd+DzmV9KFqI4EEuF7bNbDjHGuKzvEbI6CoFE2v0hH7MPS+SIokLEtCSRAShYHgiqpkxTgix87W2SCVBWOA3oA2TAndoTEaosuXvaAL2wwQkCg2dopd+pP9QWyYctwwlE+3AkHIFzNLQLiG6dVJsgmaYhVWHqAo7M401X8DM5L82hNT949b7C1AEjlizRZINbIATCWlBEJWw7aopcJPAMQUDx8lUCgitCzAXGQjW1PXnmyECSsaqECHfyRUCZ6eeV6urvxg2R+YcacAbLxtuHKgVzh37Hvhg3UfcyWTHlJaRB7Qu7ag0gOcee2aB4InEWQJCLWqvkfbIG5eUYgYRV0h0OU9iAbL8DKloQ32SBSW2GM2m8AJKyXkKON161Bckqpke7BNrRmbsi5YkqmO1r3IcOErUGwMygHATykuZNZxHbDi813iRD9XkmCWIGslrFLeJIgKnA0DQeUlwgipbrQZmSAoKvHRZBxUFwPdT7s9yA8AYTtekLPvszy7Bvyag3JIrB+BFgpeKEAVhX31yQQNgFufuCvz0+C0ORBBIkd9CZXPAgqL4HwRaJXHwf4z2MgJM0TtT0C4o/r4l4pcZciyM6/mYEQOiueBqEvBJJ2IBc0eoGAQ4well4FI4Igf78DWbqb7cpIrMaZgwQbmQPNMEMQ8ObTOAhoNuyM678GhLOc8yO0piOBiPdr+avCtAUHe4u2eyIg3k3iljQ9fJkAyZpRiYDgQ4wTAOii8MlzXg8RnKLxgFlnrHK3lgyyJYcYAxCrORkgNhdeeEOF4pnuyh4FCW5qdCqCspIMEFeYF+abu5vfGAg9xEiNcj6ILxcJo/MJkLRDTIE8epB7AALC6kwQNKKfAbJxxQaFpoaqFiDSRZktQJr60tpelZgGUc3FPrViZ08NZYHA47poa+WB+FP+QochuPGscg13E8av+eMYvYFQHQnSK6MiyP6lQDa+tfqUPzVUDgITPAJCHR1QEfcPLaoUCHgD9xU6OhK9S7MY5BkmeBjE3TJAQOjdoXEQcgnoZgy+uKuOIFDZg6LQi+dwR2cCXX/MB0Hnfq8aAdnC8lMLEMn82psoF0GsJQ0IBF9JAoTc+Lvx0fZVKxBQKsR9EjIoaEFATR6LJNikURB8mbRW93oU5D5ybVAKBB1i9G0Ed6khyUcU10ksByGXr1092RpJoi8qgjz7Q4wrkpCYvIhkiHKxoQwEF/b0bhJPjmSBcHewmp+gKtKL8xZS+adIR2g32OAl29QyyBYfYnz2l1fgMp0DQdcmtbVaG1r8TgFkgJCDbmDIH42mvvjHBXsb7K0iP7Kh3eDuIPR8mDt2MdFCcb9mQcRJhyLPvgHGcp0LUseXERB8B2v0zodXGBmucPW9LNba0E5XBkcd9MdAnsEhxsnWXznGXJXwGh3ZUAXRbwBylSOQk0hiIO64rqpvGADZKgOyi4FY15CRj7SQyCYJAu502Y7ppeVxEPS4cHQrlSG20JEMEN9W2MJstRhE5efsLaxWtFwKUSbPdhpiba8cK9tabBVlp9jpoBZ+JEcibsFslZy9TCt7M7RB61ogBcPvVsCSFYxUkbIDtZeqdc+M+WWcBa00kuovH2v5m0M2MIAP585VyvxCEP69CgAkMfo3ruCTN50J11jA0a+v1JpaL3ryjTlbi0muckGEah0Ekbq4ViO0+bu2W814qEZCOB+hb+kxNX3IwSQneSDSS1QAyA695yMM1aFOgLNPmvjbm197HEW/vTkfoH/sjrfNbt7Qes9aAQj42T8feh3xA3xE1hHd33P6AIK8obzm5m3v7vWqn3yvZbnUXroriL4XKgT5QJe83EZBYFytbgnIHg4W3rwBOwF+WMvyJgQ5aLl3B7Gp7n0cBL6VXv8qAnmDA+r1n1zA1/zswZU5AxA72NUN5OgGy+/jO+uDHD+nIHtfjbnRf3hwItij1IKCHFxvsouOnL5oXSHUzcXEzsJ3FlGQN32sxmjFHuYS0BJoBUIgB183mnYCOfp2dFQgx2qORroDkDc7NIX2Ui2CPTq4TkDgrKQEcsBpd8WCfLh2dFQgH2ike1kbBgry1nhvbMNqEKxdGOQAi9KSSA7Es09ZkA/zerw4xwcdvmdATgA31BhTkDkFQWfZYiAg1hJAPvQXmeAgIIsKgzhZNIrvM4DT93yj0Nt1MAjuP0h7i/RqJJCsBW9FWFEQczu8fcY9Suhummzi0b4mCIIccEdI2luke/beBQTOAFAQp2YtQFCPLqokbrNOzwWyt4NVViT81loyW+tA7loQQaAB7wYS2VrYjRcpOzqnE1ESZM/fO4JIyo7dOAOyr0Q/QicLRE8C/g8nkEMdDHQEgeYX5fYOBDlE0bPTMaKpKBIfhQKOE0k5iOgQ0dd/g7Wdhig01qJjRCLIQaendUjmfkWb/GkxyJFcq+hCFO77JrsNZTLtQBqSkzymSCBKtQDBt0gBEEaVcSZlO7HahrUEMQ8+JX9OOnJeSer/dUXCeOy3rZKQxKoR0NztvRYg9WZCf2gmRlqAINtzC0CA17UgONV9AydC3gKQLGVnqGo33AbkiFT2w4PgG0uxszcB2F7ZQhIWSGB+i0DSsS7/YwVm10SQmQ9bxiY5MVsNcAkOMeJHWknkqJWIFYk1Rfqvx0GagNgKYI+5/C+Rke1CkFkM5Fj/nywJ2OkfaZDTw9/AP6BA3y6faayWJPVLraTVsrv7lkG0V958iCDhs/pvn/nsYI9JaltYoiJJP0KeljCOQVIMQZYZINx6x3drxEEO+McJEF/dvOVIgDSRH3Ej01Y/bjJB0NGfhE7gMt4hGqIcQaWa33iODztEEErtFa8QHMgBjvZNoxxB8Us7ehEEOGFZjyAIql/XEcjeTrLlgLzDQbIERzNsRH4ggcCwqBikCRr9Oc0ckIM/XRgRiFzGk0FAmTYLhGRRLswKCtY8yLvnjgkEvLMqG8RH3JkgOPmAxwlyQN5tQy8mEHgN+TQH5Ihus8wD2YM3utUNBTAqnRQJsEdRTZdtmwzyUAqCY3b8f2aBaEP6HlsRb9MnCIzZ9/i9GHkgqRXz/32CwJh9j4cbU3srE6SSI7IedeTN1uF1VovnpvoBicXI/YJoZ64fmo5A9AMSqX4JIDB7KgFxQPQ+815A6JGtaR7IQ5lnJyDkqvxMkKj9PcRKFCII6CS3AnnAIDdZIHFPcsDfTh7IEVzQ/tEChF6skgWiaPEtAnLasHkgJJ/9DJCDGxLvEeTDnC2+T0bxHUB8At8IxB1nnqZBlrkgR/f1fJwLZG9uFjAgB7AJ+tMRW3xIdas/ZGVPgezdV2VAwB3fsvktBzGJeVIgrc2vOVWm85VmZ3lDORVBpHpkBOQjjyMJUv+fivuVO1/R1iDAdfEgyLMvYp69VX2eA8kIUWAOakF8KZYHyYy16hrpLEsIaRDygo9QRXD76h1mf2L9wTexV9H2tK0wJ2RCVYcHwS/4uOF/AeZIB0rGgchjBf8gDmPGb+MczfBPHGSPX/ARgthnWtsHzwCJDXp4kKOtx8TLJu7ahNsoyBs+xyeA+g5eHkhk9OYf8IQ5YeLRF5Zu4yDwZCWjIv66RZNaZIH4YahKqqIcaf9GEIgrLLlf4kH28Kwrx/mAutxToOwEBAYs8N6xdwEkK0s/giFSKzce5M38n0KBjs4dTLVjZ7bNgbRw3U1wUwEk6HEKAvFTtRZXANn78+CsCuGRlin2I/4pg2a0UVHCAUDAPY7i5O8RF3miIDaWmrGdBTpkpEHuQs9+oF2EQ9MyDJToH/SIj667LIH40fellZsEohuGUoeEgNRmS4He5BQJhGmHvMsg+NJWuJ38H/xZAO1aEyBNvsH/iAE5KPueCaAiptUZPngmyAzYW+DG6ctPEyCRxYC8+3AYqojuoheAkFfHAXkAElrk6Q1Ex00H804eqMkmjsgHUQTkFgrEhsRHehasPQhjtXSuS/pqh0wQ/4wSiJ28D0DW7rdagjyEwycHcxkdesIqvbUOYCRAADnaV14YEJLAdgAhr3OceouEnjpHRw5OkkkQk7kfq95AwljLkdA9M050f95R2C4ou62l9A8SRL/irkmcjgFmuw7/BD+iOoDsVWxcIMhHYiBxgVQwbD/CLFIEKdIRNBbEKAnJECP2KC4QlCLjEMWHvxSkwGr5po+0t5LphyWJqzp431UtAOWvTwVB41EESfmRvfEAMZG4KkrJWFO4s3DYzofxGKTEs++t4Y6IBNe1Wi68lW6FxAqb35JYi55c4ERyhyqN7dYhGN+FeY0ViFUSG2wVRL+uSCHuLVr7bQtCzS18ax8EASFKST6ivdR9DMRV4286gQQOkG/toKCxIEPcu9dmyyC0P9JqhbP6H8pF0V4gtK6Yl7OrG3nCl18dQILTE1JHBNVHs6oopy/6a0HwqCK/cupaJ5CbrK11NhCtEVGOrEpjDaKtVkLZzwiSsdwb7tzfkEHirr0/EKLsuSTparxSOQ6xJxAhtc0iyeiPGCcxy5uW7eRHhPJPzjqiP4nR79h4O3ntbaTfAsBExQe5IFe6vkUeMi4Qc2KiTYjiMscDn3/0CPKW5LDBlioPGkGZiM8/+gTRjjHOYS15cRjvTEmYj5wBJLHkpk2GQGzcMCWdrg4cbUH2Y7GNlhYIPB9aMCMXrrT5zRCIazLOS0XSHG9cNI0re353nhxaDNYRhWRtQVBGkKPerghxUOjosWvHlQoE5ygtQWAjPuO8Hmr4HOAUPhxYK95XMP1tCwKHkeU7BZwM4JCgO05g2qb2lTLFAsEjA+1A6JXySe0GRxnQgXO774pnaD78QY+uILCyFBeInUtsrAK46MjW9g7lhjfISNqB0CEyrNaEC1zCIICUH8xvBmkbI94BZC/OJ6pw+JdekNAPCHTIZwA5KKgNVkOQleoLBBxz6gLywM7wHvy8BNxZ+O6NQNlbgphbwOddQR4ZEHB//hSCoMI9Nb+tQI7wEFjvIAc4U+RVhFzrQhziF4OwOoKmvKYiiApvR+kd5KYLyAFfTiiB4KCxueGlvY6wIPG+IQK5w80vAwLvGJ/JIOial0P09HQMRLRa+1h6ixNG1iGSCzynEghMrIznuW0B4l/DS0FUpAFKGNlYixwxE0FgqmvrEOUg3kDe8iCsSPaEEfWM7FE8Or8tgsAZFZW69EHSdn8RIwXZm1kSTiRmhAFru3/no1cRb5IjIHCQwPSiykGs7w37I6AFxQgEV4Lx0VEOxEeSyCG6EMDEY4eghpu/t2qSignjbZmUBaFNXuYA8YFeci6EKGgd0NxPqUhgQygLJKxp7+Hx6ikHModxC/XkPYCQBmkmSNhlgA7pPQrC3nPmQVRbkCO2d21B9uAKgjgITawYkDY6QgY2qdXilZ1ru4NLIVgQH20pnOrilXfjZ8ZCIGOhJcLOD9hN4Q7m0/OkRA/44cbasLfxIxGQiEPkO3Hu1y0Ib36ZC9vw3mrl2WMg1H3HrJb9C/q3ORBUWgknssnPugskcilrCCJFlOaR6AsF6PfOcrSNfqMgb1Lwm2gp2u8W38M/fc9b7fKROIi44i1FZ2ZNALbKKgjD9YkgsVzFfrXjshL9V4DEZzWd3mbfTvV1INEs2IH4N7WV7axPBIkur7bu3XllAvnTQFxZYV5lHLX4k0Hybqf6G0BcelHG8eeB5Iz1/x0g6SHyvwak1fof5D8I8i8rFUXTyFgUuwAAACV0RVh0ZGF0ZTpjcmVhdGUAMjAxOC0xMS0xNFQxNTo0MDoxMy0wNjowMByhwjAAAAAldEVYdGRhdGU6bW9kaWZ5ADIwMTgtMTEtMTRUMTU6NDA6MTMtMDY6MDBt/HqMAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] @@ -330,14 +337,14 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJN\nAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAFVBMVEX///+AgIA4vPKTUVDp\ngJFyEhJNv8S4mNHKAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+EEAxMKAvx4/yYAAB1WSURBVHja7V3p\nWePKEiUFKwPLGWAywGQAZAD5h/CwpO6u5dTSLZmZud/TrxljSzq9VNd66unpv3Kd/iPX/4H8bdf/\ngfxt1z8JZJrn+b8AZF6u8z8PZJ6vP9ezQPLvAZl+cNxudyT/NpAfHC+3+3V9Pv/TQObLy+3t4+Pj\n9nplU/LvAbnePpbr9sqm5FeAIHE5eqvLy9vHhoRNyXFAprOH4ygkdUI+Pt5fnx8AZHZeFsr90efU\nCfn4YGvrICCLbH/GSKb7H5+PmZJ1ZX1+f38va+twIPPlulxw2PEBNgrkvrJ+cHx//awtuklGgbAB\nLrL9FSLZUB4D5L5F7hNyn5L31zwQa93zNX/Hcd99PxJRf30qQI5YW/P1bZ2QFQjZ7T6QydikE9vZ\nPwv3ViXiWX93XXbPBwP5uu/2LJDZA3Jud781ifhgIC/jQNDzGRByRolDagSIf+IcDWTmQO4T8rlI\nRDAlfUCCE+ehMzJXSXKfEgCkZ7MHSsCePRIDeaGSRK8tT/xO4tMpUAImLrXS4nfOSK3LSxHt8ubl\nJvaBKO8fArkkzhHw+ylxjkwEyBcCsq4tc2bP4r4BkJd2st/gyY6PjMm6aXuvZZC+4bptd352hEYn\nkLcFyZfcj0/mPbMXBfKNgDhvpx4aLS2hxgMg4R3yQPRNzFvvAoINq3kenRIB5AaAnKx3mzAQ92lF\nHXrnit3hQNCM2L9Vwx++xc+ULHrEu3CjECDXg4D0/Fg9MjTvp8v19vPAm9S0G5DLoBUXSi0fSeIT\n8fc7Em370BkZNEejc+Sgq8yVofLs3yOhinIcjoIE+QCOAHLd/AHraXt+DBByGEHR5pwj2ZOlOgS+\nPjqFVsc13TexN9I2kEXXyyChhtXDtshsK2wciBLg25bKeHHM09YB37uMQ5/SE/kq8/BcVg9PxrBr\nU/KeXFndTtRNVKWAiBFYPTwqDIGf0k7b1NtN3U7UKXSOGUCqSpN6t+WMerulUC/fjxb8cUDasr9l\nVssSDfs5bLOekmidBEAmMJ8YCPPwJIDES5h/+9rre2RAoCg2gKwTspwMqbO6a9Ub4czsSOFlhoHU\nw/orLYk6xNC8SsQuvzYRv/jnGMilqk9fubXVdW36ayeQIh/uMvKmTwYIZGIensPVJ0MRd6+2L7aD\nYZGn5PcGkGZi3DfJXwCkKo13w8pzPphAHmBiDAEpyoAV1sVAmmPygUA6g1ibemb5g/4QkG6pVa91\nQj6/ZVj3TwC5nyMvfedI++3iMwVh3V8Aonw+Ayc7AWKEdY3NfiQQZHsOB6zNsG5K/JpAMsc58Ab0\na78EiBEgSRyI9jmSclggj9l4UgcLWdH49B4V5e4LC636VcV7lrtkl8sGeQMNpbFkfHw7jpFNC41c\nnDPa2DmbfdJwe4FcaaB2th5zDay25VbzuIQCK9BMfUgYVtZev09IbFaMA9kGiq/KqXoDM5s9SmYo\n37km/CzVrTSCQ9/fjIZaQMqU2M6HkhIUrK1hIPcszNvP87lHw4yGWl6UzcPzbvse/Nym9j6DS8vw\nMZnRUNOvtWqoMimVPql94xFAZsPrt+z2rPZ7WmXry89L2m9ZgNx8IPAcyQAx/LDW52bmQ3hKpIFk\nneHi7oZnnPhnc0CKFHcelVpag5GX1SGFYhVl7wjXoZeLEgSKc5t9MIJvR4/uA/i2iKGDEpiJ+A3e\naUSvcuJ5eAB3pMuuFmvC4TuiHzoRVpzJtg/I9eX1oPRRcS3nHol5x+bMDiBQFzoMiBe8R+KjEwi7\nw3rWPWJCwnQKLT36gPDMq/7A02FAtG3WBUSqG+Ox+d1AJmlB9wDZFGuitjxoPuIEF53H0QOklGm9\njnjWjgWinZUdQLYjVWVKPeTyM3XWGMk7HdMOIFXJef+FKfEzdeZLjZGUMc0DmUik96AUIEdUuJk6\n7VXamHYAaYp1NsGhvRb6uiu9q/aLHFIoRpIHQiyazhwg/L5+SJvbI2KvoxhJN5Dv7rW1qEYaSRE8\nYr+VlGkvU2crffumaQAdQDar/zsbe684oEq22mXvAkkzXqpo0QtZlL7NnUDqr7+60rImI72DOknO\n7cuXYhfXoKd2SDGXUHmVTiDf3UBmbNrD1KiJOsFWvxbKT4JOujSQNgzfPUCqe0roA3ALFNALksn0\nZAogz78DxMiXmplQOm/fXXKTP7bZM+0dErVq6ksPkJGlZexZHtPc/vCjyzH/yJaLLHGsQL5HgQxu\ndrqEKJA6qC2mqVWHCR+Y+2akBX96xK8Vn5hBKDCtOuzbI5vT9XP1HiMgaPzWJaQdhvdR+RSrg6oO\ntw4gveI3es6EVrSoUayfMy1903jSJ+7OcySYeaw6Uc8n/RUBUopn8tJk58neamngylrPAHmAz4YL\nl7itGhCyYPzqjX26FjGsUEVuPcDZZI0B+Qr0633ab50SaOrO61km6Up6gLDTIdjte+yRzWv9DqNY\ntKSdhfw69gjJ5ImA7LIQqwWB4giW0t0htTqA7LPZT16E1jKDOs6RHiCrF+Vj0IvihNESgbLwZO/Y\nIzv9Wic7GWYbeFCbP7/AqQK6VofU2utpXO8AP70QnwcbTCOFAmi/8hxxgeiKsWMIXtoASyBtbXET\n8aLskU6DZ5c33ryc0GWzEK+BhTikXv8aECOlhNvss8Q2kMb+cCA1+Bt6Ufa4Mn8BSPFrXZFfixUP\nB1rpLwG52S5ng7JKJ5bvcvcfKLWg+F2RoNNHF+aUAEyU3vJIIEWlEkpufWeDguPC1QTm1/pDQICb\nKv6VNPIXT9ZtLHR/FHVbIpkzBW0suetAIHYQYOoJYI+TgTyBGw3cxizX7X6z0cA9BzKcy2CV65a1\n8vDgKQdSXTpwqM7eBxvFqKShW5KcwWn4YCCruQKR6DXH1z7Op6vylFmND8lgoUCcLEXMpHBm/9Wa\nyLT5VpjfZdpO+oORUCAOSyGsZkFlOqK6H+Rzu+SUhwBxonYoMovZcoCfkcuyPad3EkjLfVZ55Ej4\nANe7LAoDnge42o4G0sbvBoGIdRPKVFRMnqge2A3EqRmpsXL2YTSiqJg8U8+xE4hHbwKBhFfL8CGu\nkmOUMh/IzfBNDQIhtXPV45aqeToEyCdkFRgEclNFgKkqtL1AzLT6PUC+G5AZT9LjgIDBylXqOUC2\ndD40SX8ASKdagYAwj/uhXFy5pWXlXI0A+eQfPQCIs9mHqlnAHpl/B8hLO4dlUHVE9QYbgsR3dBHC\ngUDMA3HIBgakJL8CJPKE9KvdugabA3nQ0jpeDdKn369sdhqTOQiI0kceyIqRNKyGLq0hPpACR5q6\nbx8GIX/25alsa5HQGl6jguxBKsqpEuj1ceEIGKx2qaYRlDl+HE2UdActpfaj/rRJOCBKnKClmM4s\nn/C5od8NCTnohv2CS0EfdZAsAUQ6x8X6JbN0UKXWQS7T7ddyaao0bBAnPKh27iAn9nIRYXGu9xMJ\nCiyT9Ny+sdvNdWDXJJirLLXN+qWybUooca+b60AgONVJzvE6be8lTmhk3v1pIAkHybaRyrYxvPh/\nEgjTnh3VgNv/vTSovwLEbUNCkcyXstwaNd7etXUckBna/EgMks9ILuRORfVIIG8k+6EUIAQHUz9V\n8C8AYWrUqtiG/jCQJPQXAamenimkzyUie+duzwHJuFAuzYwtQNqZYRMp/SqQrXb87H5JA0nQOP8u\nkBwhk15aCSfAr272OUdYzjb73YrNxEJ+E0iJ+kVZVLM0/krA5cuJhfzmOdKifq/uKm7nyBafYh8Y\nrznBAtEHAUkSltMWC8s6SdnnhYZ4zabdYfWGQNKE5dLlKlkr8WzOJUfzvgP3GKghEOZl89cW39s5\n9mNqHveRtncDIX7Pd7/yiUtb7lU0HdaESinHlzgIpCPqJxpVUj+v04So2exZwq5RIHknZ6uzehbR\nKa+b0hIMu2eHkO2SASISRkIg3O3sbxKWJkVDCA0I0NqKe8KoNzNh9FUrZIJltQKCUYdBIDNql7X+\nmhqL4dqalDcsBBLHmErBkkglREvLI23simpoudAD5BsCofwG9CAAQFZWb4P80WGrAKBVrl8vEL1l\nGZ09Nce1+HWTzno6/YDtxJoNjQDZyife1RbVB+Kifr6Br/YCQarmU1kfhnIQLi2D//EEVBS3RrYH\nCPrq04bDUg6iOKxTvaL0er8+JA8EBtJXIHaTgCgOO1mtL09ajW8VtjD83QEE5Gg8rX9ZFi/iMpNp\nMGf1+Do4UoHZDKtPog67BfFpILx+iAJxrJsgfAmL/cEKWNVhtxisAwgiEXs6ccNZ3YKWOWo7r5WE\ngipIEgtd1OEskPebr/9yfw3hfGAF+yp12aVDcfkNqtKxmX8aCKWPz/shRA15A2K19Fjvz8ocnS30\nBbK87rHQt0KAtO0Z8qqsTDXth8CMgU+yYlhuZ39BzG4ROneZaiDs8Ep3jMNklE8nm1O+vKuM+ouF\n9/ltywKqJOoh4QHGbPAK9799OjnN27fftai/mvEACDMbFBBRyp3lpjZnZGZEEmj51Ki/umkEhNYC\nUO6MNVeo8ZttU5JiQjaBAGoPtT5eDMb+GAgbEX4gziJ3NRlyl6TGBYgA+G6y6cxQZX1LM05QXWsD\nIhS8ar77F6QHFUC+oH9BGcjtL674lUBIssAZAUny7s74QEwAsasDIdWVdQlqNrW0TrqWxgRC/P5d\nQOx3AwxR9htsYqNY9WKzpx+KGCZ2Ajk5SiN4g5lZC4pJo2v0RIFIZrP7g2yq8ejbLPY1jeR3L7cB\nEaSE+A1GxzSs4JQwqptR5zuK6SUOxACIYy3pl+Sp6dPYhEgvcwPiqijx6LxhnRWPN3vzaTAaAko/\nY6UxHB3DHZRr6DYCA1YGh2p8PDq43+S8qPDpnsH9DxXaDDGsvqBhlRidC3DoTu3MOOtf7M+OVRwl\noambQQKOA8JJJt/6mPRYFFbYm6IDhI/Va7WM5cB5rqGwO0TuoMHrZ9+9Yam8p8+md0UOutHRMrpN\nEXmDWUT3AdnXVxVc1W/3JTWXzTCH8TXLj55ahpETexSIpvPd3so8Pw2e4impj0VhhbFrtjp7Wu20\nTpYdmq47iAI9g0Asrcdx8KpShxXH1SQGMYAcepk+CcTp3V5Zh+XmpIvI6b63a34AVen6tqYQwF4t\ni4ovD2RULSVAkK3G+hQL/vW1ZoPrNJW3PA5ZWx0qHwTE6lMs8z/Kh3liBQfIjrVlATGlWXPMUmlm\nUbr+3UBQFY1Xmf5rQNBmNxsuQ+tVx4XVOz4aiCV+basacp66yYTiCH/MZpe8xKkZqUBqyYaXTCj9\nxI8Rv9ZecPeImhGvK5diebQPxB04TKXRFr9sj9Rv23kwa54TcQk8RkUR4fdz+9zSilHptgvkKlwe\nBwBBvgTDsOLEyxEvsADC2n7MteqP5aLsxAHccIaWa0fCgc/NAdICp1kC/NirCQ2ftuRlw2XLVdxi\n8kShEUuLALk2J00KiAqGKXJvwySjTyJ/Q80Eyn2q87B+ZO8R7S7xgejwJKrFQ0GOSQR1KhDbTaRi\nutPlzRK/2mJ2gUxKNmixDHvDnoB8dGfqRCrk6sf2gdiaK+SAgFwETaRs+arxH7S44T9gy5irKONA\npjaBTW1IA7HCCnaynhIbqLRxCEjGlepEDyYk7zSDNEfC7m3VBHUCSeW0xWEQ8X4dnier6Qfd7NcY\nSJmQrw8vPhgC0RHPFccUGwp23ZyWGS4Qj3q8vZUfLLfKdVMORLOSUXtePSBXS8Nj3/KD5YJQgHye\nSTsxa0s1fb4DxDYeBBBvyQuKB4Yjwb0BnUSnEwiIeUBAGzP0NQ9IG1LR7OKSizBanRm0FHeAyIZA\n5rOcxW4UMebLqeh5T4P0VXlKGFaXZCKBI36sIkbc+dRAUoINE6NultIyDWSIeqkUMS4SnPYacsI9\naJw25hGqwcq9GQAZzlHZXpgXMdY7jAQtZRWNWNIHLC1nMK0iRiduZU+M6kXOVvT+ze4CwdSljlVi\nT0hxzBv5zfvFrwfEKGK0DGF3TF6MQz4HJD4QIyCfDcgrARKr1RKIl6KPgDRddWZZQyPUElYR4wCQ\n1b/0bVM3SyDzSuxzPp32Zg0JINR7wJZWDkir8jciDIrxjBykl60b6cfHIO1zCoi874QUhbiIW3LQ\n0dDqzqwhB4g28DgO0NY2KuLmQHgHbtsJtRMIMPAaDtg5JOYO5fSG3JNM6pnHuD7MimXTKWR3c4+K\nuDnhpKj7LdljqQpzcJlFjFvoEyRCGgkCcRH3k/w6XUv7yAu4q5DHSYxEQaube6L2mXOZUjLpUhjx\n8jqcJSaKGMkW1f3Rto9x98ROIIB8bSetZSti/FLxDeSTsJ0NfUAgZ9me+C7iYN3+grxIjvuna4/g\naPCeMLUuYizjAo8LmzapS2rlqGX6gIhCN5Lkj0z9qlApPaTrHPGiwaOXqA+rZRfnE/CHO5kOXSc7\nFZbfWSCBE5dzsEaNhi6W971P18pSy7AXjSQBKWI8E8oxrCg4KRtzh/brxB6dCYmmpOgc8zm0Vj0i\nq3o0JOyRgRmZQiCE0yW0Vu0IFU/diizEgT2SANLCCsxa7QYy5232AamVqy3awDJrFQ2qy5GG1iXd\nofvOkcyM1Itaq7j20gMCWpKxM5X3VjBnFgxBvVcSSIJy7ML8NrI9r5TdXMvBula7z0TfU6sVHYpY\ngnIsSP2jp2mJTFQFmmu/zVNbGmjxajvdcii/shKUYz6zIy9irCbNZrlwe0QywGw+fYpj2DhJNFGZ\nWJ8PvYvosJLapwWxMKy4Fs15vWp525D9ngKS7wQ3Sf4eabNzXZXyeq1z+XEbNeBjyjFB8uIOF9Or\n7+8nvCgs+ZbLCZ9UKj5SkkCKqRsV+RNLZ/nqk/wzSVDnrBUBqVQYbs4AaYnl0bSr8vMndaPqNil3\nXadgktQTYj5EbNIHYp63q4fwPeRq1JxFT/JGq9tk8/0SW8JlrMn4jVJAzHC0/N71JlK6tTe+Jk21\ntVQ5Vi0gOlcSPZsBsdhmVfoZ/hojzbgLhid9pxpFZWtJ7S4+jrFvNUkmmmuPp5PvQcSqHjkmEEV2\nZNeyCSC+rtVGMsAhgNxH5Uneggxz4w74GeeZddfhD8rFH8JYTcdFWRG+JBDBDi8MutWa+KRtgsSE\nRBGhMHrWcdEcAAlEJlb1AEnF6Hg8c9eEeEBqYlUteoBLC71Gkmsh4KDrupylpcmn0pud1ek4i4bL\n851ArM3Oj3EAxGOWTGYWgGeMln7b4hfY9vJANBdGOo0IeBqnDhuTTa51IKpGWY6KIhdGOo1oMc1u\nItd+qAhK0yqWMnAw3mIl2Eyk+XwoZSxnq9U1EMkiVYGArcwtqRKJ1WpIR2LXJGZguI5+3o6uL6HG\n83O7bBJhWK0TpLte9GWoiSK8mjHSCUTKngqEnLrX6m1gFNa4IuS0K2fw7t0YohZRW7YAgYyl3Pmw\nOdZ1sq5VPpkCEuu6+M8zyV0zgTzT5zT7FXLrnE576KDiGbF8AdOW31mSJ3wg0kHHcJHBeRvNUGvJ\nsDaOyxUiESs9AMJdpsbo2GWS4RVKrbK6z+CnjPLG3SPoxuhtSpnk10fnXg/PEY+gh7MjUiAfCSAY\nncmtAV6N/zkA4no3UVgBpUzn5XsHnbV8bz/CIj3V+nb1p/xAJO/Sozyk019Baos7XPYhrC6ta62r\n3MhkM6dEVpyaOLqSjWy1KAGk5JyluHfLyM2ZFqb96V++mxYBAbkWXc1VU7h5z6okEEPjtoFw46O7\nJ2mCzLC/HVrA/Q+B8I5uc5rDnSCJ9PGeBiMVyFsvEK6z9z8zk5bScdrsAEJ1dtgXY+81cP4PAaE6\n+wQ7y+4FkiprlEC698iJ1qKaVct7rgEdeRBIY3WJUyBGgAzU1RDraYQobIJ85nuvkZKtjlosDOTW\nnxN4ChJTppEiutmPJMdA3jiQzDOjk0SMTlZstdMt+CoEYhOruDgu3tkuydBzQEzPzeOAGEni+4CU\n2uKEmjkORBjwER3qGKt7NDwHAJk4H0q4CIb2SEoZdYBcYiByqEyO73oNlv7OSUt1VPxuVWVVLMaC\nUlbMJoFk6fAiIF9GM8dSx32rjuHo6GInu9sYY8g/D0tcYxWF2qAFCHPFop/kdK27EHnuj5lAIKqx\nN5gQrl8TP4rd+JvxmJpDnqW1EEsOA0ENPtQXqI2UsByyaQU5VmvIG4+e+VGeCY258k417StjAqUS\nPao/xnebYCZ/8D0/SVJH8DJAUqk3OTUR91bAg+c5Hwjd4mZapIC0ZCjzJadYaKzvVxxLAU8jYk1k\n93kTUW7TlmOJ3HEPtJyVTpNIA6aa7ZmGg07nHRAWB/oC08xDuNdr4Gm0WQPxt0JWQEy13YDIlBZ2\njrTnq2B01DpsJjyOphbjND/FSExtTScZtRQQerIrj/60hgy9NfNWs2BtAimS+l4FuUleMTkRcAFk\nbWLMY5NtQjrCIbK/Fw8/tH/WtCpytNosHE66CwCydTHmCedhqDMA0vQx9jJm89PuCwBR6vAKJJMx\nSoGI1nF0SBoSo2ndfiDL5E5bTx66kxn9UuaajR54LPhr9gztv4DUQukEUxJIe0cLCMu8N1ttjgAB\n3N2TJqPjfcSMi6YEGEB4ywtswI4Bge0cZQpTbo+wpks+kG1usUthCIjWtSoSDjiRxMTUdmOz84P0\nSCDJLoiJEHcT28+zeY7MO4BM7kgqe8QD4k/IlartzIo0gXTtEf8VcrpdwgMiTGSjubgE0iG1Am9S\nOiYWVomTfleLptPoU2n/FRNIdI6EJ0BHJN+95NTiqRZac8fJHvYqUn6twUsS5eOpFt2pOnQtixKW\nDuVLT2aEPR5ic1t2DR3XDu1XOyn0nGXd5D4QKW5x1z5eHNhhjyyn1KvrWxptistvog5A7NXnGYN5\nC9Frmy1uvgcG8mLUzpaCBIk9LGezz2c7w/foC2mfOCLCy78yXpR7PccfBZLpYJnxa/0AOaeW1sOA\nJNrRZDyNdyArDcHrnwKS+eFFdrizgXRaqcNAxGbP/jL2xs9z5kA8CEhPn3r5mon4SK28eTQQ0Xyo\npzzg58fsf6b2mzFSh0pYyk/vv3Qccr2X08cqOO7mwP0ZDcFcwyrK/jgQSMIkCh3SHo4iSbD9cSSQ\nqNBxDxclESVZU3MHkOhl4qCNMwZVb9jN+7oXCOFr7bZIZGKuDkek7xSK38TL1CDjrXdK1vLGd071\nMUAbw6XmKJC+xGCqPQkC+BqO650QUdU0BoQG4hMHMvPJ8BLgcVNT1KIMAqHJyGFiMK/QreUE74UU\ndMjUlOJ/DEgfpTxXQVXB+ZipKQs9xoEoz5I5diUv8bnQrbwJgTv1C151jo0BsdsugMijIGEAQAau\nSVbqDHrjrfzEWdupkiDhGCDqQD4WCIrO81qhw4DIMqdBIG8whxfmS0jujYjRNAtkYwG/7QXyCYBw\n/vz2SKZNSfE7dOkisCOBwJwiVRgoDsQ/DATuEZjlpYHMgB3laCDJG2MgmJxQ83swpXGY4cUDEjvG\n2suB4BfnGH9u35VAOCfGmAe50f5pIK6qwMcNHoiYwFMD4Q7Z2D9qABFteFWxGL6nwAh1LVxiBqhj\nqNesp0ibvZAUkLDqDc0Hw8hiRiX4hfO3ARAaoezOhCpvJNtiNCAllwTdc0thILfRPR8N3mNI5tOG\nbItF9QMpZ6+OjzghKBXk5TEjBKRpkoBeqYmVSflws9dK93sFavxk10jpIC8o5rVIzl16sDB87E4J\nfd8UEO3TnnR5tVUHNHsn+TgQyZ6aA6KjDLrg3QICec7aQ4cbnAh5NwpEGTZ2ZZZmt5NAxtz+/CgV\nUgtvdhR2V+QIZkeDmZu64ta9+ZvWxYBcjLGB+QNyUZiNACYv/nUX7CPniAPEORDxm4iv240AnEU7\nfLJ7QOw0MSMSN7FvO40AvKSpQV3LA2LHRDhhmHU5DQUcFXdQ+3WBnKxb5kKKPg+/eY3ZIz4Q+x0z\nsqX2OvlKOYQPvrJAElvyIdwdBwNJ0iIPsKn8NpCUFUw6tf36yhqOIWK0tXfer0/IoUAIB9Sjk1ge\nDKSTlesvBrLDvPi7gBynOv1pID1cgn85kD91/R/I33b9D2BqKIQldvOpAAAAJXRFWHRkYXRlOmNy\nZWF0ZQAyMDE3LTA0LTAzVDE0OjEwOjAyLTA1OjAwNSCOuQAAACV0RVh0ZGF0ZTptb2RpZnkAMjAx\nNy0wNC0wM1QxNDoxMDowMi0wNTowMER9NgUAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQBAMAAABykSv/AAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAFVBMVEX///+AgIA4vPKTUVDpgJFyEhJNv8S4mNHKAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+ILDgkoD8Jk3r8AAB1WSURBVHja7V3pWePKEiUFKwPLGWAywGQAZAD5h/CwpO6u5dTSLZmZud/TrxljSzq9VNd66unpv3Kd/iPX/4H8bdf/gfxt1z8JZJrn+b8AZF6u8z8PZJ6vP9ezQPLvAZl+cNxudyT/NpAfHC+3+3V9Pv/TQObLy+3t4+Pj9nplU/LvAbnePpbr9sqm5FeAIHE5eqvLy9vHhoRNyXFAprOH4ygkdUI+Pt5fnx8AZHZeFsr90efUCfn4YGvrICCLbH/GSKb7H5+PmZJ1ZX1+f38va+twIPPlulxw2PEBNgrkvrJ+cHx//awtuklGgbABLrL9FSLZUB4D5L5F7hNyn5L31zwQa93zNX/Hcd99PxJRf30qQI5YW/P1bZ2QFQjZ7T6QydikE9vZPwv3ViXiWX93XXbPBwP5uu/2LJDZA3Jud781ifhgIC/jQNDzGRByRolDagSIf+IcDWTmQO4T8rlIRDAlfUCCE+ehMzJXSXKfEgCkZ7MHSsCePRIDeaGSRK8tT/xO4tMpUAImLrXS4nfOSK3LSxHt8ublJvaBKO8fArkkzhHw+ylxjkwEyBcCsq4tc2bP4r4BkJd2st/gyY6PjMm6aXuvZZC+4bptd352hEYnkLcFyZfcj0/mPbMXBfKNgDhvpx4aLS2hxgMg4R3yQPRNzFvvAoINq3kenRIB5AaAnKx3mzAQ92lFHXrnit3hQNCM2L9Vwx++xc+ULHrEu3CjECDXg4D0/Fg9MjTvp8v19vPAm9S0G5DLoBUXSi0fSeIT8fc7Em370BkZNEejc+Sgq8yVofLs3yOhinIcjoIE+QCOAHLd/AHraXt+DBByGEHR5pwj2ZOlOgS+PjqFVsc13TexN9I2kEXXyyChhtXDtshsK2wciBLg25bKeHHM09YB37uMQ5/SE/kq8/BcVg9PxrBrU/KeXFndTtRNVKWAiBFYPTwqDIGf0k7b1NtN3U7UKXSOGUCqSpN6t+WMerulUC/fjxb8cUDasr9lVssSDfs5bLOekmidBEAmMJ8YCPPwJIDES5h/+9rre2RAoCg2gKwTspwMqbO6a9Ub4czsSOFlhoHUw/orLYk6xNC8SsQuvzYRv/jnGMilqk9fubXVdW36ayeQIh/uMvKmTwYIZGIensPVJ0MRd6+2L7aDYZGn5PcGkGZi3DfJXwCkKo13w8pzPphAHmBiDAEpyoAV1sVAmmPygUA6g1ibemb5g/4QkG6pVa91Qj6/ZVj3TwC5nyMvfedI++3iMwVh3V8Aonw+Ayc7AWKEdY3NfiQQZHsOB6zNsG5K/JpAMsc58Ab0a78EiBEgSRyI9jmSclggj9l4UgcLWdH49B4V5e4LC636VcV7lrtkl8sGeQMNpbFkfHw7jpFNC41cnDPa2DmbfdJwe4FcaaB2th5zDay25VbzuIQCK9BMfUgYVtZev09IbFaMA9kGiq/KqXoDM5s9SmYo37km/CzVrTSCQ9/fjIZaQMqU2M6HkhIUrK1hIPcszNvP87lHw4yGWl6UzcPzbvse/Nym9j6DS8vwMZnRUNOvtWqoMimVPql94xFAZsPrt+z2rPZ7WmXry89L2m9ZgNx8IPAcyQAx/LDW52bmQ3hKpIFkneHi7oZnnPhnc0CKFHcelVpag5GX1SGFYhVl7wjXoZeLEgSKc5t9MIJvR4/uA/i2iKGDEpiJ+A3eaUSvcuJ5eAB3pMuuFmvC4TuiHzoRVpzJtg/I9eX1oPRRcS3nHol5x+bMDiBQFzoMiBe8R+KjEwi7w3rWPWJCwnQKLT36gPDMq/7A02FAtG3WBUSqG+Ox+d1AJmlB9wDZFGuitjxoPuIEF53H0QOklGm9jnjWjgWinZUdQLYjVWVKPeTyM3XWGMk7HdMOIFXJef+FKfEzdeZLjZGUMc0DmUik96AUIEdUuJk67VXamHYAaYp1NsGhvRb6uiu9q/aLHFIoRpIHQiyazhwg/L5+SJvbI2KvoxhJN5Dv7rW1qEYaSRE8Yr+VlGkvU2crffumaQAdQDar/zsbe684oEq22mXvAkkzXqpo0QtZlL7NnUDqr7+60rImI72DOknO7cuXYhfXoKd2SDGXUHmVTiDf3UBmbNrD1KiJOsFWvxbKT4JOujSQNgzfPUCqe0roA3ALFNALksn0ZAogz78DxMiXmplQOm/fXXKTP7bZM+0dErVq6ksPkJGlZexZHtPc/vCjyzH/yJaLLHGsQL5HgQxudrqEKJA6qC2mqVWHCR+Y+2akBX96xK8Vn5hBKDCtOuzbI5vT9XP1HiMgaPzWJaQdhvdR+RSrg6oOtw4gveI3es6EVrSoUayfMy1903jSJ+7OcySYeaw6Uc8n/RUBUopn8tJk58neamngylrPAHmAz4YLl7itGhCyYPzqjX26FjGsUEVuPcDZZI0B+Qr0633ab50SaOrO61km6Up6gLDTIdjte+yRzWv9DqNYtKSdhfw69gjJ5ImA7LIQqwWB4giW0t0htTqA7LPZT16E1jKDOs6RHiCrF+Vj0IvihNESgbLwZO/YIzv9Wic7GWYbeFCbP7/AqQK6VofU2utpXO8AP70QnwcbTCOFAmi/8hxxgeiKsWMIXtoASyBtbXET8aLskU6DZ5c33ryc0GWzEK+BhTikXv8aECOlhNvss8Q2kMb+cCA1+Bt6Ufa4Mn8BSPFrXZFfixUPB1rpLwG52S5ng7JKJ5bvcvcfKLWg+F2RoNNHF+aUAEyU3vJIIEWlEkpufWeDguPC1QTm1/pDQICbKv6VNPIXT9ZtLHR/FHVbIpkzBW0suetAIHYQYOoJYI+TgTyBGw3cxizX7X6z0cA9BzKcy2CV65a18vDgKQdSXTpwqM7eBxvFqKShW5KcwWn4YCCruQKR6DXH1z7Op6vylFmND8lgoUCcLEXMpHBm/9WayLT5VpjfZdpO+oORUCAOSyGsZkFlOqK6H+Rzu+SUhwBxonYoMovZcoCfkcuyPad3EkjLfVZ55Ej4ANe7LAoDnge42o4G0sbvBoGIdRPKVFRMnqge2A3EqRmpsXL2YTSiqJg8U8+xE4hHbwKBhFfL8CGukmOUMh/IzfBNDQIhtXPV45aqeToEyCdkFRgEclNFgKkqtL1AzLT6PUC+G5AZT9LjgIDBylXqOUC2dD40SX8ASKdagYAwj/uhXFy5pWXlXI0A+eQfPQCIs9mHqlnAHpl/B8hLO4dlUHVE9QYbgsR3dBHCgUDMA3HIBgakJL8CJPKE9KvdugabA3nQ0jpeDdKn369sdhqTOQiI0kceyIqRNKyGLq0hPpACR5q6bx8GIX/25alsa5HQGl6jguxBKsqpEuj1ceEIGKx2qaYRlDl+HE2UdActpfaj/rRJOCBKnKClmM4sn/C5od8NCTnohv2CS0EfdZAsAUQ6x8X6JbN0UKXWQS7T7ddyaao0bBAnPKh27iAn9nIRYXGu9xMJCiyT9Ny+sdvNdWDXJJirLLXN+qWybUooca+b60AgONVJzvE6be8lTmhk3v1pIAkHybaRyrYxvPh/EgjTnh3VgNv/vTSovwLEbUNCkcyXstwaNd7etXUckBna/EgMks9ILuRORfVIIG8k+6EUIAQHUz9V8C8AYWrUqtiG/jCQJPQXAamenimkzyUie+duzwHJuFAuzYwtQNqZYRMp/SqQrXb87H5JA0nQOP8ukBwhk15aCSfAr272OUdYzjb73YrNxEJ+E0iJ+kVZVLM0/krA5cuJhfzmOdKifq/uKm7nyBafYh8YrznBAtEHAUkSltMWC8s6SdnnhYZ4zabdYfWGQNKE5dLlKlkr8WzOJUfzvgP3GKghEOZl89cW39s59mNqHveRtncDIX7Pd7/yiUtb7lU0HdaESinHlzgIpCPqJxpVUj+v04So2exZwq5RIHknZ6uzehbRKa+b0hIMu2eHkO2SASISRkIg3O3sbxKWJkVDCA0I0NqKe8KoNzNh9FUrZIJltQKCUYdBIDNql7X+mhqL4dqalDcsBBLHmErBkkglREvLI23simpoudAD5BsCofwG9CAAQFZWb4P80WGrAKBVrl8vEL1lGZ09Nce1+HWTzno6/YDtxJoNjQDZyife1RbVB+Kifr6Br/YCQarmU1kfhnIQLi2D//EEVBS3RrYHCPrq04bDUg6iOKxTvaL0er8+JA8EBtJXIHaTgCgOO1mtL09ajW8VtjD83QEE5Gg8rX9ZFi/iMpNpMGf1+Do4UoHZDKtPog67BfFpILx+iAJxrJsgfAmL/cEKWNVhtxisAwgiEXs6ccNZ3YKWOWo7r5WEgipIEgtd1OEskPebr/9yfw3hfGAF+yp12aVDcfkNqtKxmX8aCKWPz/shRA15A2K19Fjvz8ocnS30BbK87rHQt0KAtO0Z8qqsTDXth8CMgU+yYlhuZ39BzG4ROneZaiDs8Ep3jMNklE8nm1O+vKuM+ouF9/ltywKqJOoh4QHGbPAK9799OjnN27fftai/mvEACDMbFBBRyp3lpjZnZGZEEmj51Ki/umkEhNYCUO6MNVeo8ZttU5JiQjaBAGoPtT5eDMb+GAgbEX4gziJ3NRlyl6TGBYgA+G6y6cxQZX1LM05QXWsDIhS8ar77F6QHFUC+oH9BGcjtL674lUBIssAZAUny7s74QEwAsasDIdWVdQlqNrW0TrqWxgRC/P5dQOx3AwxR9htsYqNY9WKzpx+KGCZ2Ajk5SiN4g5lZC4pJo2v0RIFIZrP7g2yq8ejbLPY1jeR3L7cBEaSE+A1GxzSs4JQwqptR5zuK6SUOxACIYy3pl+Sp6dPYhEgvcwPiqijx6LxhnRWPN3vzaTAaAko/Y6UxHB3DHZRr6DYCA1YGh2p8PDq43+S8qPDpnsH9DxXaDDGsvqBhlRidC3DoTu3MOOtf7M+OVRwloambQQKOA8JJJt/6mPRYFFbYm6IDhI/Va7WM5cB5rqGwO0TuoMHrZ9+9Yam8p8+md0UOutHRMrpNEXmDWUT3AdnXVxVc1W/3JTWXzTCH8TXLj55ahpETexSIpvPd3so8Pw2e4impj0VhhbFrtjp7Wu20TpYdmq47iAI9g0Asrcdx8KpShxXH1SQGMYAcepk+CcTp3V5Zh+XmpIvI6b63a34AVen6tqYQwF4ti4ovD2RULSVAkK3G+hQL/vW1ZoPrNJW3PA5ZWx0qHwTE6lMs8z/Kh3liBQfIjrVlATGlWXPMUmlmUbr+3UBQFY1Xmf5rQNBmNxsuQ+tVx4XVOz4aiCV+basacp66yYTiCH/MZpe8xKkZqUBqyYaXTCj9xI8Rv9ZecPeImhGvK5diebQPxB04TKXRFr9sj9Rv23kwa54TcQk8RkUR4fdz+9zSilHptgvkKlweBwBBvgTDsOLEyxEvsADC2n7MteqP5aLsxAHccIaWa0fCgc/NAdICp1kC/NirCQ2ftuRlw2XLVdxi8kShEUuLALk2J00KiAqGKXJvwySjTyJ/Q80Eyn2q87B+ZO8R7S7xgejwJKrFQ0GOSQR1KhDbTaRiutPlzRK/2mJ2gUxKNmixDHvDnoB8dGfqRCrk6sf2gdiaK+SAgFwETaRs+arxH7S44T9gy5irKONApjaBTW1IA7HCCnaynhIbqLRxCEjGlepEDyYk7zSDNEfC7m3VBHUCSeW0xWEQ8X4dnier6Qfd7NcYSJmQrw8vPhgC0RHPFccUGwp23ZyWGS4Qj3q8vZUfLLfKdVMORLOSUXtePSBXS8Nj3/KD5YJQgHyeSTsxa0s1fb4DxDYeBBBvyQuKB4Yjwb0BnUSnEwiIeUBAGzP0NQ9IG1LR7OKSizBanRm0FHeAyIZA5rOcxW4UMebLqeh5T4P0VXlKGFaXZCKBI36sIkbc+dRAUoINE6NultIyDWSIeqkUMS4SnPYacsI9aJw25hGqwcq9GQAZzlHZXpgXMdY7jAQtZRWNWNIHLC1nMK0iRiduZU+M6kXOVvT+ze4CwdSljlViT0hxzBv5zfvFrwfEKGK0DGF3TF6MQz4HJD4QIyCfDcgrARKr1RKIl6KPgDRddWZZQyPUElYR4wCQ1b/0bVM3SyDzSuxzPp32Zg0JINR7wJZWDkir8jciDIrxjBykl60b6cfHIO1zCoi874QUhbiIW3LQ0dDqzqwhB4g28DgO0NY2KuLmQHgHbtsJtRMIMPAaDtg5JOYO5fSG3JNM6pnHuD7MimXTKWR3c4+KuDnhpKj7LdljqQpzcJlFjFvoEyRCGgkCcRH3k/w6XUv7yAu4q5DHSYxEQaube6L2mXOZUjLpUhjx8jqcJSaKGMkW1f3Rto9x98ROIIB8bSetZSti/FLxDeSTsJ0NfUAgZ9me+C7iYN3+grxIjvuna4/gaPCeMLUuYizjAo8LmzapS2rlqGX6gIhCN5Lkj0z9qlApPaTrHPGiwaOXqA+rZRfnE/CHO5kOXSc7FZbfWSCBE5dzsEaNhi6W971P18pSy7AXjSQBKWI8E8oxrCg4KRtzh/brxB6dCYmmpOgc8zm0Vj0iq3o0JOyRgRmZQiCE0yW0Vu0IFU/diizEgT2SANLCCsxa7QYy5232AamVqy3awDJrFQ2qy5GG1iXdofvOkcyM1Itaq7j20gMCWpKxM5X3VjBnFgxBvVcSSIJy7ML8NrI9r5TdXMvBula7z0TfU6sVHYpYgnIsSP2jp2mJTFQFmmu/zVNbGmjxajvdcii/shKUYz6zIy9irCbNZrlwe0QywGw+fYpj2DhJNFGZWJ8PvYvosJLapwWxMKy4Fs15vWp525D9ngKS7wQ3Sf4eabNzXZXyeq1z+XEbNeBjyjFB8uIOF9Or7+8nvCgs+ZbLCZ9UKj5SkkCKqRsV+RNLZ/nqk/wzSVDnrBUBqVQYbs4AaYnl0bSr8vMndaPqNil3XadgktQTYj5EbNIHYp63q4fwPeRq1JxFT/JGq9tk8/0SW8JlrMn4jVJAzHC0/N71JlK6tTe+Jk21tVQ5Vi0gOlcSPZsBsdhmVfoZ/hojzbgLhid9pxpFZWtJ7S4+jrFvNUkmmmuPp5PvQcSqHjkmEEV2ZNeyCSC+rtVGMsAhgNxH5Uneggxz4w74GeeZddfhD8rFH8JYTcdFWRG+JBDBDi8MutWa+KRtgsSERBGhMHrWcdEcAAlEJlb1AEnF6Hg8c9eEeEBqYlUteoBLC71Gkmsh4KDrupylpcmn0pud1ek4i4bL851ArM3Oj3EAxGOWTGYWgGeMln7b4hfY9vJANBdGOo0IeBqnDhuTTa51IKpGWY6KIhdGOo1oMc1uItd+qAhK0yqWMnAw3mIl2Eyk+XwoZSxnq9U1EMkiVYGArcwtqRKJ1WpIR2LXJGZguI5+3o6uL6HG83O7bBJhWK0TpLte9GWoiSK8mjHSCUTKngqEnLrX6m1gFNa4IuS0K2fw7t0YohZRW7YAgYyl3PmwOdZ1sq5VPpkCEuu6+M8zyV0zgTzT5zT7FXLrnE576KDiGbF8AdOW31mSJ3wg0kHHcJHBeRvNUGvJsDaOyxUiESs9AMJdpsbo2GWS4RVKrbK6z+CnjPLG3SPoxuhtSpnk10fnXg/PEY+gh7MjUiAfCSAYncmtAV6N/zkA4no3UVgBpUzn5XsHnbV8bz/CIj3V+nb1p/xAJO/Sozyk019Baos7XPYhrC6ta62r3MhkM6dEVpyaOLqSjWy1KAGk5JyluHfLyM2ZFqb96V++mxYBAbkWXc1VU7h5z6okEEPjtoFw46O7J2mCzLC/HVrA/Q+B8I5uc5rDnSCJ9PGeBiMVyFsvEK6z9z8zk5bScdrsAEJ1dtgXY+81cP4PAaE6+wQ7y+4FkiprlEC698iJ1qKaVct7rgEdeRBIY3WJUyBGgAzU1RDraYQobIJ85nuvkZKtjlosDOTWnxN4ChJTppEiutmPJMdA3jiQzDOjk0SMTlZstdMt+CoEYhOruDgu3tkuydBzQEzPzeOAGEni+4CU2uKEmjkORBjwER3qGKt7NDwHAJk4H0q4CIb2SEoZdYBcYiByqEyO73oNlv7OSUt1VPxuVWVVLMaCUlbMJoFk6fAiIF9GM8dSx32rjuHo6GInu9sYY8g/D0tcYxWF2qAFCHPFop/kdK27EHnuj5lAIKqxN5gQrl8TP4rd+JvxmJpDnqW1EEsOA0ENPtQXqI2UsByyaQU5VmvIG4+e+VGeCY258k417StjAqUSPao/xnebYCZ/8D0/SVJH8DJAUqk3OTUR91bAg+c5Hwjd4mZapIC0ZCjzJadYaKzvVxxLAU8jYk1k93kTUW7TlmOJ3HEPtJyVTpNIA6aa7ZmGg07nHRAWB/oC08xDuNdr4Gm0WQPxt0JWQEy13YDIlBZ2jrTnq2B01DpsJjyOphbjND/FSExtTScZtRQQerIrj/60hgy9NfNWs2BtAimS+l4FuUleMTkRcAFkbWLMY5NtQjrCIbK/Fw8/tH/WtCpytNosHE66CwCydTHmCedhqDMA0vQx9jJm89PuCwBR6vAKJJMxSoGI1nF0SBoSo2ndfiDL5E5bTx66kxn9UuaajR54LPhr9gztv4DUQukEUxJIe0cLCMu8N1ttjgAB3N2TJqPjfcSMi6YEGEB4ywtswI4Bge0cZQpTbo+wpks+kG1usUthCIjWtSoSDjiRxMTUdmOz84P0SCDJLoiJEHcT28+zeY7MO4BM7kgqe8QD4k/IlartzIo0gXTtEf8VcrpdwgMiTGSjubgE0iG1Am9SOiYWVomTfleLptPoU2n/FRNIdI6EJ0BHJN+95NTiqRZac8fJHvYqUn6twUsS5eOpFt2pOnQtixKWDuVLT2aEPR5ic1t2DR3XDu1XOyn0nGXd5D4QKW5x1z5eHNhhjyyn1KvrWxptistvog5A7NXnGYN5C9Frmy1uvgcG8mLUzpaCBIk9LGezz2c7w/foC2mfOCLCy78yXpR7PccfBZLpYJnxa/0AOaeW1sOAJNrRZDyNdyArDcHrnwKS+eFFdrizgXRaqcNAxGbP/jL2xs9z5kA8CEhPn3r5mon4SK28eTQQ0Xyopzzg58fsf6b2mzFSh0pYyk/vv3Qccr2X08cqOO7mwP0ZDcFcwyrK/jgQSMIkCh3SHo4iSbD9cSSQqNBxDxclESVZU3MHkOhl4qCNMwZVb9jN+7oXCOFr7bZIZGKuDkek7xSK38TL1CDjrXdK1vLGd071MUAbw6XmKJC+xGCqPQkC+BqO650QUdU0BoQG4hMHMvPJ8BLgcVNT1KIMAqHJyGFiMK/QreUE74UUdMjUlOJ/DEgfpTxXQVXB+ZipKQs9xoEoz5I5diUv8bnQrbwJgTv1C151jo0BsdsugMijIGEAQAauSVbqDHrjrfzEWdupkiDhGCDqQD4WCIrO81qhw4DIMqdBIG8whxfmS0jujYjRNAtkYwG/7QXyCYBw/vz2SKZNSfE7dOkisCOBwJwiVRgoDsQ/DATuEZjlpYHMgB3laCDJG2MgmJxQ83swpXGY4cUDEjvG2suB4BfnGH9u35VAOCfGmAe50f5pIK6qwMcNHoiYwFMD4Q7Z2D9qABFteFWxGL6nwAh1LVxiBqhjqNesp0ibvZAUkLDqDc0Hw8hiRiX4hfO3ARAaoezOhCpvJNtiNCAllwTdc0thILfRPR8N3mNI5tOGbItF9QMpZ6+OjzghKBXk5TEjBKRpkoBeqYmVSflws9dK93sFavxk10jpIC8o5rVIzl16sDB87E4Jfd8UEO3TnnR5tVUHNHsn+TgQyZ6aA6KjDLrg3QICec7aQ4cbnAh5NwpEGTZ2ZZZmt5NAxtz+/CgVUgtvdhR2V+QIZkeDmZu64ta9+ZvWxYBcjLGB+QNyUZiNACYv/nUX7CPniAPEORDxm4iv240AnEU7fLJ7QOw0MSMSN7FvO40AvKSpQV3LA2LHRDhhmHU5DQUcFXdQ+3WBnKxb5kKKPg+/eY3ZIz4Q+x0zsqX2OvlKOYQPvrJAElvyIdwdBwNJ0iIPsKn8NpCUFUw6tf36yhqOIWK0tXfer0/IoUAIB9Sjk1geDKSTlesvBrLDvPi7gBynOv1pID1cgn85kD91/R/I33b9D2BqKIQldvOpAAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDE4LTExLTE0VDE1OjQwOjE1LTA2OjAwf3H3CgAAACV0RVh0ZGF0ZTptb2RpZnkAMjAxOC0xMS0xNFQxNTo0MDoxNS0wNjowMA4sT7YAAAAASUVORK5CYII=\n", "text/plain": [ "" ] @@ -356,7 +363,7 @@ "metadata": { "anaconda-cloud": {}, "kernelspec": { - "display_name": "Python [default]", + "display_name": "Python 3", "language": "python", "name": "python3" }, @@ -370,9 +377,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.6.5" } }, "nbformat": 4, - "nbformat_minor": 0 + "nbformat_minor": 1 } From fc5bbd6a9bcb138c48e7114e64325b29119ba4e2 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 15 Nov 2018 10:46:18 -0600 Subject: [PATCH 7/9] Fix TRISO regression test --- tests/regression_tests/triso/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/triso/results_true.dat b/tests/regression_tests/triso/results_true.dat index 2bbecb4a9..dfc700c05 100644 --- a/tests/regression_tests/triso/results_true.dat +++ b/tests/regression_tests/triso/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.683227E+00 7.383566E-02 +1.683226E+00 7.383559E-02 From 9b00c61a03fa4cdc76870ab2a34d94f2e6329ba7 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 19 Nov 2018 21:55:53 -0600 Subject: [PATCH 8/9] Address @paulromano comments on #1121 --- openmc/model/triso.py | 593 ++++++++++++++++++------------------------ 1 file changed, 251 insertions(+), 342 deletions(-) diff --git a/openmc/model/triso.py b/openmc/model/triso.py index a99f136bb..e04d7f589 100644 --- a/openmc/model/triso.py +++ b/openmc/model/triso.py @@ -8,7 +8,6 @@ from collections.abc import Iterable from heapq import heappush, heappop from math import pi, sin, cos, floor, log10, sqrt from numbers import Real -from operator import attrgetter from random import uniform, gauss import numpy as np @@ -158,10 +157,7 @@ class _Container(metaclass=ABCMeta): @center.setter def center(self, center): - if np.asarray(center).size != 3: - raise ValueError('Unable to set domain center to {} since it must ' - 'be of length 3'.format(center)) - self._center = [float(x) for x in center] + self._center = center def mesh_cell(self, p): """Calculate the index of the cell in a mesh overlaid on the domain in @@ -200,6 +196,20 @@ class _Container(metaclass=ABCMeta): for i in range(3)] return list(itertools.product(*({int(x) for x in y} for y in r))) + @abstractmethod + def from_region(self, region, sphere_radius): + """Create a container to pack spheres in based on a region. + + Parameters + ---------- + region : openmc.Region + Region to create container from. + sphere_radius : float + Outer radius of spheres. + + """ + pass + @abstractmethod def random_point(self): """Generate Cartesian coordinates of center of a sphere that is @@ -268,7 +278,7 @@ class _RectangularPrism(_Container): Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum distance from center in x-, y-, or z-direction where sphere + Minimum and maximum distance in x-, y-, and z-direction where sphere center can be placed. volume : float Volume of the container. @@ -296,9 +306,11 @@ class _RectangularPrism(_Container): @property def limits(self): if self._limits is None: - self._limits = [self.width/2 - self.sphere_radius, - self.depth/2 - self.sphere_radius, - self.height/2 - self.sphere_radius] + c = self.center + r = self.sphere_radius + x, y, z = self.width/2, self.depth/2, self.height/2 + self._limits = [[c[0] - x + r, c[1] - y + r, c[2] - z + r], + [c[0] + x - r, c[1] + y - r, c[2] + z - r]] return self._limits @property @@ -335,11 +347,50 @@ class _RectangularPrism(_Container): def limits(self, limits): self._limits = limits + @classmethod + def from_region(self, region, sphere_radius): + check_type('region', region, openmc.Region) + + # Assume the simplest case where the prism volume is the intersection + # of the half-spaces of six planes + if not isinstance(region, openmc.Intersection): + raise ValueError + + if any(not isinstance(node, openmc.Halfspace) for node in region): + raise ValueError + + if len(region) != 6: + raise ValueError + + # Sort half-spaces by surface type + px1, px2, py1, py2, pz1, pz2 = sorted(region, key=lambda x: x.surface.type) + + # Make sure the region consists of the correct surfaces + if (not isinstance(px1.surface, openmc.XPlane) or + not isinstance(px2.surface, openmc.XPlane) or + not isinstance(py1.surface, openmc.YPlane) or + not isinstance(py2.surface, openmc.YPlane) or + not isinstance(pz1.surface, openmc.ZPlane) or + not isinstance(pz2.surface, openmc.ZPlane)): + raise ValueError + + # Make sure the half-spaces are on the correct side of the surfaces + ll, ur = region.bounding_box + if any(x in ll or x in ur for x in (-np.inf, np.inf)): + raise ValueError + + # Calculate the parameters for the container + width, depth, height = ur - ll + center = ll + [width/2, depth/2, height/2] + + # The region is the volume of a rectangular prism, so create container + return _RectangularPrism(width, depth, height, sphere_radius, center) + def random_point(self): - x_max, y_max, z_max = self.limits - return [uniform(-x_max, x_max), - uniform(-y_max, y_max), - uniform(-z_max, z_max)] + ll, ul = self.limits + return [uniform(ll[0], ul[0]), + uniform(ll[1], ul[1]), + uniform(ll[2], ul[2])] def repel_spheres(self, p, q, d, d_new): # Moving each sphere distance 's' away from the other along the line @@ -354,8 +405,8 @@ class _RectangularPrism(_Container): # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - p[:] = np.clip(p, [-x for x in self.limits], self.limits) - q[:] = np.clip(q, [-x for x in self.limits], self.limits) + p[:] = np.clip(p, self.limits[0], self.limits[1]) + q[:] = np.clip(q, self.limits[0], self.limits[1]) class _Cylinder(_Container): @@ -396,8 +447,8 @@ class _Cylinder(_Container): Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum radial distance and maximum distance from center in z-direction - where sphere center can be placed. + Maximum radial distance and minimum and maximum distance in the + direction parallel to the axis where sphere center can be placed. volume : float Volume of the container. @@ -436,8 +487,10 @@ class _Cylinder(_Container): @property def limits(self): if self._limits is None: - self._limits = [self.radius - self.sphere_radius, - self.length/2 - self.sphere_radius] + z0 = self.center[self.shift[2]] + z = self.length/2 + r = self.sphere_radius + self._limits = [[z0 - z + r], [z0 + z - r, self.radius - r]] return self._limits @property @@ -476,15 +529,72 @@ class _Cylinder(_Container): def limits(self, limits): self._limits = limits + @classmethod + def from_region(self, region, sphere_radius): + check_type('region', region, openmc.Region) + + # Assume the simplest case where the cylinder volume is the + # intersection of the half-spaces of a cylinder and two planes + if not isinstance(region, openmc.Intersection): + raise ValueError + + if any(not isinstance(node, openmc.Halfspace) for node in region): + raise ValueError + + if len(region) != 3: + raise ValueError + + # Identify the axis that the cylinder lies along + axis = region[0].surface.type[0] + + # Make sure the region is composed of a cylinder and two planes on the + # same axis + count = Counter(node.surface.type for node in region) + if count[axis + '-cylinder'] != 1 or count[axis + '-plane'] != 2: + raise ValueError + + # Sort the half-spaces by surface type + cyl, p1, p2 = sorted(region, key=lambda x: x.surface.type) + + # Calculate the parameters for a cylinder along the x-axis + if axis == 'x': + if p1.surface.x0 > p2.surface.x0: + p1, p2 = p2, p1 + length = p2.surface.x0 - p1.surface.x0 + center = (p1.surface.x0 + length/2, cyl.surface.y0, cyl.surface.z0) + + # Calculate the parameters for a cylinder along the y-axis + elif axis == 'y': + if p1.surface.y0 > p2.surface.y0: + p1, p2 = p2, p1 + length = p2.surface.y0 - p1.surface.y0 + center = (cyl.surface.x0, p1.surface.y0 + length/2, cyl.surface.z0) + + # Calculate the parameters for a cylinder along the z-axis + else: + if p1.surface.z0 > p2.surface.z0: + p1, p2 = p2, p1 + length = p2.surface.z0 - p1.surface.z0 + center = (cyl.surface.x0, cyl.surface.y0, p1.surface.z0 + length/2) + + # Make sure the half-spaces are on the correct side of the surfaces + if cyl.side != '-' or p1.side != '+' or p2.side != '-': + raise ValueError + + radius = cyl.surface.r + + # The region is the volume of a cylinder, so create container + return _Cylinder(length, radius, axis, sphere_radius, center) + def random_point(self): - r_max, z_max = self.limits - r = sqrt(uniform(0, r_max**2)) + ll, ul = self.limits + r = sqrt(uniform(0, ul[1]**2)) t = uniform(0, 2*pi) i, j, k = self.shift p = [None]*3 - p[i] = r*cos(t) - p[j] = r*sin(t) - p[k] = uniform(-z_max, z_max) + p[i] = r*cos(t) + self.center[i] + p[j] = r*sin(t) + self.center[j] + p[k] = uniform(ll[0], ul[0]) return p def repel_spheres(self, p, q, d, d_new): @@ -500,119 +610,24 @@ class _Cylinder(_Container): # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - r_max, z_max = self.limits + ll, ul = self.limits + c = self.center i, j, k = self.shift - r = sqrt(p[i]**2 + p[j]**2) - if r > r_max: - p[i] *= r_max/r - p[j] *= r_max/r - p[k] = np.clip(p[k], -z_max, z_max) + r = sqrt((p[i] - c[i])**2 + (p[j] - c[j])**2) + if r > ul[1]: + p[i] = (p[i] - c[i])*ul[1]/r + c[i] + p[j] = (p[j] - c[j])*ul[1]/r + c[j] + p[k] = np.clip(p[k], ll[0], ul[0]) - r = sqrt(q[i]**2 + q[j]**2) - if r > r_max: - q[i] *= r_max/r - q[j] *= r_max/r - q[k] = np.clip(q[k], -z_max, z_max) + r = sqrt((q[i] - c[i])**2 + (q[j] - c[j])**2) + if r > ul[1]: + q[i] = (q[i] - c[i])*ul[1]/r + c[i] + q[j] = (q[j] - c[j])*ul[1]/r + c[j] + q[k] = np.clip(q[k], ll[0], ul[0]) -class _Sphere(_Container): - """Spherical container in which to pack spheres. - - Parameters - ---------- - radius : float - Radius of the spherical container. - center : Iterable of float - Cartesian coordinates of the center of the container. Default is - (0., 0., 0.) - - Attributes - ---------- - radius : float - Radius of the spherical container. - sphere_radius : float - Radius of spheres to be packed in container. - center : list of float - Cartesian coordinates of the center of the container. Default is - (0., 0., 0.) - cell_length : list of float - Length in x-, y-, and z- directions of each cell in mesh overlaid on - domain. - limits : list of float - Maximum radial distance where sphere center can be placed. - volume : float - Volume of the container. - - """ - - def __init__(self, radius, sphere_radius, center=(0., 0., 0.)): - super().__init__(sphere_radius, center) - self.radius = radius - - @property - def radius(self): - return self._radius - - @property - def limits(self): - if self._limits is None: - self._limits = [self.radius - self.sphere_radius] - return self._limits - - @property - def cell_length(self): - if self._cell_length is None: - mesh_length = 3*[2*self.radius] - self._cell_length = [x/int(x/(4*self.sphere_radius)) - for x in mesh_length] - return self._cell_length - - @property - def volume(self): - return 4/3*pi*self.radius**3 - - @radius.setter - def radius(self, radius): - self._radius = float(radius) - self._limits = None - self._cell_length = None - - @limits.setter - def limits(self, limits): - self._limits = limits - - def random_point(self): - r_max = self.limits[0] - x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) - r = (uniform(0, r_max**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2)) - return [r*s for s in x] - - def repel_spheres(self, p, q, d, d_new): - # Moving each sphere distance 's' away from the other along the line - # joining the sphere centers will ensure their final distance is - # equal to the outer diameter - s = (d_new - d)/2 - - v = (p - q)/d - p += s*v - q -= s*v - - # Enforce the rigid boundary by moving each sphere back along the - # surface normal until it is completely within the container if it - # overlaps the surface - r_max = self.limits[0] - - r = sqrt(p[0]**2 + p[1]**2 + p[2]**2) - if r > r_max: - p *= r_max/r - - r = sqrt(q[0]**2 + q[1]**2 + q[2]**2) - if r > r_max: - q *= r_max/r - - -class _SphericalShell(_Sphere): +class _SphericalShell(_Container): """Spherical shell container in which to pack spheres. Parameters @@ -640,8 +655,7 @@ class _SphericalShell(_Sphere): Length in x-, y-, and z- directions of each cell in mesh overlaid on domain. limits : list of float - Maximum radial distance and minimum radial distance where sphere - center can be placed. + Minimum and maximum radial distance where sphere center can be placed. volume : float Volume of the container. @@ -649,9 +663,14 @@ class _SphericalShell(_Sphere): def __init__(self, radius, inner_radius, sphere_radius, center=(0., 0., 0.)): - super().__init__(radius, sphere_radius, center) + super().__init__(sphere_radius, center) + self.radius = radius self.inner_radius = inner_radius + @property + def radius(self): + return self._radius + @property def inner_radius(self): return self._inner_radius @@ -659,14 +678,32 @@ class _SphericalShell(_Sphere): @property def limits(self): if self._limits is None: - self._limits = [self.radius - self.sphere_radius, - self.inner_radius + self.sphere_radius] + r_max = self.radius - self.sphere_radius + if self.inner_radius == 0: + r_min = 0 + else: + r_min = self.inner_radius + self.sphere_radius + self._limits = [[r_min], [r_max]] return self._limits + @property + def cell_length(self): + if self._cell_length is None: + mesh_length = 3*[2*self.radius] + self._cell_length = [x/int(x/(4*self.sphere_radius)) + for x in mesh_length] + return self._cell_length + @property def volume(self): return 4/3*pi*(self.radius**3 - self.inner_radius**3) + @radius.setter + def radius(self, radius): + self._radius = float(radius) + self._limits = None + self._cell_length = None + @inner_radius.setter def inner_radius(self, inner_radius): self._inner_radius = float(inner_radius) @@ -676,11 +713,58 @@ class _SphericalShell(_Sphere): def limits(self, limits): self._limits = limits + @classmethod + def from_region(self, region, sphere_radius): + check_type('region', region, openmc.Region) + + # First check if the region is the volume inside a sphere. Assume the + # simplest case where the sphere volume is the negative half-space of a + # sphere. + if (isinstance(region, openmc.Halfspace) + and isinstance(region.surface, openmc.Sphere) + and region.side == '-'): + + # The region is the volume of a sphere, so create container + radius = region.surface.r + center = (region.surface.x0, region.surface.y0, region.surface.z0) + + return _SphericalShell(radius, 0., sphere_radius, center) + + # Next check for a spherical shell volume. Assume the simplest case + # where the spherical shell volume is the intersection of the + # half-spaces of two spheres. + if not isinstance(region, openmc.Intersection): + raise ValueError + + if any(not isinstance(node, openmc.Halfspace) for node in region): + raise ValueError + + if len(region) != 2: + raise ValueError + + if any(not isinstance(node.surface, openmc.Sphere) for node in region): + raise ValueError + + s1, s2 = sorted(region, key=lambda x: x.surface.r) + radius = s2.surface.r + inner_radius = s1.surface.r + center = (s1.surface.x0, s1.surface.y0, s1.surface.z0) + + if center != (s2.surface.x0, s2.surface.y0, s2.surface.z0): + raise ValueError + + if s1.side != '+' or s2.side != '-': + raise ValueError + + # The region is the volume of a spherical shell, so create container + return _SphericalShell(radius, inner_radius, sphere_radius, center) + def random_point(self): - r_max, r_min = self.limits - x = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) - r = (uniform(r_min**3, r_max**3)**(1/3)/sqrt(x[0]**2 + x[1]**2 + x[2]**2)) - return [r*s for s in x] + c = self.center + ll, ul = self.limits + x, y, z = (gauss(0, 1), gauss(0, 1), gauss(0, 1)) + r = (uniform(ll[0]**3, ul[0]**3)**(1/3)/sqrt(x**2 + y**2 + z**2)) + return [r*x + c[0], r*y + c[1], r*z + c[2]] def repel_spheres(self, p, q, d, d_new): # Moving each sphere distance 's' away from the other along the line @@ -695,19 +779,20 @@ class _SphericalShell(_Sphere): # Enforce the rigid boundary by moving each sphere back along the # surface normal until it is completely within the container if it # overlaps the surface - r_max, r_min = self.limits + c = self.center + ll, ul = self.limits - r = sqrt(p[0]**2 + p[1]**2 + p[2]**2) - if r > r_max: - p *= r_max/r - elif r < r_min: - p *= r_min/r + r = sqrt((p[0] - c[0])**2 + (p[1] - c[1])**2 + (p[2] - c[2])**2) + if r > ul[0]: + p[:] = (p - c)*ul[0]/r + c + elif r < ll[0]: + p[:] = (p - c)*ll[0]/r + c - r = sqrt(q[0]**2 + q[1]**2 + q[2]**2) - if r > r_max: - q *= r_max/r - elif r < r_min: - q *= r_min/r + r = sqrt((q[0] - c[0])**2 + (q[1] - c[1])**2 + (q[2] - c[2])**2) + if r > ul[0]: + q[:] = (q - c)*ul[0]/r + c + elif r < ll[0]: + q[:] = (q - c)*ll[0]/r + c def create_triso_lattice(trisos, lower_left, pitch, shape, background): @@ -783,193 +868,6 @@ def create_triso_lattice(trisos, lower_left, pitch, shape, background): return lattice -def _create_container(region, sphere_radius): - """Create a container to pack spheres in using the region to determine the - container shape. - - Parameters - ---------- - region : openmc.Region - Container in which the spheres are packed. Supported shapes are - rectangular_prism, cylinder, sphere, and spherical shell. - sphere_radius : float - Outer radius of spheres. - - Returns - ------- - domain : openmc.model._Container - Container in which to pack spheres. - - """ - - def rectangular_prism(): - # Assume the simplest case where the prism volume is the intersection - # of the half-spaces of six planes - if not isinstance(region, openmc.Intersection): - return None - - if any(not isinstance(node, openmc.Halfspace) for node in region): - return None - - if len(region) != 6: - return None - - # Sort half-spaces by surface type - px1, px2, py1, py2, pz1, pz2 = sorted(region, key=attrgetter('surface.type')) - - # Make sure the region consists of the correct surfaces - if (not isinstance(px1.surface, openmc.XPlane) or - not isinstance(px2.surface, openmc.XPlane) or - not isinstance(py1.surface, openmc.YPlane) or - not isinstance(py2.surface, openmc.YPlane) or - not isinstance(pz1.surface, openmc.ZPlane) or - not isinstance(pz2.surface, openmc.ZPlane)): - return None - - # Secondary sorting by location of the plane - if px1.surface.x0 > px2.surface.x0: - px1, px2 = px2, px1 - - if py1.surface.y0 > py2.surface.y0: - py1, py2 = py2, py1 - - if pz1.surface.z0 > pz2.surface.z0: - pz1, pz2 = pz2, pz1 - - # Make sure the half-spaces are on the correct side of the surfaces - if (px1.side != '+' or px2.side != '-' or - py1.side != '+' or py2.side != '-' or - pz1.side != '+' or pz2.side != '-'): - return None - - # Calculate the parameters for the container - width = px2.surface.x0 - px1.surface.x0 - depth = py2.surface.y0 - py1.surface.y0 - height = pz2.surface.z0 - pz1.surface.z0 - center = (px1.surface.x0 + width/2, - py1.surface.y0 + depth/2, - pz1.surface.z0 + height/2) - - # The region is the volume of a rectangular prism, so create container - return _RectangularPrism(width, depth, height, sphere_radius, center) - - def cylinder(): - # Assume the simplest case where the cylinder volume is the - # intersection of the half-spaces of a cylinder and two planes - if not isinstance(region, openmc.Intersection): - return None - - if any(not isinstance(node, openmc.Halfspace) for node in region): - return None - - if len(region) != 3: - return None - - # Identify the axis that the cylinder lies along - axis = region[0].surface.type[0] - - # Make sure the region is composed of a cylinder and two planes on the - # same axis - count = Counter((node.surface.type for node in region)) - if count[axis + '-cylinder'] != 1 or count[axis + '-plane'] != 2: - return None - - # Sort the half-spaces by surface type - cyl, p1, p2 = sorted(region, key=attrgetter('surface.type')) - - # Calculate the parameters for a cylinder along the x-axis - if axis == 'x': - if p1.surface.x0 > p2.surface.x0: - p1, p2 = p2, p1 - length = p2.surface.x0 - p1.surface.x0 - center = (p1.surface.x0 + length/2, cyl.surface.y0, cyl.surface.z0) - - # Calculate the parameters for a cylinder along the y-axis - elif axis == 'y': - if p1.surface.y0 > p2.surface.y0: - p1, p2 = p2, p1 - length = p2.surface.y0 - p1.surface.y0 - center = (cyl.surface.x0, p1.surface.y0 + length/2, cyl.surface.z0) - - # Calculate the parameters for a cylinder along the z-axis - else: - if p1.surface.z0 > p2.surface.z0: - p1, p2 = p2, p1 - length = p2.surface.z0 - p1.surface.z0 - center = (cyl.surface.x0, cyl.surface.y0, p1.surface.z0 + length/2) - - # Make sure the half-spaces are on the correct side of the surfaces - if cyl.side != '-' or p1.side != '+' or p2.side != '-': - return None - - radius = cyl.surface.r - - # The region is the volume of a cylinder, so create container - return _Cylinder(length, radius, axis, sphere_radius, center) - - def sphere(): - # Assume the simplest case where the sphere volume is the negative - # half-space of a sphere - if not isinstance(region, openmc.Halfspace): - return None - - if not isinstance(region.surface, openmc.Sphere): - return None - - if region.side != '-': - return None - - radius = region.surface.r - center = (region.surface.x0, region.surface.y0, region.surface.z0) - - # The region is the volume of a sphere, so create container - return _Sphere(radius, sphere_radius, center) - - def spherical_shell(): - # Assume the simplest case where the spherical shell volume is the - # intersection of the half-spaces of two spheres - if not isinstance(region, openmc.Intersection): - return None - - if any(not isinstance(node, openmc.Halfspace) for node in region): - return None - - if len(region) != 2: - return None - - if any(not isinstance(node.surface, openmc.Sphere) for node in region): - return None - - s1, s2 = sorted(region, key=attrgetter('surface.r')) - radius = s2.surface.r - inner_radius = s1.surface.r - center = (s1.surface.x0, s1.surface.y0, s1.surface.z0) - - if center != (s2.surface.x0, s2.surface.y0, s2.surface.z0): - return None - - if s1.side != '+' or s2.side != '-': - return None - - # The region is the volume of a spherical shell, so create container - return _SphericalShell(radius, inner_radius, sphere_radius, center) - - check_type('region', region, openmc.Region) - - # Check whether the region matches any of the supported container shapes - # and create the container if it does - shapes = [rectangular_prism, cylinder, sphere, spherical_shell] - for shape in shapes: - container = shape() - if container: - return container - - msg = ('Could not translate region {} into a container: supported ' - 'container shapes are rectangular prism, cylinder, sphere, ' - 'and spherical shell.'.format(region)) - raise ValueError(msg) - - def _random_sequential_pack(domain, num_spheres): """Random sequential packing of spheres within a container. @@ -1374,7 +1272,17 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, random.seed(seed) # Create container with the correct shape based on the supplied region - domain = _create_container(region, radius) + domain = None + for cls in _Container.__subclasses__(): + try: + domain = cls.from_region(region, radius) + except ValueError: + pass + + if not domain: + raise ValueError('Could not map region {} to a container: supported ' + 'container shapes are rectangular prism, cylinder, ' + 'sphere, and spherical shell.'.format(region)) # Determine the packing fraction/number of spheres volume = 4/3*pi*radius**3 @@ -1410,7 +1318,8 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, # Recalculate the limits for the initial random sequential packing using # the desired final sphere radius to ensure spheres are fully contained # within the domain during the close random pack - domain.limits = [x + initial_radius - radius for x in domain.limits] + domain.limits = [[x - initial_radius + radius for x in domain.limits[0]], + [x + initial_radius - radius for x in domain.limits[1]]] # Generate non-overlapping spheres for an initial inner radius using # random sequential packing algorithm @@ -1423,4 +1332,4 @@ def pack_spheres(radius, region, pf=None, num_spheres=None, initial_pf=0.3, domain.sphere_radius = radius _close_random_pack(domain, spheres, contraction_rate) - return spheres + domain.center + return spheres From b7b410975b5254e7b775f4722bade0a322cfbe31 Mon Sep 17 00:00:00 2001 From: amandalund Date: Mon, 19 Nov 2018 21:56:36 -0600 Subject: [PATCH 9/9] Update TRISO unit tests --- tests/unit_tests/test_model_triso.py | 197 +++++++++++++++++---------- 1 file changed, 127 insertions(+), 70 deletions(-) diff --git a/tests/unit_tests/test_model_triso.py b/tests/unit_tests/test_model_triso.py index 3c261ed6f..1d2e0905c 100644 --- a/tests/unit_tests/test_model_triso.py +++ b/tests/unit_tests/test_model_triso.py @@ -12,52 +12,84 @@ import scipy.spatial _RADIUS = 0.1 _PACKING_FRACTION = 0.35 -_SHAPES = ['rectangular_prism', 'cylinder', 'sphere', 'spherical_shell'] -_VOLUMES = [1.**3, 1.*pi*1.**2, 4/3*pi*1.**3, 4/3*pi*(1.**3 - 0.5**3)] +_PARAMS = [ + {'shape': 'rectangular_prism', 'volume': 1**3}, + {'shape': 'x_cylinder', 'volume': 1*pi*1**2}, + {'shape': 'y_cylinder', 'volume': 1*pi*1**2}, + {'shape': 'z_cylinder', 'volume': 1*pi*1**2}, + {'shape': 'sphere', 'volume': 4/3*pi*1**3}, + {'shape': 'spherical_shell', 'volume': 4/3*pi*(1**3 - 0.5**3)} +] -@pytest.fixture(scope='module') +@pytest.fixture(scope='module', params=_PARAMS) def container(request): - return request.getfixturevalue(request.param) + return request.param @pytest.fixture(scope='module') -def centers(request): - container = request.getfixturevalue(request.param) - return openmc.model.pack_spheres(radius=_RADIUS, region=container, +def centers(request, container): + return request.getfixturevalue('centers_' + container['shape']) + + +@pytest.fixture(scope='module') +def centers_rectangular_prism(): + min_x = openmc.XPlane(x0=0) + max_x = openmc.XPlane(x0=1) + min_y = openmc.YPlane(y0=0) + max_y = openmc.YPlane(y0=1) + min_z = openmc.ZPlane(z0=0) + max_z = openmc.ZPlane(z0=1) + region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z + return openmc.model.pack_spheres(radius=_RADIUS, region=region, pf=_PACKING_FRACTION, initial_pf=0.2) @pytest.fixture(scope='module') -def rectangular_prism(): - min_x = openmc.XPlane(x0=-0.5) - max_x = openmc.XPlane(x0=0.5) - min_y = openmc.YPlane(y0=-0.5) - max_y = openmc.YPlane(y0=0.5) - min_z = openmc.ZPlane(z0=-0.5) - max_z = openmc.ZPlane(z0=0.5) - return +min_x & -max_x & +min_y & -max_y & +min_z & -max_z +def centers_x_cylinder(): + cylinder = openmc.XCylinder(R=1, y0=1, z0=2) + min_x = openmc.XPlane(x0=0) + max_x = openmc.XPlane(x0=1) + region = +min_x & -max_x & -cylinder + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) @pytest.fixture(scope='module') -def cylinder(): - cylinder = openmc.ZCylinder(R=1.) - min_z = openmc.ZPlane(z0=-0.5) - max_z = openmc.ZPlane(z0=0.5) - return +min_z & -max_z & -cylinder +def centers_y_cylinder(): + cylinder = openmc.YCylinder(R=1, x0=1, z0=2) + min_y = openmc.YPlane(y0=0) + max_y = openmc.YPlane(y0=1) + region = +min_y & -max_y & -cylinder + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) @pytest.fixture(scope='module') -def sphere(): - sphere = openmc.Sphere(R=1.) - return -sphere +def centers_z_cylinder(): + cylinder = openmc.ZCylinder(R=1, x0=1, y0=2) + min_z = openmc.ZPlane(z0=0) + max_z = openmc.ZPlane(z0=1) + region = +min_z & -max_z & -cylinder + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) @pytest.fixture(scope='module') -def spherical_shell(): - sphere = openmc.Sphere(R=1.) - inner_sphere = openmc.Sphere(R=0.5) - return -sphere & +inner_sphere +def centers_sphere(): + sphere = openmc.Sphere(R=1, x0=1, y0=2, z0=3) + region = -sphere + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) + + +@pytest.fixture(scope='module') +def centers_spherical_shell(): + sphere = openmc.Sphere(R=1, x0=1, y0=2, z0=3) + inner_sphere = openmc.Sphere(R=0.5, x0=1, y0=2, z0=3) + region = -sphere & +inner_sphere + return openmc.model.pack_spheres(radius=_RADIUS, region=region, + pf=_PACKING_FRACTION, initial_pf=0.2) @pytest.fixture(scope='module') @@ -68,7 +100,6 @@ def triso_universe(): return univ -@pytest.mark.parametrize('centers', _SHAPES, indirect=True) def test_overlap(centers): """Check that none of the spheres in the packed configuration overlap.""" # Create KD tree for quick nearest neighbor search @@ -82,57 +113,83 @@ def test_overlap(centers): assert d_min > 2*_RADIUS or d_min == pytest.approx(2*_RADIUS) -@pytest.mark.parametrize('centers,shape', zip(_SHAPES, _SHAPES), - indirect=['centers']) -def test_contained(centers, shape): +def test_contained_rectangular_prism(centers_rectangular_prism): """Make sure all spheres are entirely contained within the domain.""" - if shape == 'rectangular_prism': - x = np.amax(abs(centers)) + _RADIUS - assert x < 0.5 or x == pytest.approx(0.5) - - elif shape == 'cylinder': - r = max(np.linalg.norm(centers[:,0:2], axis=1)) + _RADIUS - z = max(abs(centers[:,2])) + _RADIUS - assert r < 1. or r == pytest.approx(1.) - assert z < 0.5 or z == pytest.approx(0.5) - - elif shape == 'sphere': - r = max(np.linalg.norm(centers, axis=1)) + _RADIUS - assert r < 1. or r == pytest.approx(1.) - - elif shape == 'spherical_shell': - d = np.linalg.norm(centers, axis=1) - r_max = max(d) + _RADIUS - r_min = min(d) - _RADIUS - assert r_max < 1. or r_max == pytest.approx(1.) - assert r_min > 0.5 or r_min == pytest.approx(0.5) + d_max = np.amax(centers_rectangular_prism) + _RADIUS + d_min = np.amin(centers_rectangular_prism) - _RADIUS + assert d_max < 1 or d_max == pytest.approx(1) + assert d_min > 0 or d_min == pytest.approx(0) -@pytest.mark.parametrize('centers,volume', zip(_SHAPES, _VOLUMES), - indirect=['centers']) -def test_packing_fraction(centers, volume): +def test_contained_x_cylinder(centers_x_cylinder): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_x_cylinder[:,[1,2]] - [1, 2], axis=1) + r_max = max(d) + _RADIUS + x_max = max(centers_x_cylinder[:,0]) + _RADIUS + x_min = min(centers_x_cylinder[:,0]) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert x_max < 1 or x_max == pytest.approx(1) + assert x_min > 0 or x_min == pytest.approx(0) + + +def test_contained_y_cylinder(centers_y_cylinder): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_y_cylinder[:,[0,2]] - [1, 2], axis=1) + r_max = max(d) + _RADIUS + y_max = max(centers_y_cylinder[:,1]) + _RADIUS + y_min = min(centers_y_cylinder[:,1]) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert y_max < 1 or y_max == pytest.approx(1) + assert y_min > 0 or y_min == pytest.approx(0) + + +def test_contained_z_cylinder(centers_z_cylinder): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_z_cylinder[:,[0,1]] - [1, 2], axis=1) + r_max = max(d) + _RADIUS + z_max = max(centers_z_cylinder[:,2]) + _RADIUS + z_min = min(centers_z_cylinder[:,2]) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert z_max < 1 or z_max == pytest.approx(1) + assert z_min > 0 or z_min == pytest.approx(0) + + +def test_contained_sphere(centers_sphere): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_sphere - [1, 2, 3], axis=1) + r_max = max(d) + _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + + +def test_contained_spherical_shell(centers_spherical_shell): + """Make sure all spheres are entirely contained within the domain.""" + d = np.linalg.norm(centers_spherical_shell - [1, 2, 3], axis=1) + r_max = max(d) + _RADIUS + r_min = min(d) - _RADIUS + assert r_max < 1 or r_max == pytest.approx(1) + assert r_min > 0.5 or r_min == pytest.approx(0.5) + + +def test_packing_fraction(container, centers): """Check that the actual PF is close to the requested PF.""" - pf = len(centers) * 4/3 * pi *_RADIUS**3 / volume + pf = len(centers) * 4/3 * pi *_RADIUS**3 / container['volume'] assert pf == pytest.approx(_PACKING_FRACTION, rel=1e-2) -@pytest.mark.parametrize('container', _SHAPES, indirect=True) -def test_num_spheres(container): +def test_num_spheres(): """Check that the function returns the correct number of spheres""" centers = openmc.model.pack_spheres( - radius=_RADIUS, region=container, num_spheres=50 + radius=_RADIUS, region=-openmc.Sphere(R=1), num_spheres=50 ) assert len(centers) == 50 -def test_triso_lattice(triso_universe, rectangular_prism): - centers = openmc.model.pack_spheres( - radius=_RADIUS, region=rectangular_prism, pf=0.2 - ) - trisos = [openmc.model.TRISO(_RADIUS, triso_universe, c) for c in centers] +def test_triso_lattice(triso_universe, centers_rectangular_prism): + trisos = [openmc.model.TRISO(_RADIUS, triso_universe, c) + for c in centers_rectangular_prism] - lower_left = np.array((-.5, -.5, -.5)) - upper_right = np.array((.5, .5, .5)) + lower_left = np.array((0, 0, 0)) + upper_right = np.array((1, 1, 1)) shape = (3, 3, 3) pitch = (upper_right - lower_left)/shape background = openmc.Material() @@ -146,25 +203,25 @@ def test_container_input(triso_universe): # Invalid container shape with pytest.raises(ValueError): centers = openmc.model.pack_spheres( - radius=_RADIUS, region=+openmc.Sphere(R=1.), num_spheres=100 + radius=_RADIUS, region=+openmc.Sphere(R=1), num_spheres=100 ) -def test_packing_fraction_input(sphere): +def test_packing_fraction_input(): # Provide neither packing fraction nor number of spheres with pytest.raises(ValueError): centers = openmc.model.pack_spheres( - radius=_RADIUS, region=sphere + radius=_RADIUS, region=-openmc.Sphere(R=1) ) # Specify a packing fraction that is too high for CRP with pytest.raises(ValueError): centers = openmc.model.pack_spheres( - radius=_RADIUS, region=sphere, pf=1. + radius=_RADIUS, region=-openmc.Sphere(R=1), pf=1 ) # Specify a packing fraction that is too high for RSP with pytest.raises(ValueError): centers = openmc.model.pack_spheres( - radius=_RADIUS, region=sphere, pf=0.5, initial_pf=0.4 + radius=_RADIUS, region=-openmc.Sphere(R=1), pf=0.5, initial_pf=0.4 )