From 81b859ad4ff3995ac84e27b8bf15ba2f87cc4cc4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Feb 2018 07:33:18 -0600 Subject: [PATCH] Get rid of example_run.py and example_plot.py --- scripts/example_geometry.py | 378 -------------------- scripts/example_plot.py | 43 --- scripts/example_run.py | 33 -- tests/regression_tests/example_geometry.py | 379 ++++++++++++++++++++- 4 files changed, 378 insertions(+), 455 deletions(-) delete mode 100644 scripts/example_geometry.py delete mode 100644 scripts/example_plot.py delete mode 100644 scripts/example_run.py mode change 120000 => 100644 tests/regression_tests/example_geometry.py diff --git a/scripts/example_geometry.py b/scripts/example_geometry.py deleted file mode 100644 index ca10c1f725..0000000000 --- a/scripts/example_geometry.py +++ /dev/null @@ -1,378 +0,0 @@ -"""An example file showing how to make a geometry. - -This particular example creates a 3x3 geometry, with 8 regular pins and one -Gd-157 2 wt-percent enriched. All pins are segmented. -""" - -from collections import OrderedDict -import math - -import numpy as np -import openmc - - -def density_to_mat(dens_dict): - """Generates an OpenMC material from a cell ID and self.number_density. - - Parameters - ---------- - dens_dict : dict - Dictionary mapping nuclide names to densities - - Returns - ------- - openmc.Material - The OpenMC material filled with nuclides. - - """ - mat = openmc.Material() - for key in dens_dict: - mat.add_nuclide(key, 1.0e-24*dens_dict[key]) - mat.set_density('sum') - - return mat - - -def generate_initial_number_density(): - """ Generates initial number density. - - These results were from a CASMO5 run in which the gadolinium pin was - loaded with 2 wt percent of Gd-157. - """ - - # Concentration to be used for all fuel pins - fuel_dict = OrderedDict() - fuel_dict['U235'] = 1.05692e21 - fuel_dict['U234'] = 1.00506e19 - fuel_dict['U238'] = 2.21371e22 - fuel_dict['O16'] = 4.62954e22 - fuel_dict['O17'] = 1.127684e20 - fuel_dict['I135'] = 1.0e10 - fuel_dict['Xe135'] = 1.0e10 - fuel_dict['Xe136'] = 1.0e10 - fuel_dict['Cs135'] = 1.0e10 - fuel_dict['Gd156'] = 1.0e10 - fuel_dict['Gd157'] = 1.0e10 - # fuel_dict['O18'] = 9.51352e19 # Does not exist in ENDF71, merged into 17 - - # Concentration to be used for the gadolinium fuel pin - fuel_gd_dict = OrderedDict() - fuel_gd_dict['U235'] = 1.03579e21 - fuel_gd_dict['U238'] = 2.16943e22 - fuel_gd_dict['Gd156'] = 3.95517E+10 - fuel_gd_dict['Gd157'] = 1.08156e20 - fuel_gd_dict['O16'] = 4.64035e22 - fuel_dict['I135'] = 1.0e10 - fuel_dict['Xe136'] = 1.0e10 - fuel_dict['Xe135'] = 1.0e10 - fuel_dict['Cs135'] = 1.0e10 - # There are a whole bunch of 1e-10 stuff here. - - # Concentration to be used for cladding - clad_dict = OrderedDict() - clad_dict['O16'] = 3.07427e20 - clad_dict['O17'] = 7.48868e17 - clad_dict['Cr50'] = 3.29620e18 - clad_dict['Cr52'] = 6.35639e19 - clad_dict['Cr53'] = 7.20763e18 - clad_dict['Cr54'] = 1.79413e18 - clad_dict['Fe54'] = 5.57350e18 - clad_dict['Fe56'] = 8.74921e19 - clad_dict['Fe57'] = 2.02057e18 - clad_dict['Fe58'] = 2.68901e17 - clad_dict['Cr50'] = 3.29620e18 - clad_dict['Cr52'] = 6.35639e19 - clad_dict['Cr53'] = 7.20763e18 - clad_dict['Cr54'] = 1.79413e18 - clad_dict['Ni58'] = 2.51631e19 - clad_dict['Ni60'] = 9.69278e18 - clad_dict['Ni61'] = 4.21338e17 - clad_dict['Ni62'] = 1.34341e18 - clad_dict['Ni64'] = 3.43127e17 - clad_dict['Zr90'] = 2.18320e22 - clad_dict['Zr91'] = 4.76104e21 - clad_dict['Zr92'] = 7.27734e21 - clad_dict['Zr94'] = 7.37494e21 - clad_dict['Zr96'] = 1.18814e21 - clad_dict['Sn112'] = 4.67352e18 - clad_dict['Sn114'] = 3.17992e18 - clad_dict['Sn115'] = 1.63814e18 - clad_dict['Sn116'] = 7.00546e19 - clad_dict['Sn117'] = 3.70027e19 - clad_dict['Sn118'] = 1.16694e20 - clad_dict['Sn119'] = 4.13872e19 - clad_dict['Sn120'] = 1.56973e20 - clad_dict['Sn122'] = 2.23076e19 - clad_dict['Sn124'] = 2.78966e19 - - # Gap concentration - # Funny enough, the example problem uses air. - gap_dict = OrderedDict() - gap_dict['O16'] = 7.86548e18 - gap_dict['O17'] = 2.99548e15 - gap_dict['N14'] = 3.38646e19 - gap_dict['N15'] = 1.23717e17 - - # Concentration to be used for coolant - # No boron - cool_dict = OrderedDict() - cool_dict['H1'] = 4.68063e22 - cool_dict['O16'] = 2.33427e22 - cool_dict['O17'] = 8.89086e18 - - # Store these dictionaries in the initial conditions dictionary - initial_density = OrderedDict() - initial_density['fuel_gd'] = fuel_gd_dict - initial_density['fuel'] = fuel_dict - initial_density['gap'] = gap_dict - initial_density['clad'] = clad_dict - initial_density['cool'] = cool_dict - - # Set up libraries to use - temperature = OrderedDict() - sab = OrderedDict() - - # Toggle betweeen MCNP and NNDC data - MCNP = False - - if MCNP: - temperature['fuel_gd'] = 900.0 - temperature['fuel'] = 900.0 - # We approximate temperature of everything as 600K, even though it was - # actually 580K. - temperature['gap'] = 600.0 - temperature['clad'] = 600.0 - temperature['cool'] = 600.0 - else: - temperature['fuel_gd'] = 293.6 - temperature['fuel'] = 293.6 - temperature['gap'] = 293.6 - temperature['clad'] = 293.6 - temperature['cool'] = 293.6 - - sab['cool'] = 'c_H_in_H2O' - - # Set up burnable materials - burn = OrderedDict() - burn['fuel_gd'] = True - burn['fuel'] = True - burn['gap'] = False - burn['clad'] = False - burn['cool'] = False - - return temperature, sab, initial_density, burn - -def segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad): - """ Calculates a segmented pin. - - Separates a pin with n_rings and n_wedges. All cells have equal volume. - Pin is centered at origin. - """ - - # Calculate all the volumes of interest - v_fuel = math.pi * r_fuel**2 - v_gap = math.pi * r_gap**2 - v_fuel - v_clad = math.pi * r_clad**2 - v_fuel - v_gap - v_ring = v_fuel / n_rings - v_segment = v_ring / n_wedges - - # Compute ring radiuses - r_rings = np.zeros(n_rings) - - for i in range(n_rings): - r_rings[i] = math.sqrt(1.0/(math.pi) * v_ring * (i+1)) - - # Compute thetas - theta = np.linspace(0, 2*math.pi, n_wedges + 1) - - # Compute surfaces - fuel_rings = [openmc.ZCylinder(x0=0, y0=0, R=r_rings[i]) - for i in range(n_rings)] - - fuel_wedges = [openmc.Plane(A=math.cos(theta[i]), B=math.sin(theta[i])) - for i in range(n_wedges)] - - gap_ring = openmc.ZCylinder(x0=0, y0=0, R=r_gap) - clad_ring = openmc.ZCylinder(x0=0, y0=0, R=r_clad) - - # Create cells - fuel_cells = [] - if n_wedges == 1: - for i in range(n_rings): - cell = openmc.Cell(name='fuel') - if i == 0: - cell.region = -fuel_rings[0] - else: - cell.region = +fuel_rings[i-1] & -fuel_rings[i] - fuel_cells.append(cell) - else: - for i in range(n_rings): - for j in range(n_wedges): - cell = openmc.Cell(name='fuel') - if i == 0: - if j != n_wedges-1: - cell.region = (-fuel_rings[0] - & +fuel_wedges[j] - & -fuel_wedges[j+1]) - else: - cell.region = (-fuel_rings[0] - & +fuel_wedges[j] - & -fuel_wedges[0]) - else: - if j != n_wedges-1: - cell.region = (+fuel_rings[i-1] - & -fuel_rings[i] - & +fuel_wedges[j] - & -fuel_wedges[j+1]) - else: - cell.region = (+fuel_rings[i-1] - & -fuel_rings[i] - & +fuel_wedges[j] - & -fuel_wedges[0]) - fuel_cells.append(cell) - - # Gap ring - gap_cell = openmc.Cell(name='gap') - gap_cell.region = +fuel_rings[-1] & -gap_ring - fuel_cells.append(gap_cell) - - # Clad ring - clad_cell = openmc.Cell(name='clad') - clad_cell.region = +gap_ring & -clad_ring - fuel_cells.append(clad_cell) - - # Moderator - mod_cell = openmc.Cell(name='cool') - mod_cell.region = +clad_ring - fuel_cells.append(mod_cell) - - # Form universe - fuel_u = openmc.Universe() - fuel_u.add_cells(fuel_cells) - - return fuel_u, v_segment, v_gap, v_clad - -def generate_geometry(n_rings, n_wedges): - """ Generates example geometry. - - This function creates the initial geometry, a 9 pin reflective problem. - One pin, containing gadolinium, is discretized into sectors. - - In addition to what one would do with the general OpenMC geometry code, it - is necessary to create a dictionary, volume, that maps a cell ID to a - volume. Further, by naming cells the same as the above materials, the code - can automatically handle the mapping. - - Parameters - ---------- - n_rings : int - Number of rings to generate for the geometry - n_wedges : int - Number of wedges to generate for the geometry - """ - - pitch = 1.26197 - r_fuel = 0.412275 - r_gap = 0.418987 - r_clad = 0.476121 - - n_pin = 3 - - # This table describes the 'fuel' to actual type mapping - # It's not necessary to do it this way. Just adjust the initial conditions - # below. - mapping = ['fuel', 'fuel', 'fuel', - 'fuel', 'fuel_gd', 'fuel', - 'fuel', 'fuel', 'fuel'] - - # Form pin cell - fuel_u, v_segment, v_gap, v_clad = segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad) - - # Form lattice - all_water_c = openmc.Cell(name='cool') - all_water_u = openmc.Universe(cells=(all_water_c, )) - - lattice = openmc.RectLattice() - lattice.pitch = [pitch]*2 - lattice.lower_left = [-pitch*n_pin/2, -pitch*n_pin/2] - lattice_array = [[fuel_u for i in range(n_pin)] for j in range(n_pin)] - lattice.universes = lattice_array - lattice.outer = all_water_u - - # Bound universe - x_low = openmc.XPlane(x0=-pitch*n_pin/2, boundary_type='reflective') - x_high = openmc.XPlane(x0=pitch*n_pin/2, boundary_type='reflective') - y_low = openmc.YPlane(y0=-pitch*n_pin/2, boundary_type='reflective') - y_high = openmc.YPlane(y0=pitch*n_pin/2, boundary_type='reflective') - z_low = openmc.ZPlane(z0=-10, boundary_type='reflective') - z_high = openmc.ZPlane(z0=10, boundary_type='reflective') - - # Compute bounding box - lower_left = [-pitch*n_pin/2, -pitch*n_pin/2, -10] - upper_right = [pitch*n_pin/2, pitch*n_pin/2, 10] - - root_c = openmc.Cell(fill=lattice) - root_c.region = (+x_low & -x_high - & +y_low & -y_high - & +z_low & -z_high) - root_u = openmc.Universe(universe_id=0, cells=(root_c, )) - geometry = openmc.Geometry(root_u) - - v_cool = pitch**2 - (v_gap + v_clad + n_rings * n_wedges * v_segment) - - # Store volumes for later usage - volume = {'fuel': v_segment, 'gap':v_gap, 'clad':v_clad, 'cool':v_cool} - - return geometry, volume, mapping, lower_left, upper_right - -def generate_problem(n_rings=5, n_wedges=8): - """ Merges geometry and materials. - - This function initializes the materials for each cell using the dictionaries - provided by generate_initial_number_density. It is assumed a cell named - 'fuel' will have further region differentiation (see mapping). - - Parameters - ---------- - n_rings : int, optional - Number of rings to generate for the geometry - n_wedges : int, optional - Number of wedges to generate for the geometry - """ - - # Get materials dictionary, geometry, and volumes - temperature, sab, initial_density, burn = generate_initial_number_density() - geometry, volume, mapping, lower_left, upper_right = generate_geometry(n_rings, n_wedges) - - # Apply distribmats, fill geometry - cells = geometry.root_universe.get_all_cells() - for cell_id in cells: - cell = cells[cell_id] - if cell.name == 'fuel': - - omc_mats = [] - - for cell_type in mapping: - omc_mat = density_to_mat(initial_density[cell_type]) - - if cell_type in sab: - omc_mat.add_s_alpha_beta(sab[cell_type]) - omc_mat.temperature = temperature[cell_type] - omc_mat.depletable = burn[cell_type] - omc_mat.volume = volume['fuel'] - - omc_mats.append(omc_mat) - - cell.fill = omc_mats - elif cell.name != '': - omc_mat = density_to_mat(initial_density[cell.name]) - - if cell.name in sab: - omc_mat.add_s_alpha_beta(sab[cell.name]) - omc_mat.temperature = temperature[cell.name] - omc_mat.depletable = burn[cell.name] - omc_mat.volume = volume[cell.name] - - cell.fill = omc_mat - - return geometry, lower_left, upper_right diff --git a/scripts/example_plot.py b/scripts/example_plot.py deleted file mode 100644 index ab5ac204d8..0000000000 --- a/scripts/example_plot.py +++ /dev/null @@ -1,43 +0,0 @@ -"""An example file showing how to plot data from a simulation.""" - -import matplotlib.pyplot as plt -from openmc.deplete import (read_results, evaluate_single_nuclide, - evaluate_reaction_rate, evaluate_eigenvalue) - -# Set variables for where the data is, and what we want to read out. -result_folder = "test" - -# Load data -results = read_results(result_folder + "/deplete_results.h5") - -cell = "5" -nuc = "Gd157" -rxn = "(n,gamma)" - -# Total number of nuclides -plt.figure() -# Pointwise data -x, y = evaluate_single_nuclide(results, cell, nuc) -plt.semilogy(x, y) - -plt.xlabel("Time, s") -plt.ylabel("Total Number") -plt.savefig("number.pdf") - -# Reaction rate -plt.figure() -x, y = evaluate_reaction_rate(results, cell, nuc, rxn) -plt.plot(x, y) -plt.xlabel("Time, s") -plt.ylabel("Reaction Rate, 1/s") - -plt.savefig("rate.pdf") - -# Eigenvalue -plt.figure() -x, y = evaluate_eigenvalue(results) -plt.plot(x, y) -plt.xlabel("Time, s") -plt.ylabel("Eigenvalue") - -plt.savefig("eigvl.pdf") diff --git a/scripts/example_run.py b/scripts/example_run.py deleted file mode 100644 index 36b6cce1a4..0000000000 --- a/scripts/example_run.py +++ /dev/null @@ -1,33 +0,0 @@ -"""An example file showing how to run a simulation.""" - -import numpy as np -import openmc -from openmc.data import JOULE_PER_EV -import openmc.deplete - -import example_geometry - -# Load geometry from example -geometry, lower_left, upper_right = example_geometry.generate_problem() - -# Create dt vector for 5.5 months with 15 day timesteps -dt1 = 15*24*60*60 # 15 days -dt2 = 5.5*30*24*60*60 # 5.5 months -N = np.floor(dt2/dt1) -dt = np.repeat([dt1], N) - -# Power for simulation -power = 2.337e15*4*JOULE_PER_EV*1e6 # MeV/second cm from CASMO - -# OpenMC settings -settings = openmc.Settings() -settings.particles = 1000 -settings.batches = 100 -settings.inactive = 40 -settings.source = openmc.Source(space=openmc.stats.Box(lower_left, upper_right)) - -op = openmc.deplete.Operator(geometry, settings) -op.output_dir = 'test' - -# Perform simulation using the MCNPX/MCNP6 algorithm -openmc.deplete.integrator.cecm(op, dt, power) diff --git a/tests/regression_tests/example_geometry.py b/tests/regression_tests/example_geometry.py deleted file mode 120000 index 1071aabc05..0000000000 --- a/tests/regression_tests/example_geometry.py +++ /dev/null @@ -1 +0,0 @@ -../../scripts/example_geometry.py \ No newline at end of file diff --git a/tests/regression_tests/example_geometry.py b/tests/regression_tests/example_geometry.py new file mode 100644 index 0000000000..ca10c1f725 --- /dev/null +++ b/tests/regression_tests/example_geometry.py @@ -0,0 +1,378 @@ +"""An example file showing how to make a geometry. + +This particular example creates a 3x3 geometry, with 8 regular pins and one +Gd-157 2 wt-percent enriched. All pins are segmented. +""" + +from collections import OrderedDict +import math + +import numpy as np +import openmc + + +def density_to_mat(dens_dict): + """Generates an OpenMC material from a cell ID and self.number_density. + + Parameters + ---------- + dens_dict : dict + Dictionary mapping nuclide names to densities + + Returns + ------- + openmc.Material + The OpenMC material filled with nuclides. + + """ + mat = openmc.Material() + for key in dens_dict: + mat.add_nuclide(key, 1.0e-24*dens_dict[key]) + mat.set_density('sum') + + return mat + + +def generate_initial_number_density(): + """ Generates initial number density. + + These results were from a CASMO5 run in which the gadolinium pin was + loaded with 2 wt percent of Gd-157. + """ + + # Concentration to be used for all fuel pins + fuel_dict = OrderedDict() + fuel_dict['U235'] = 1.05692e21 + fuel_dict['U234'] = 1.00506e19 + fuel_dict['U238'] = 2.21371e22 + fuel_dict['O16'] = 4.62954e22 + fuel_dict['O17'] = 1.127684e20 + fuel_dict['I135'] = 1.0e10 + fuel_dict['Xe135'] = 1.0e10 + fuel_dict['Xe136'] = 1.0e10 + fuel_dict['Cs135'] = 1.0e10 + fuel_dict['Gd156'] = 1.0e10 + fuel_dict['Gd157'] = 1.0e10 + # fuel_dict['O18'] = 9.51352e19 # Does not exist in ENDF71, merged into 17 + + # Concentration to be used for the gadolinium fuel pin + fuel_gd_dict = OrderedDict() + fuel_gd_dict['U235'] = 1.03579e21 + fuel_gd_dict['U238'] = 2.16943e22 + fuel_gd_dict['Gd156'] = 3.95517E+10 + fuel_gd_dict['Gd157'] = 1.08156e20 + fuel_gd_dict['O16'] = 4.64035e22 + fuel_dict['I135'] = 1.0e10 + fuel_dict['Xe136'] = 1.0e10 + fuel_dict['Xe135'] = 1.0e10 + fuel_dict['Cs135'] = 1.0e10 + # There are a whole bunch of 1e-10 stuff here. + + # Concentration to be used for cladding + clad_dict = OrderedDict() + clad_dict['O16'] = 3.07427e20 + clad_dict['O17'] = 7.48868e17 + clad_dict['Cr50'] = 3.29620e18 + clad_dict['Cr52'] = 6.35639e19 + clad_dict['Cr53'] = 7.20763e18 + clad_dict['Cr54'] = 1.79413e18 + clad_dict['Fe54'] = 5.57350e18 + clad_dict['Fe56'] = 8.74921e19 + clad_dict['Fe57'] = 2.02057e18 + clad_dict['Fe58'] = 2.68901e17 + clad_dict['Cr50'] = 3.29620e18 + clad_dict['Cr52'] = 6.35639e19 + clad_dict['Cr53'] = 7.20763e18 + clad_dict['Cr54'] = 1.79413e18 + clad_dict['Ni58'] = 2.51631e19 + clad_dict['Ni60'] = 9.69278e18 + clad_dict['Ni61'] = 4.21338e17 + clad_dict['Ni62'] = 1.34341e18 + clad_dict['Ni64'] = 3.43127e17 + clad_dict['Zr90'] = 2.18320e22 + clad_dict['Zr91'] = 4.76104e21 + clad_dict['Zr92'] = 7.27734e21 + clad_dict['Zr94'] = 7.37494e21 + clad_dict['Zr96'] = 1.18814e21 + clad_dict['Sn112'] = 4.67352e18 + clad_dict['Sn114'] = 3.17992e18 + clad_dict['Sn115'] = 1.63814e18 + clad_dict['Sn116'] = 7.00546e19 + clad_dict['Sn117'] = 3.70027e19 + clad_dict['Sn118'] = 1.16694e20 + clad_dict['Sn119'] = 4.13872e19 + clad_dict['Sn120'] = 1.56973e20 + clad_dict['Sn122'] = 2.23076e19 + clad_dict['Sn124'] = 2.78966e19 + + # Gap concentration + # Funny enough, the example problem uses air. + gap_dict = OrderedDict() + gap_dict['O16'] = 7.86548e18 + gap_dict['O17'] = 2.99548e15 + gap_dict['N14'] = 3.38646e19 + gap_dict['N15'] = 1.23717e17 + + # Concentration to be used for coolant + # No boron + cool_dict = OrderedDict() + cool_dict['H1'] = 4.68063e22 + cool_dict['O16'] = 2.33427e22 + cool_dict['O17'] = 8.89086e18 + + # Store these dictionaries in the initial conditions dictionary + initial_density = OrderedDict() + initial_density['fuel_gd'] = fuel_gd_dict + initial_density['fuel'] = fuel_dict + initial_density['gap'] = gap_dict + initial_density['clad'] = clad_dict + initial_density['cool'] = cool_dict + + # Set up libraries to use + temperature = OrderedDict() + sab = OrderedDict() + + # Toggle betweeen MCNP and NNDC data + MCNP = False + + if MCNP: + temperature['fuel_gd'] = 900.0 + temperature['fuel'] = 900.0 + # We approximate temperature of everything as 600K, even though it was + # actually 580K. + temperature['gap'] = 600.0 + temperature['clad'] = 600.0 + temperature['cool'] = 600.0 + else: + temperature['fuel_gd'] = 293.6 + temperature['fuel'] = 293.6 + temperature['gap'] = 293.6 + temperature['clad'] = 293.6 + temperature['cool'] = 293.6 + + sab['cool'] = 'c_H_in_H2O' + + # Set up burnable materials + burn = OrderedDict() + burn['fuel_gd'] = True + burn['fuel'] = True + burn['gap'] = False + burn['clad'] = False + burn['cool'] = False + + return temperature, sab, initial_density, burn + +def segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad): + """ Calculates a segmented pin. + + Separates a pin with n_rings and n_wedges. All cells have equal volume. + Pin is centered at origin. + """ + + # Calculate all the volumes of interest + v_fuel = math.pi * r_fuel**2 + v_gap = math.pi * r_gap**2 - v_fuel + v_clad = math.pi * r_clad**2 - v_fuel - v_gap + v_ring = v_fuel / n_rings + v_segment = v_ring / n_wedges + + # Compute ring radiuses + r_rings = np.zeros(n_rings) + + for i in range(n_rings): + r_rings[i] = math.sqrt(1.0/(math.pi) * v_ring * (i+1)) + + # Compute thetas + theta = np.linspace(0, 2*math.pi, n_wedges + 1) + + # Compute surfaces + fuel_rings = [openmc.ZCylinder(x0=0, y0=0, R=r_rings[i]) + for i in range(n_rings)] + + fuel_wedges = [openmc.Plane(A=math.cos(theta[i]), B=math.sin(theta[i])) + for i in range(n_wedges)] + + gap_ring = openmc.ZCylinder(x0=0, y0=0, R=r_gap) + clad_ring = openmc.ZCylinder(x0=0, y0=0, R=r_clad) + + # Create cells + fuel_cells = [] + if n_wedges == 1: + for i in range(n_rings): + cell = openmc.Cell(name='fuel') + if i == 0: + cell.region = -fuel_rings[0] + else: + cell.region = +fuel_rings[i-1] & -fuel_rings[i] + fuel_cells.append(cell) + else: + for i in range(n_rings): + for j in range(n_wedges): + cell = openmc.Cell(name='fuel') + if i == 0: + if j != n_wedges-1: + cell.region = (-fuel_rings[0] + & +fuel_wedges[j] + & -fuel_wedges[j+1]) + else: + cell.region = (-fuel_rings[0] + & +fuel_wedges[j] + & -fuel_wedges[0]) + else: + if j != n_wedges-1: + cell.region = (+fuel_rings[i-1] + & -fuel_rings[i] + & +fuel_wedges[j] + & -fuel_wedges[j+1]) + else: + cell.region = (+fuel_rings[i-1] + & -fuel_rings[i] + & +fuel_wedges[j] + & -fuel_wedges[0]) + fuel_cells.append(cell) + + # Gap ring + gap_cell = openmc.Cell(name='gap') + gap_cell.region = +fuel_rings[-1] & -gap_ring + fuel_cells.append(gap_cell) + + # Clad ring + clad_cell = openmc.Cell(name='clad') + clad_cell.region = +gap_ring & -clad_ring + fuel_cells.append(clad_cell) + + # Moderator + mod_cell = openmc.Cell(name='cool') + mod_cell.region = +clad_ring + fuel_cells.append(mod_cell) + + # Form universe + fuel_u = openmc.Universe() + fuel_u.add_cells(fuel_cells) + + return fuel_u, v_segment, v_gap, v_clad + +def generate_geometry(n_rings, n_wedges): + """ Generates example geometry. + + This function creates the initial geometry, a 9 pin reflective problem. + One pin, containing gadolinium, is discretized into sectors. + + In addition to what one would do with the general OpenMC geometry code, it + is necessary to create a dictionary, volume, that maps a cell ID to a + volume. Further, by naming cells the same as the above materials, the code + can automatically handle the mapping. + + Parameters + ---------- + n_rings : int + Number of rings to generate for the geometry + n_wedges : int + Number of wedges to generate for the geometry + """ + + pitch = 1.26197 + r_fuel = 0.412275 + r_gap = 0.418987 + r_clad = 0.476121 + + n_pin = 3 + + # This table describes the 'fuel' to actual type mapping + # It's not necessary to do it this way. Just adjust the initial conditions + # below. + mapping = ['fuel', 'fuel', 'fuel', + 'fuel', 'fuel_gd', 'fuel', + 'fuel', 'fuel', 'fuel'] + + # Form pin cell + fuel_u, v_segment, v_gap, v_clad = segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad) + + # Form lattice + all_water_c = openmc.Cell(name='cool') + all_water_u = openmc.Universe(cells=(all_water_c, )) + + lattice = openmc.RectLattice() + lattice.pitch = [pitch]*2 + lattice.lower_left = [-pitch*n_pin/2, -pitch*n_pin/2] + lattice_array = [[fuel_u for i in range(n_pin)] for j in range(n_pin)] + lattice.universes = lattice_array + lattice.outer = all_water_u + + # Bound universe + x_low = openmc.XPlane(x0=-pitch*n_pin/2, boundary_type='reflective') + x_high = openmc.XPlane(x0=pitch*n_pin/2, boundary_type='reflective') + y_low = openmc.YPlane(y0=-pitch*n_pin/2, boundary_type='reflective') + y_high = openmc.YPlane(y0=pitch*n_pin/2, boundary_type='reflective') + z_low = openmc.ZPlane(z0=-10, boundary_type='reflective') + z_high = openmc.ZPlane(z0=10, boundary_type='reflective') + + # Compute bounding box + lower_left = [-pitch*n_pin/2, -pitch*n_pin/2, -10] + upper_right = [pitch*n_pin/2, pitch*n_pin/2, 10] + + root_c = openmc.Cell(fill=lattice) + root_c.region = (+x_low & -x_high + & +y_low & -y_high + & +z_low & -z_high) + root_u = openmc.Universe(universe_id=0, cells=(root_c, )) + geometry = openmc.Geometry(root_u) + + v_cool = pitch**2 - (v_gap + v_clad + n_rings * n_wedges * v_segment) + + # Store volumes for later usage + volume = {'fuel': v_segment, 'gap':v_gap, 'clad':v_clad, 'cool':v_cool} + + return geometry, volume, mapping, lower_left, upper_right + +def generate_problem(n_rings=5, n_wedges=8): + """ Merges geometry and materials. + + This function initializes the materials for each cell using the dictionaries + provided by generate_initial_number_density. It is assumed a cell named + 'fuel' will have further region differentiation (see mapping). + + Parameters + ---------- + n_rings : int, optional + Number of rings to generate for the geometry + n_wedges : int, optional + Number of wedges to generate for the geometry + """ + + # Get materials dictionary, geometry, and volumes + temperature, sab, initial_density, burn = generate_initial_number_density() + geometry, volume, mapping, lower_left, upper_right = generate_geometry(n_rings, n_wedges) + + # Apply distribmats, fill geometry + cells = geometry.root_universe.get_all_cells() + for cell_id in cells: + cell = cells[cell_id] + if cell.name == 'fuel': + + omc_mats = [] + + for cell_type in mapping: + omc_mat = density_to_mat(initial_density[cell_type]) + + if cell_type in sab: + omc_mat.add_s_alpha_beta(sab[cell_type]) + omc_mat.temperature = temperature[cell_type] + omc_mat.depletable = burn[cell_type] + omc_mat.volume = volume['fuel'] + + omc_mats.append(omc_mat) + + cell.fill = omc_mats + elif cell.name != '': + omc_mat = density_to_mat(initial_density[cell.name]) + + if cell.name in sab: + omc_mat.add_s_alpha_beta(sab[cell.name]) + omc_mat.temperature = temperature[cell.name] + omc_mat.depletable = burn[cell.name] + omc_mat.volume = volume[cell.name] + + cell.fill = omc_mat + + return geometry, lower_left, upper_right