Add initial LWR AMR test problems.

lwr/3x3_lattice/make_openmc_model.py

@@ -0,0 +1,171 @@
+#********************************************************************/
+#*                  SOFTWARE COPYRIGHT NOTIFICATION                 */
+#*                             Cardinal                             */
+#*                                                                  */
+#*                  (c) 2021 UChicago Argonne, LLC                  */
+#*                        ALL RIGHTS RESERVED                       */
+#*                                                                  */
+#*                 Prepared by UChicago Argonne, LLC                */
+#*               Under Contract No. DE-AC02-06CH11357               */
+#*                With the U. S. Department of Energy               */
+#*                                                                  */
+#*             Prepared by Battelle Energy Alliance, LLC            */
+#*               Under Contract No. DE-AC07-05ID14517               */
+#*                With the U. S. Department of Energy               */
+#*                                                                  */
+#*                 See LICENSE for full restrictions                */
+#********************************************************************/
+
+import openmc
+import numpy as np
+from argparse import ArgumentParser
+
+ap = ArgumentParser()
+ap.add_argument('-n', dest='n_axial', type=int, default=40,
+                help='Number of axial cell divisions')
+ap.add_argument('-s', '--entropy', action='store_true',
+                help='Whether to add a Shannon entropy mesh')
+args = ap.parse_args()
+
+N = args.n_axial # Number of axial cells to build in the solid to receive feedback
+height = 30.0   # Total height of pincell
+
+###############################################################################
+# Create materials for the problem
+
+uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
+uo2.set_density('g/cm3', 10.29769)
+uo2.add_element('U', 1., enrichment=2.4)
+uo2.add_element('O', 2.)
+
+helium = openmc.Material(name='Helium for gap')
+helium.set_density('g/cm3', 0.001598)
+helium.add_element('He', 2.4044e-4)
+
+zircaloy = openmc.Material(name='Zircaloy 4')
+zircaloy.set_density('g/cm3', 6.55)
+zircaloy.add_element('Sn', 0.014  , 'wo')
+zircaloy.add_element('Fe', 0.00165, 'wo')
+zircaloy.add_element('Cr', 0.001  , 'wo')
+zircaloy.add_element('Zr', 0.98335, 'wo')
+
+borated_water = openmc.Material(name='Borated water')
+borated_water.set_density('g/cm3', 0.740582)
+borated_water.add_element('B', 4.0e-5)
+borated_water.add_element('H', 5.0e-2)
+borated_water.add_element('O', 2.4e-2)
+borated_water.add_s_alpha_beta('c_H_in_H2O')
+
+# Collect the materials together and export to XML
+materials = openmc.Materials([uo2, helium, zircaloy, borated_water])
+materials.export_to_xml()
+
+###############################################################################
+# Define problem geometry
+
+# Create cylindrical surfaces
+fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
+clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
+clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
+
+# Create a region represented as the inside of a rectangular prism
+pitch = 1.25984
+box = openmc.model.RectangularPrism(pitch, pitch)
+
+# Create cells, mapping materials to regions - split up the axial height
+planes = np.linspace(0.0, height, N + 1)
+plane_surfaces = []
+for i in range(N + 1):
+  plane_surfaces.append(openmc.ZPlane(z0=planes[i]))
+
+# set the boundary condition on the topmost and bottommost planes to vacuum
+plane_surfaces[0].boundary_type = 'reflective'
+plane_surfaces[-1].boundary_type = 'vacuum'
+
+fuel_cells = []
+clad_cells = []
+gap_cells = []
+water_cells = []
+all_cells = []
+for i in range(N):
+  layer = +plane_surfaces[i] & -plane_surfaces[i + 1]
+  fuel_cells.append(openmc.Cell(fill=uo2, region=-fuel_or & layer, name='Fuel{:n}'.format(i)))
+  gap_cells.append(openmc.Cell(fill=helium, region=+fuel_or & -clad_ir & layer, name='Gap{:n}'.format(i)))
+  clad_cells.append(openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or & layer, name='Clad{:n}'.format(i)))
+  water_cells.append(openmc.Cell(fill=borated_water, region=+clad_or & layer & -box, name='Water{:n}'.format(i)))
+  all_cells.append(fuel_cells[i])
+  all_cells.append(gap_cells[i])
+  all_cells.append(clad_cells[i])
+  all_cells.append(water_cells[i])
+
+fuel_rod_universe = openmc.Universe(cells=all_cells)
+
+assembly_cells = [
+  [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],
+  [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],
+  [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe]
+]
+
+assembly = openmc.RectLattice(name='Fuel Assembly')
+assembly.pitch = (pitch, pitch)
+assembly.lower_left = (-3.0 * pitch / 2.0, -3.0 * pitch / 2.0)
+assembly.universes = assembly_cells
+
+assembly_region = openmc.model.RectangularPrism(width = 3.0 * pitch, height = 3.0 * pitch, origin = (0.0, 0.0), boundary_type = 'reflective')
+full_assembly_cell = openmc.Cell(name='Assembly Cell', fill = assembly, region=-assembly_region & -plane_surfaces[-1] & +plane_surfaces[0])
+
+# Create a geometry and export to XML
+geometry = openmc.Geometry([full_assembly_cell])
+geometry.export_to_xml()
+
+###############################################################################
+# Define problem settings
+
+# Indicate how many particles to run
+settings = openmc.Settings()
+settings.batches = 1500
+settings.inactive = 500
+settings.particles = 20000
+
+# Create an initial uniform spatial source distribution over fissionable zones
+lower_left = (-3.0 * pitch / 2.0, -3.0 * pitch / 2.0, 0.0)
+upper_right = (3.0 * pitch / 2.0, 3.0 * pitch / 2.0, height)
+uniform_dist = openmc.stats.Box(lower_left, upper_right)
+settings.source = openmc.IndependentSource(space=uniform_dist)
+
+if (args.entropy):
+  entropy_mesh = openmc.RegularMesh()
+  entropy_mesh.lower_left = lower_left
+  entropy_mesh.upper_right = upper_right
+  entropy_mesh.dimension = (10, 10, 20)
+  settings.entropy_mesh = entropy_mesh
+
+settings.temperature = {'default': 280.0 + 273.15,
+                        'method': 'interpolation',
+                        'range': (294.0, 3000.0),
+                        'tolerance': 1000.0}
+
+settings.export_to_xml()
+
+quit()
+
+# create some plots to look at the geometry for the sake of the tutorial
+plot1          = openmc.Plot()
+plot1.filename = 'plot1'
+plot1.width    = (pitch, pitch)
+plot1.basis    = 'xy'
+plot1.origin   = (0.0, 0.0, height/2.0)
+plot1.pixels   = (1000, 1000)
+plot1.color_by = 'cell'
+
+plot2          = openmc.Plot()
+plot2.filename = 'plot2'
+plot2.width    = (pitch, height)
+plot2.basis    = 'xz'
+plot2.origin   = (0.0, 0.0, height/2.0)
+plot2.pixels   = (100, int(100 * (height/2.0/pitch)))
+plot2.color_by = 'cell'
+
+plots = openmc.Plots([plot1, plot2])
+plots.export_to_xml()
+

lwr/3x3_lattice/mesh.i

@@ -0,0 +1,84 @@
+#----------------------------------------------------------------------------------------
+# Assembly geometrical information
+#----------------------------------------------------------------------------------------
+pitch        = 1.25984
+height       = 30.0
+r_fuel       = 0.39218
+t_gap        = 0.00787
+r_clad_inner = 0.40005
+t_clad       = 0.05715
+#----------------------------------------------------------------------------------------
+
+#----------------------------------------------------------------------------------------
+# Meshing parameters
+#----------------------------------------------------------------------------------------
+NUM_SECTORS              = 6
+FUEL_RADIAL_DIVISIONS    = 4
+BACKGROUND_DIVISIONS     = 2
+AXIAL_DIVISIONS          = 10
+#----------------------------------------------------------------------------------------
+
+[Mesh]
+  [Pin]
+    type = PolygonConcentricCircleMeshGenerator
+    num_sides = 4
+    num_sectors_per_side = '${NUM_SECTORS} ${NUM_SECTORS} ${NUM_SECTORS} ${NUM_SECTORS}'
+    ring_radii = '${r_fuel} ${fparse r_fuel + t_gap} ${fparse r_fuel + t_gap + t_clad}'
+    ring_intervals = '${FUEL_RADIAL_DIVISIONS} 1 1'
+    polygon_size = ${fparse pitch / 2.0}
+
+    ring_block_ids = '0 1 2 3'
+    ring_block_names = 'fuel_center fuel gap cladding'
+    background_block_ids = '4'
+    background_block_names = 'water'
+    background_intervals = ${BACKGROUND_DIVISIONS}
+
+    flat_side_up = true
+    quad_center_elements = false
+    preserve_volumes = true
+
+    create_outward_interface_boundaries = false
+  []
+  [Core_2D]
+    type = PatternedCartesianMeshGenerator
+    inputs = 'Pin'
+    pattern = '0 0 0;
+               0 0 0;
+               0 0 0'
+
+    pattern_boundary = 'none'
+
+    external_boundary_id = 0
+    external_boundary_name = 'vacuum'
+
+    assign_type = 'cell'
+    id_name = 'pin_id'
+    generate_core_metadata = false
+  []
+  [Core_3D]
+    type = AdvancedExtruderGenerator
+    input = 'Core_2D'
+    heights = '${fparse height}'
+    num_layers = '${AXIAL_DIVISIONS}'
+    direction = '0.0 0.0 1.0'
+
+    bottom_boundary = '10001'
+    top_boundary = '10000'
+  []
+  [To_Origin]
+    type = TransformGenerator
+    input = 'Core_3D'
+    transform = TRANSLATE_CENTER_ORIGIN
+  []
+  [Down]
+    type = TransformGenerator
+    input = 'To_Origin'
+    transform = TRANSLATE
+    vector_value = '0.0 0.0 ${fparse height / 2.0}'
+  []
+  [Delete_Gap]
+    type = BlockDeletionGenerator
+    input = Down
+    block = '2'
+  []
+[]

lwr/3x3_lattice/openmc.i

@@ -0,0 +1,89 @@
+[Mesh]
+  [file]
+    type = FileMeshGenerator
+    file = mesh_in.e
+  []
+[]
+
+[Adaptivity]
+  marker = error_combo
+  steps = 10
+
+  [Indicators/error]
+    type = ValueJumpIndicator
+    variable = heat_source
+  []
+  [Markers]
+    [error_frac]
+      type = ErrorFractionMarker
+      indicator = error
+      refine = 0.3
+      coarsen = 0.0
+    []
+    [rel_error]
+      type = ValueThresholdMarker
+      invert = true
+      coarsen = 1e-1
+      refine = 1e-2
+      variable = heat_source_rel_error
+      third_state = DO_NOTHING
+    []
+    [error_combo]
+      type = BooleanComboMarker
+      # Only refine iff the relative error is sufficiently low AND there is a large enough
+      # jump discontinuity in the solution.
+      refine_markers = 'rel_error error_frac'
+      # Coarsen based exclusively on relative error. Jump discontinuities in the solution
+      # from large relative errors causes the 'error_frac' marker to erroneously mark elements
+      # for refinement.
+      coarsen_markers = 'rel_error'
+      boolean_operator = and
+    []
+  []
+[]
+
+[Problem]
+  type = OpenMCCellAverageProblem
+  particles = 20000
+  inactive_batches = 500
+  batches = 10000
+
+  verbose = true
+  power = ${fparse 3000e6 / 273 / (17 * 17) * 9}
+  cell_level = 1
+  normalize_by_global_tally = false
+
+  [Tallies]
+    [heat_source]
+      type = MeshTally
+      score = 'kappa_fission'
+      name = heat_source
+      output = 'unrelaxed_tally_rel_error'
+    []
+  []
+[]
+
+[Executioner]
+  type = Steady
+[]
+
+[Postprocessors]
+  [num_active]
+    type = NumElems
+    elem_filter = active
+  []
+  [num_total]
+    type = NumElems
+    elem_filter = total
+  []
+  [max_rel_err]
+    type = TallyRelativeError
+    value_type = max
+    tally_score = kappa_fission
+  []
+[]
+
+[Outputs]
+  exodus = true
+  csv = true
+[]