Reopening of Issue #971: 3D Tally-Score Capability

SterlingButters commented 5 years ago

Greetings @paulromano,

I wanted to revisit my issue from long ago. Points of Interest:

"running openmc-plot-mesh-tally on the statepoint file even lets me choose the axial level for a specified score which is exactly what I want. (With returned values/plot obviously)"

"The purpose of the now-defunct openmc-statepoint-3d script was to create SILO or VTK files based on a statepoint containing a 3D mesh tally, so that particularly capability is not available at the moment. However, you can definitely still do tallies containing 3D mesh filters without issue."

So my new question is if 3D statepoint evaluation is deprecated then how is openmc-plot-mesh-tally able to parse axial slices?

Adapting one of the example models into a more "complete" function and running the simulation:

def pwr_assembly():
    """Create a PWR assembly model.

    This model is a reflected 17x17 fuel assembly from the the `BEAVRS
    <http://crpg.mit.edu/research/beavrs>`_ benchmark. The fuel is 2.4 w/o
    enriched UO2 corresponding to a beginning-of-cycle condition. Note that the
    number of particles/batches is initially set very low for testing purposes.

    Returns
    -------
    model : openmc.model.Model
        A PWR assembly model

    """

    model = openmc.model.Model()

    ############################################################################################################
    # MATERIALS

    # Define materials.
    fuel = openmc.Material(name='Fuel')
    fuel.set_density('g/cm3', 10.29769)
    fuel.add_nuclide("U234", 4.4843e-6)
    fuel.add_nuclide("U235", 5.5815e-4)
    fuel.add_nuclide("U238", 2.2408e-2)
    fuel.add_nuclide("O16", 4.5829e-2)

    clad = openmc.Material(name='Cladding')
    clad.set_density('g/cm3', 6.55)
    clad.add_nuclide("Zr90", 2.1827e-2)
    clad.add_nuclide("Zr91", 4.7600e-3)
    clad.add_nuclide("Zr92", 7.2758e-3)
    clad.add_nuclide("Zr94", 7.3734e-3)
    clad.add_nuclide("Zr96", 1.1879e-3)

    hot_water = openmc.Material(name='Hot borated water')
    hot_water.set_density('g/cm3', 0.740582)
    hot_water.add_nuclide("H1", 4.9457e-2)
    hot_water.add_nuclide("O16", 2.4672e-2)
    hot_water.add_nuclide("B10", 8.0042e-6)
    hot_water.add_nuclide("B11", 3.2218e-5)
    hot_water.add_s_alpha_beta('c_H_in_H2O')

    # Define the materials file.
    model.materials = (fuel, clad, hot_water)

    ############################################################################################################
    # GEOMETRY

    # Instantiate ZCylinder surfaces
    fuel_or = openmc.ZCylinder(x0=0, y0=0, R=0.39218, name='Fuel OR')
    clad_or = openmc.ZCylinder(x0=0, y0=0, R=0.45720, name='Clad OR')

    # Create boundary planes to surround the geometry
    pitch = 21.42
    height = 200
    min_x = openmc.XPlane(x0=-pitch/2, boundary_type='reflective')
    max_x = openmc.XPlane(x0=+pitch/2, boundary_type='reflective')
    min_y = openmc.YPlane(y0=-pitch/2, boundary_type='reflective')
    max_y = openmc.YPlane(y0=+pitch/2, boundary_type='reflective')
    min_z = openmc.ZPlane(z0=-height/2, boundary_type='reflective')
    max_z = openmc.ZPlane(z0=+height/2, boundary_type='reflective')

    # Create a control rod guide tube universe
    guide_tube_universe = openmc.Universe(name='Guide Tube')
    gt_inner_cell = openmc.Cell(name='guide tube inner water', fill=hot_water,
                                region=-fuel_or)
    gt_clad_cell = openmc.Cell(name='guide tube clad', fill=clad,
                               region=+fuel_or & -clad_or)
    gt_outer_cell = openmc.Cell(name='guide tube outer water', fill=hot_water,
                                region=+clad_or)
    guide_tube_universe.add_cells([gt_inner_cell, gt_clad_cell, gt_outer_cell])

    # Create fuel assembly Lattice
    assembly = openmc.RectLattice(name='Fuel Assembly')
    assembly.pitch = (pitch/17, pitch/17)
    assembly.lower_left = (-pitch/2, -pitch/2)

    # Create array indices for guide tube locations in lattice
    template_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8,
                           11, 14, 2, 5, 8, 11, 14, 3, 13, 5, 8, 11])
    template_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8,
                           8, 11, 11, 11, 11, 11, 13, 13, 14, 14, 14])

    # Create 17x17 array of universes
    fuel_pin_universe = openmc.Universe(name='Fuel Pin')
    fuel_cell = openmc.Cell(name='fuel', fill=fuel, region=-fuel_or)
    clad_cell = openmc.Cell(name='clad', fill=clad, region=+fuel_or & -clad_or)
    hot_water_cell = openmc.Cell(name='hot water', fill=hot_water, region=+clad_or)
    fuel_pin_universe.add_cells([fuel_cell, clad_cell, hot_water_cell])

    assembly.universes = np.tile(fuel_pin_universe, (17, 17))
    assembly.universes[template_x, template_y] = guide_tube_universe

    # Create root Cell
    root_cell = openmc.Cell(name='root cell')
    root_cell.fill = assembly
    root_cell.region = +min_x & -max_x & \
                       +min_y & -max_y & \
                       +min_z & -max_z

    # Create root Universe
    model.geometry.root_universe = openmc.Universe(name='root universe')
    model.geometry.root_universe.add_cell(root_cell)

    ############################################################################################################
    # MESH

    mesh = openmc.Mesh(mesh_id=1)
    mesh.type = 'regular'
    dim_x, dim_y, dim_z = [100, 100, 10]
    mesh.dimension = [dim_x, dim_y, dim_z]
    mesh.lower_left = [*assembly.lower_left, -height / 2]
    mesh.width = [*assembly.pitch, height / dim_z]

    # Create a mesh filter
    mesh_filter = openmc.MeshFilter(mesh)

    ############################################################################################################
    # CROSS-SECTION LIBRARY

    energy_groups = openmc.mgxs.EnergyGroups()
    energy_groups.group_edges = np.logspace(-3, 7.3, 21)

    # Instantiate a 1-group EnergyGroups object
    one_group = openmc.mgxs.EnergyGroups()
    one_group.group_edges = np.array([energy_groups.group_edges[0], energy_groups.group_edges[-1]])

    # Initialize an 20-energy-group and 6-delayed-group MGXS Library
    mgxs_lib = openmc.mgxs.Library(model.geometry)
    mgxs_lib.energy_groups = one_group
    mgxs_lib.num_delayed_groups = 6

    # Specify multi-group cross section types to compute
    mgxs_lib.mgxs_types = ['total', 'transport', 'nu-scatter matrix', 'kappa-fission', 'inverse-velocity', 'chi-prompt',
                           'prompt-nu-fission', 'chi-delayed', 'delayed-nu-fission', 'beta']
    # Specify a "mesh" domain type for the cross section tally filters
    mgxs_lib.domain_type = 'mesh'
    # Specify the mesh domain over which to compute multi-group cross sections
    mgxs_lib.domains = [mesh]

    # Construct all tallies needed for the multi-group cross section library
    mgxs_lib.build_library()

    ############################################################################################################
    # TALLIES

    # Create a "tallies.xml" file for the MGXS Library
    model.tallies = openmc.Tallies()
    mgxs_lib.add_to_tallies_file(model.tallies, merge=True)

    # Instantiate a flux tally; Other valid options: 'current', 'fission', etc
    flux_tally = openmc.Tally(name='Flux')
    flux_tally.filters = [mesh_filter]
    flux_tally.scores = ['flux']

    # Add tallies to the tallies file
    model.tallies.append(flux_tally)

    ############################################################################################################
    # SETTINGS

    model.settings.batches = 10
    model.settings.inactive = 5
    model.settings.particles = 100
    model.settings.source = openmc.Source(space=openmc.stats.Box(
        [-pitch/2, -pitch/2, -height/2], [pitch/2, pitch/2, height/2], only_fissionable=True))

    ############################################################################################################
    # PLOT

    plot = openmc.Plot()
    plot.origin = (0.0, 0.0, 0)
    plot.width = (21.42, 21.42)
    plot.pixels = (300, 300)
    plot.color_by = 'material'
    model.plots.append(plot)

    openmc.plot_geometry(output=False)

    return model

pwr_assembly().export_to_xml()
openmc.run(output=True)

yields the statepoint file statepoint.10.h5.

Then using openmc-plot-mesh-tally, I am able to get what appears to be acurate axial slices of the assembly (see attachments).

SterlingButters commented 5 years ago

The real puzzler that initiated this question is as follows:

I have the following, more "complete" model of the full core example:

def pwr_core():
    """Create a PWR full-core model.

    This model is the OECD/NEA Monte Carlo Performance benchmark which is a
    grossly simplified pressurized water reactor (PWR) with 241 fuel
    assemblies. Note that the number of particles/batches is initially set very
    low for testing purposes.

    Returns
    -------
    model : openmc.model.Model
        Full-core PWR model

    """
    model = openmc.model.Model()

    ############################################################################################################
    # MATERIALS

    # Define materials.

    # 1 Fuel
    fuel = openmc.Material(1, name='UOX fuel')
    fuel.depletable = True
    fuel.set_density('g/cm3', 10.062)
    fuel.add_nuclide("U234", 4.9476e-6)
    fuel.add_nuclide("U235", 4.8218e-4)
    fuel.add_nuclide("U238", 2.1504e-2)
    fuel.add_nuclide("Xe135", 1.0801e-8)
    fuel.add_nuclide("O16", 4.5737e-2)

    clad = openmc.Material(2, name='Zircaloy')
    clad.set_density('g/cm3', 5.77)
    clad.add_nuclide("Zr90", 0.5145)
    clad.add_nuclide("Zr91", 0.1122)
    clad.add_nuclide("Zr92", 0.1715)
    clad.add_nuclide("Zr94", 0.1738)
    clad.add_nuclide("Zr96", 0.0280)

    cold_water = openmc.Material(3, name='Cold borated water')
    cold_water.set_density('atom/b-cm', 0.07416)
    cold_water.add_nuclide("H1", 2.0)
    cold_water.add_nuclide("O16", 1.0)
    cold_water.add_nuclide("B10", 6.490e-4)
    cold_water.add_nuclide("B11", 2.689e-3)
    cold_water.add_s_alpha_beta('c_H_in_H2O')

    hot_water = openmc.Material(4, name='Hot borated water')
    hot_water.set_density('atom/b-cm', 0.06614)
    hot_water.add_nuclide("H1", 2.0)
    hot_water.add_nuclide("O16", 1.0)
    hot_water.add_nuclide("B10", 6.490e-4)
    hot_water.add_nuclide("B11", 2.689e-3)
    hot_water.add_s_alpha_beta('c_H_in_H2O')

    rpv_steel = openmc.Material(5, name='Reactor pressure vessel steel')
    rpv_steel.set_density('g/cm3', 7.9)
    rpv_steel.add_nuclide("Fe54", 0.05437098, 'wo')
    rpv_steel.add_nuclide("Fe56", 0.88500663, 'wo')
    rpv_steel.add_nuclide("Fe57", 0.0208008, 'wo')
    rpv_steel.add_nuclide("Fe58", 0.00282159, 'wo')
    rpv_steel.add_nuclide("Ni58", 0.0067198, 'wo')
    rpv_steel.add_nuclide("Ni60", 0.0026776, 'wo')
    rpv_steel.add_nuclide("Mn55", 0.01, 'wo')
    rpv_steel.add_nuclide("Cr52", 0.002092475, 'wo')
    rpv_steel.add_nuclide("C0", 0.0025, 'wo')
    rpv_steel.add_nuclide("Cu63", 0.0013696, 'wo')

    lower_rad_ref = openmc.Material(6, name='Lower radial reflector')
    lower_rad_ref.set_density('g/cm3', 4.32)
    lower_rad_ref.add_nuclide("H1", 0.0095661, 'wo')
    lower_rad_ref.add_nuclide("O16", 0.0759107, 'wo')
    lower_rad_ref.add_nuclide("B10", 3.08409e-5, 'wo')
    lower_rad_ref.add_nuclide("B11", 1.40499e-4, 'wo')
    lower_rad_ref.add_nuclide("Fe54", 0.035620772088, 'wo')
    lower_rad_ref.add_nuclide("Fe56", 0.579805982228, 'wo')
    lower_rad_ref.add_nuclide("Fe57", 0.01362750048, 'wo')
    lower_rad_ref.add_nuclide("Fe58", 0.001848545204, 'wo')
    lower_rad_ref.add_nuclide("Ni58", 0.055298376566, 'wo')
    lower_rad_ref.add_nuclide("Mn55", 0.0182870, 'wo')
    lower_rad_ref.add_nuclide("Cr52", 0.145407678031, 'wo')
    lower_rad_ref.add_s_alpha_beta('c_H_in_H2O')

    upper_rad_ref = openmc.Material(7, name='Upper radial reflector / Top plate region')
    upper_rad_ref.set_density('g/cm3', 4.28)
    upper_rad_ref.add_nuclide("H1", 0.0086117, 'wo')
    upper_rad_ref.add_nuclide("O16", 0.0683369, 'wo')
    upper_rad_ref.add_nuclide("B10", 2.77638e-5, 'wo')
    upper_rad_ref.add_nuclide("B11", 1.26481e-4, 'wo')
    upper_rad_ref.add_nuclide("Fe54", 0.035953677186, 'wo')
    upper_rad_ref.add_nuclide("Fe56", 0.585224740891, 'wo')
    upper_rad_ref.add_nuclide("Fe57", 0.01375486056, 'wo')
    upper_rad_ref.add_nuclide("Fe58", 0.001865821363, 'wo')
    upper_rad_ref.add_nuclide("Ni58", 0.055815129186, 'wo')
    upper_rad_ref.add_nuclide("Mn55", 0.0184579, 'wo')
    upper_rad_ref.add_nuclide("Cr52", 0.146766614995, 'wo')
    upper_rad_ref.add_s_alpha_beta('c_H_in_H2O')

    bot_plate = openmc.Material(8, name='Bottom plate region')
    bot_plate.set_density('g/cm3', 7.184)
    bot_plate.add_nuclide("H1", 0.0011505, 'wo')
    bot_plate.add_nuclide("O16", 0.0091296, 'wo')
    bot_plate.add_nuclide("B10", 3.70915e-6, 'wo')
    bot_plate.add_nuclide("B11", 1.68974e-5, 'wo')
    bot_plate.add_nuclide("Fe54", 0.03855611055, 'wo')
    bot_plate.add_nuclide("Fe56", 0.627585036425, 'wo')
    bot_plate.add_nuclide("Fe57", 0.014750478, 'wo')
    bot_plate.add_nuclide("Fe58", 0.002000875025, 'wo')
    bot_plate.add_nuclide("Ni58", 0.059855207342, 'wo')
    bot_plate.add_nuclide("Mn55", 0.0197940, 'wo')
    bot_plate.add_nuclide("Cr52", 0.157390026871, 'wo')
    bot_plate.add_s_alpha_beta('c_H_in_H2O')

    bot_nozzle = openmc.Material(9, name='Bottom nozzle region')
    bot_nozzle.set_density('g/cm3', 2.53)
    bot_nozzle.add_nuclide("H1", 0.0245014, 'wo')
    bot_nozzle.add_nuclide("O16", 0.1944274, 'wo')
    bot_nozzle.add_nuclide("B10", 7.89917e-5, 'wo')
    bot_nozzle.add_nuclide("B11", 3.59854e-4, 'wo')
    bot_nozzle.add_nuclide("Fe54", 0.030411411144, 'wo')
    bot_nozzle.add_nuclide("Fe56", 0.495012237964, 'wo')
    bot_nozzle.add_nuclide("Fe57", 0.01163454624, 'wo')
    bot_nozzle.add_nuclide("Fe58", 0.001578204652, 'wo')
    bot_nozzle.add_nuclide("Ni58", 0.047211231662, 'wo')
    bot_nozzle.add_nuclide("Mn55", 0.0156126, 'wo')
    bot_nozzle.add_nuclide("Cr52", 0.124142524198, 'wo')
    bot_nozzle.add_s_alpha_beta('c_H_in_H2O')

    top_nozzle = openmc.Material(10, name='Top nozzle region')
    top_nozzle.set_density('g/cm3', 1.746)
    top_nozzle.add_nuclide("H1", 0.0358870, 'wo')
    top_nozzle.add_nuclide("O16", 0.2847761, 'wo')
    top_nozzle.add_nuclide("B10", 1.15699e-4, 'wo')
    top_nozzle.add_nuclide("B11", 5.27075e-4, 'wo')
    top_nozzle.add_nuclide("Fe54", 0.02644016154, 'wo')
    top_nozzle.add_nuclide("Fe56", 0.43037146399, 'wo')
    top_nozzle.add_nuclide("Fe57", 0.0101152584, 'wo')
    top_nozzle.add_nuclide("Fe58", 0.00137211607, 'wo')
    top_nozzle.add_nuclide("Ni58", 0.04104621835, 'wo')
    top_nozzle.add_nuclide("Mn55", 0.0135739, 'wo')
    top_nozzle.add_nuclide("Cr52", 0.107931450781, 'wo')
    top_nozzle.add_s_alpha_beta('c_H_in_H2O')

    top_fa = openmc.Material(11, name='Top of fuel assemblies')
    top_fa.set_density('g/cm3', 3.044)
    top_fa.add_nuclide("H1", 0.0162913, 'wo')
    top_fa.add_nuclide("O16", 0.1292776, 'wo')
    top_fa.add_nuclide("B10", 5.25228e-5, 'wo')
    top_fa.add_nuclide("B11", 2.39272e-4, 'wo')
    top_fa.add_nuclide("Zr90", 0.43313403903, 'wo')
    top_fa.add_nuclide("Zr91", 0.09549277374, 'wo')
    top_fa.add_nuclide("Zr92", 0.14759527104, 'wo')
    top_fa.add_nuclide("Zr94", 0.15280552077, 'wo')
    top_fa.add_nuclide("Zr96", 0.02511169542, 'wo')
    top_fa.add_s_alpha_beta('c_H_in_H2O')

    bot_fa = openmc.Material(12, name='Bottom of fuel assemblies')
    bot_fa.set_density('g/cm3', 1.762)
    bot_fa.add_nuclide("H1", 0.0292856, 'wo')
    bot_fa.add_nuclide("O16", 0.2323919, 'wo')
    bot_fa.add_nuclide("B10", 9.44159e-5, 'wo')
    bot_fa.add_nuclide("B11", 4.30120e-4, 'wo')
    bot_fa.add_nuclide("Zr90", 0.3741373658, 'wo')
    bot_fa.add_nuclide("Zr91", 0.0824858164, 'wo')
    bot_fa.add_nuclide("Zr92", 0.1274914944, 'wo')
    bot_fa.add_nuclide("Zr94", 0.1319920622, 'wo')
    bot_fa.add_nuclide("Zr96", 0.0216912612, 'wo')
    bot_fa.add_s_alpha_beta('c_H_in_H2O')

    # Define the materials file.
    model.materials = (fuel, clad, cold_water, hot_water, rpv_steel,
                       lower_rad_ref, upper_rad_ref, bot_plate,
                       bot_nozzle, top_nozzle, top_fa, bot_fa)

    ############################################################################################################
    # GEOMETRY

    # Define surfaces.
    s1 = openmc.ZCylinder(R=0.410, surface_id=1)
    s2 = openmc.ZCylinder(R=0.475, surface_id=2)
    s3 = openmc.ZCylinder(R=0.560, surface_id=3)
    s4 = openmc.ZCylinder(R=0.620, surface_id=4)
    s5 = openmc.ZCylinder(R=187.6, surface_id=5)
    s6 = openmc.ZCylinder(R=209.0, surface_id=6)
    s7 = openmc.ZCylinder(R=229.0, surface_id=7)
    s8 = openmc.ZCylinder(R=249.0, surface_id=8, boundary_type='vacuum')

    s31 = openmc.ZPlane(z0=-229.0, surface_id=31, boundary_type='vacuum')
    s32 = openmc.ZPlane(z0=-199.0, surface_id=32)
    s33 = openmc.ZPlane(z0=-193.0, surface_id=33)
    s34 = openmc.ZPlane(z0=-183.0, surface_id=34)
    s35 = openmc.ZPlane(z0=0.000, surface_id=35)
    s36 = openmc.ZPlane(z0=183.0, surface_id=36)
    s37 = openmc.ZPlane(z0=203.0, surface_id=37)
    s38 = openmc.ZPlane(z0=215.0, surface_id=38)
    s39 = openmc.ZPlane(z0=223.0, surface_id=39, boundary_type='vacuum')

    # Define pin cells.
    fuel_cold = openmc.Universe(name='Fuel pin, cladding, cold water', universe_id=1)
    c21 = openmc.Cell(cell_id=21, fill=fuel,       region=-s1)
    c22 = openmc.Cell(cell_id=22, fill=clad,       region=+s1 & -s2)
    c23 = openmc.Cell(cell_id=23, fill=cold_water, region=+s2)
    fuel_cold.add_cells((c21, c22, c23))

    tube_cold = openmc.Universe(name='Instrumentation guide tube, '
                                'cold water', universe_id=2)
    c24 = openmc.Cell(cell_id=24, fill=cold_water, region=-s3)
    c25 = openmc.Cell(cell_id=25, fill=clad,       region=+s3 & -s4)
    c26 = openmc.Cell(cell_id=26, fill=cold_water, region=+s4)
    tube_cold.add_cells((c24, c25, c26))

    fuel_hot = openmc.Universe(name='Fuel pin, cladding, hot water', universe_id=3)
    c27 = openmc.Cell(cell_id=27, fill=fuel,      region=-s1)
    c28 = openmc.Cell(cell_id=28, fill=clad,      region=+s1 & -s2)
    c29 = openmc.Cell(cell_id=29, fill=hot_water, region=+s2)
    fuel_hot.add_cells((c27, c28, c29))

    tube_hot = openmc.Universe(name='Instrumentation guide tube, hot water', universe_id=4)
    c30 = openmc.Cell(cell_id=30, fill=hot_water, region=-s3)
    c31 = openmc.Cell(cell_id=31, fill=clad,      region=+s3 & -s4)
    c32 = openmc.Cell(cell_id=32, fill=hot_water, region=+s4)
    tube_hot.add_cells((c30, c31, c32))

    # Set positions occupied by guide tubes
    tube_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8, 11, 14,
                       2, 5, 8, 11, 14, 3, 13, 5, 8, 11])
    tube_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8, 8,
                       11, 11, 11, 11, 11, 13, 13, 14, 14, 14])

    # Define fuel lattices.
    NUM_PINS = 17
    PITCH_PIN = 1.26

    l100 = openmc.RectLattice(name='Fuel assembly (lower half)', lattice_id=100)
    l100.lower_left = (-PITCH_PIN*NUM_PINS/2, -PITCH_PIN*NUM_PINS/2)
    l100.pitch = (PITCH_PIN, PITCH_PIN)
    l100.universes = np.tile(fuel_cold, (NUM_PINS, NUM_PINS))
    l100.universes[tube_x, tube_y] = tube_cold

    l101 = openmc.RectLattice(name='Fuel assembly (upper half)', lattice_id=101)
    l101.lower_left = (-PITCH_PIN*NUM_PINS/2, -PITCH_PIN*NUM_PINS/2)
    l101.pitch = (PITCH_PIN, PITCH_PIN)
    l101.universes = np.tile(fuel_hot, (NUM_PINS, NUM_PINS))
    l101.universes[tube_x, tube_y] = tube_hot

    # Define assemblies.
    fa_cw = openmc.Universe(name='Water assembly (cold)', universe_id=5)
    c50 = openmc.Cell(cell_id=50, fill=cold_water, region=+s34 & -s35)
    fa_cw.add_cell(c50)

    fa_hw = openmc.Universe(name='Water assembly (hot)', universe_id=7)
    c70 = openmc.Cell(cell_id=70, fill=hot_water, region=+s35 & -s36)
    fa_hw.add_cell(c70)

    fa_cold = openmc.Universe(name='Fuel assembly (cold)', universe_id=6)
    c60 = openmc.Cell(cell_id=60, fill=l100, region=+s34 & -s35)
    fa_cold.add_cell(c60)

    fa_hot = openmc.Universe(name='Fuel assembly (hot)', universe_id=8)
    c80 = openmc.Cell(cell_id=80, fill=l101, region=+s35 & -s36)
    fa_hot.add_cell(c80)

    # Define core lattices
    PITCH_ASSEMBLY = 21.42

    l200 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=200)
    l200.lower_left = (-PITCH_ASSEMBLY*21/2, -PITCH_ASSEMBLY*21/2)
    l200.pitch = (PITCH_ASSEMBLY, PITCH_ASSEMBLY)
    l200.universes = [
                    [fa_cw]*21,
                    [fa_cw]*21,
        [fa_cw]*7 + [fa_cold]*7  + [fa_cw]*7,
        [fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
        [fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
        [fa_cw]*2 + [fa_cold]*17 + [fa_cw]*2,
        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
        [fa_cw]*3 + [fa_cold]*15 + [fa_cw]*3,
        [fa_cw]*4 + [fa_cold]*13 + [fa_cw]*4,
        [fa_cw]*5 + [fa_cold]*11 + [fa_cw]*5,
        [fa_cw]*7 + [fa_cold]*7 + [fa_cw]*7,
                    [fa_cw]*21,
                    [fa_cw]*21]

    l201 = openmc.RectLattice(name='Core lattice (upper half)', lattice_id=201)
    l201.lower_left = (-PITCH_ASSEMBLY*21/2, -PITCH_ASSEMBLY*21/2)
    l201.pitch = (PITCH_ASSEMBLY, PITCH_ASSEMBLY)
    l201.universes = [
                    [fa_hw]*21,
                    [fa_hw]*21,
        [fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
        [fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
        [fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
        [fa_hw]*2 + [fa_hot]*17 + [fa_hw]*2,
        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
        [fa_hw]*3 + [fa_hot]*15 + [fa_hw]*3,
        [fa_hw]*4 + [fa_hot]*13 + [fa_hw]*4,
        [fa_hw]*5 + [fa_hot]*11 + [fa_hw]*5,
        [fa_hw]*7 + [fa_hot]*7 + [fa_hw]*7,
                    [fa_hw]*21,
                    [fa_hw]*21]

    # Define root universe.
    root = openmc.Universe(universe_id=0, name='root universe')
    c1 = openmc.Cell(cell_id=1, fill=l200, region=-s6 & +s34 & -s35)
    c2 = openmc.Cell(cell_id=2, fill=l201, region=-s6 & +s35 & -s36)
    c3 = openmc.Cell(cell_id=3, fill=bot_plate, region=-s7 & +s31 & -s32)
    c4 = openmc.Cell(cell_id=4, fill=bot_nozzle, region=-s5 & +s32 & -s33)
    c5 = openmc.Cell(cell_id=5, fill=bot_fa, region=-s5 & +s33 & -s34)
    c6 = openmc.Cell(cell_id=6, fill=top_fa, region=-s5 & +s36 & -s37)
    c7 = openmc.Cell(cell_id=7, fill=top_nozzle, region=-s5 & +s37 & -s38)
    c8 = openmc.Cell(cell_id=8, fill=upper_rad_ref, region=-s7 & +s38 & -s39)
    c9 = openmc.Cell(cell_id=9, fill=bot_nozzle, region=+s6 & -s7 & +s32 & -s38)
    c10 = openmc.Cell(cell_id=10, fill=rpv_steel, region=+s7 & -s8 & +s31 & -s39)
    c11 = openmc.Cell(cell_id=11, fill=lower_rad_ref, region=+s5 & -s6 & +s32 & -s34)
    c12 = openmc.Cell(cell_id=12, fill=upper_rad_ref, region=+s5 & -s6 & +s36 & -s38)
    root.add_cells((c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12))

    # Assign root universe to geometry
    model.geometry.root_universe = root

    ############################################################################################################
    # MESH

    mesh = openmc.Mesh(mesh_id=1)
    mesh.type = 'regular'
    mesh.dimension   = (500, 500, 4)
    mesh.lower_left  = (-229.0, -229.0, -229.0)
    mesh.upper_right = (229.0, 229.0, 229.0)

    # Create a mesh filter
    mesh_filter = openmc.MeshFilter(mesh)

    ############################################################################################################
    # CROSS-SECTION LIBRARY

    energy_groups = openmc.mgxs.EnergyGroups()
    energy_groups.group_edges = np.logspace(-3, 7.3, 21)

    # Instantiate a 1-group EnergyGroups object
    one_group = openmc.mgxs.EnergyGroups()
    one_group.group_edges = np.array([energy_groups.group_edges[0], energy_groups.group_edges[-1]])

    # Initialize an 20-energy-group and 6-delayed-group MGXS Library
    mgxs_lib = openmc.mgxs.Library(model.geometry)
    mgxs_lib.energy_groups = one_group
    mgxs_lib.num_delayed_groups = 6

    # Specify multi-group cross section types to compute
    mgxs_lib.mgxs_types = ['total', 'transport', 'nu-scatter matrix', 'kappa-fission', 'inverse-velocity', 'chi-prompt',
                           'prompt-nu-fission', 'chi-delayed', 'delayed-nu-fission', 'beta']
    # Specify a "mesh" domain type for the cross section tally filters
    mgxs_lib.domain_type = 'mesh'
    # Specify the mesh domain over which to compute multi-group cross sections
    mgxs_lib.domains = [mesh]

    # Construct all tallies needed for the multi-group cross section library
    mgxs_lib.build_library()

    ############################################################################################################
    # TALLIES

    # Create a "tallies.xml" file for the MGXS Library
    model.tallies = openmc.Tallies()
    mgxs_lib.add_to_tallies_file(model.tallies, merge=True)

    # Instantiate a flux tally; Other valid options: 'current', 'fission', etc
    flux_tally = openmc.Tally(name='Flux')
    flux_tally.filters = [mesh_filter]
    flux_tally.scores = ['flux']

    # Add tallies to the tallies file
    model.tallies.append(flux_tally)

    ############################################################################################################
    # SETTINGS

    model.settings.batches = 20
    model.settings.inactive = 5
    model.settings.particles = 1000
    uniform_dist = openmc.stats.Box([-160, -160, -183], [160, 160, 183], only_fissionable=False)
    model.settings.source = openmc.Source(space=uniform_dist)
    # to track particles 3 and 4 from batch 1 and generation 2:
    # model.settings.track = (1, 2, 3, 1, 2, 4)

    ############################################################################################################
    # PLOT

    plot = openmc.Plot().from_geometry(model.geometry)
    plot.filename = 'Full-Core'
    plot.pixels = (3000, 3000)
    plot.basis = 'xy'
    plot.color_by = 'material'
    model.plots.append(plot)

    openmc.plot_geometry(output=True)

    ############################################################################################################

    return model

pwr_core().export_to_xml()
openmc.run(output=True)

I then have a custom code which I am trying to use to retrieve axial slices of the tallied scores. Please take special note of the mesh configuration and filter in the model and see my custom parsing logic below:

def statepoint_evaluation(statepoint_file, desired_score=''):
    sp = openmc.StatePoint(filename=statepoint_file)

    score_df = sp.get_tally(scores=['flux']).get_slice(scores=['flux']).get_pandas_dataframe()

    axial_dfs = []
    for i in range(4):
        temp_df = score_df[score_df['mesh 1']['z'] == i+1]
        x = temp_df['mesh 1']['x']
        y = temp_df['mesh 1']['y']
        mean = temp_df['mean']

        new_df = pd.DataFrame(list(zip(x, y, mean)), columns=['x', 'y', 'mean'])
        new_df = new_df.pivot(index='y', columns='x', values='mean')

        axial_dfs.append(np.array(new_df.values))

    return axial_dfs

Oddly enough, what ends up happening is my plot "duplicates" itself in the y direction of the x-y plane for as many iterations as there are mesh elements in the z direction (mesh.dimension = (500, 500, 4) = 4) [see image below]

One could argue that maybe it is plotting logic but I will upload text output of axial_dfs and one should see that the data is in fact strange before it is plotted.

SterlingButters commented 5 years ago

SterlingButters commented 5 years ago

customparse.txt

paulromano commented 5 years ago

Hi @SterlingButters. OpenMC fully supports using 3D meshes and writing them to statepoint files. openmc-plot-mesh-tally should in theory be able to read data from a statepoint file and give you axial slices. The only thing that was deprecated/not working was the openmc-statepoint-3d script which specifically converted 3D mesh tally data to SILO/VTK files. That script has since been removed, although we are going to look at reimplementing that kind of capability. That being said, it's quite possible that there is a bug in openmc-plot-mesh-tally, so we'll have to look specifically at your case here to see if something is wrong.

@pshriwise would you be interested in looking into this?

SterlingButters commented 5 years ago

@paulromano Ahh, I see. Thank you and if there is any additional detail you need for reproduction of the issue let me know.

pshriwise commented 5 years ago

@paulromano @SterlingButters I'll take a crack at it and let you know what I find out.

pshriwise commented 5 years ago

@SterlingButters to make sure we're seeing the same things. What version of OpenMC are you using?

SterlingButters commented 5 years ago

Sterlings-MacBook-Pro:~ sterlingbutters$ openmc -v
 OpenMC version 0.9.0
 Copyright (c) 2011-2015 Massachusetts Institute of Technology
 MIT/X license at <http://openmc.readthedocs.io/en/latest/license.html>

SterlingButters commented 5 years ago

@pshriwise Were you ever able to triangulate the issue?

pshriwise commented 5 years ago

@SterlingButters apologies for letting this fall off my radar. I'm going to be looking into this soon.

SterlingButters commented 5 years ago

No problem, looking forward to it :)

SterlingButters commented 5 years ago

@pshriwise I thought these might help:

materials.xml:

<?xml version='1.0' encoding='utf-8'?>
<materials>
    <material depletable="true" id="1" name="UOX fuel">
        <density units="g/cm3" value="10.062" />
        <nuclide ao="4.9476e-06" name="U234" />
        <nuclide ao="0.00048218" name="U235" />
        <nuclide ao="0.021504" name="U238" />
        <nuclide ao="1.0801e-08" name="Xe135" />
        <nuclide ao="0.045737" name="O16" />
    </material>
    <material id="10" name="Top nozzle region">
        <density units="g/cm3" value="1.746" />
        <nuclide name="H1" wo="0.035887" />
        <nuclide name="O16" wo="0.2847761" />
        <nuclide name="B10" wo="0.000115699" />
        <nuclide name="B11" wo="0.000527075" />
        <nuclide name="Fe54" wo="0.02644016154" />
        <nuclide name="Fe56" wo="0.43037146399" />
        <nuclide name="Fe57" wo="0.0101152584" />
        <nuclide name="Fe58" wo="0.00137211607" />
        <nuclide name="Ni58" wo="0.04104621835" />
        <nuclide name="Mn55" wo="0.0135739" />
        <nuclide name="Cr52" wo="0.107931450781" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="11" name="Top of fuel assemblies">
        <density units="g/cm3" value="3.044" />
        <nuclide name="H1" wo="0.0162913" />
        <nuclide name="O16" wo="0.1292776" />
        <nuclide name="B10" wo="5.25228e-05" />
        <nuclide name="B11" wo="0.000239272" />
        <nuclide name="Zr90" wo="0.43313403903" />
        <nuclide name="Zr91" wo="0.09549277374" />
        <nuclide name="Zr92" wo="0.14759527104" />
        <nuclide name="Zr94" wo="0.15280552077" />
        <nuclide name="Zr96" wo="0.02511169542" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="12" name="Bottom of fuel assemblies">
        <density units="g/cm3" value="1.762" />
        <nuclide name="H1" wo="0.0292856" />
        <nuclide name="O16" wo="0.2323919" />
        <nuclide name="B10" wo="9.44159e-05" />
        <nuclide name="B11" wo="0.00043012" />
        <nuclide name="Zr90" wo="0.3741373658" />
        <nuclide name="Zr91" wo="0.0824858164" />
        <nuclide name="Zr92" wo="0.1274914944" />
        <nuclide name="Zr94" wo="0.1319920622" />
        <nuclide name="Zr96" wo="0.0216912612" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="2" name="Zircaloy">
        <density units="g/cm3" value="5.77" />
        <nuclide ao="0.5145" name="Zr90" />
        <nuclide ao="0.1122" name="Zr91" />
        <nuclide ao="0.1715" name="Zr92" />
        <nuclide ao="0.1738" name="Zr94" />
        <nuclide ao="0.028" name="Zr96" />
    </material>
    <material id="3" name="Cold borated water">
        <density units="atom/b-cm" value="0.07416" />
        <nuclide ao="2.0" name="H1" />
        <nuclide ao="1.0" name="O16" />
        <nuclide ao="0.000649" name="B10" />
        <nuclide ao="0.002689" name="B11" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="4" name="Hot borated water">
        <density units="atom/b-cm" value="0.06614" />
        <nuclide ao="2.0" name="H1" />
        <nuclide ao="1.0" name="O16" />
        <nuclide ao="0.000649" name="B10" />
        <nuclide ao="0.002689" name="B11" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="5" name="Reactor pressure vessel steel">
        <density units="g/cm3" value="7.9" />
        <nuclide name="Fe54" wo="0.05437098" />
        <nuclide name="Fe56" wo="0.88500663" />
        <nuclide name="Fe57" wo="0.0208008" />
        <nuclide name="Fe58" wo="0.00282159" />
        <nuclide name="Ni58" wo="0.0067198" />
        <nuclide name="Ni60" wo="0.0026776" />
        <nuclide name="Mn55" wo="0.01" />
        <nuclide name="Cr52" wo="0.002092475" />
        <nuclide name="C0" wo="0.0025" />
        <nuclide name="Cu63" wo="0.0013696" />
    </material>
    <material id="6" name="Lower radial reflector">
        <density units="g/cm3" value="4.32" />
        <nuclide name="H1" wo="0.0095661" />
        <nuclide name="O16" wo="0.0759107" />
        <nuclide name="B10" wo="3.08409e-05" />
        <nuclide name="B11" wo="0.000140499" />
        <nuclide name="Fe54" wo="0.035620772088" />
        <nuclide name="Fe56" wo="0.579805982228" />
        <nuclide name="Fe57" wo="0.01362750048" />
        <nuclide name="Fe58" wo="0.001848545204" />
        <nuclide name="Ni58" wo="0.055298376566" />
        <nuclide name="Mn55" wo="0.018287" />
        <nuclide name="Cr52" wo="0.145407678031" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="7" name="Upper radial reflector / Top plate region">
        <density units="g/cm3" value="4.28" />
        <nuclide name="H1" wo="0.0086117" />
        <nuclide name="O16" wo="0.0683369" />
        <nuclide name="B10" wo="2.77638e-05" />
        <nuclide name="B11" wo="0.000126481" />
        <nuclide name="Fe54" wo="0.035953677186" />
        <nuclide name="Fe56" wo="0.585224740891" />
        <nuclide name="Fe57" wo="0.01375486056" />
        <nuclide name="Fe58" wo="0.001865821363" />
        <nuclide name="Ni58" wo="0.055815129186" />
        <nuclide name="Mn55" wo="0.0184579" />
        <nuclide name="Cr52" wo="0.146766614995" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="8" name="Bottom plate region">
        <density units="g/cm3" value="7.184" />
        <nuclide name="H1" wo="0.0011505" />
        <nuclide name="O16" wo="0.0091296" />
        <nuclide name="B10" wo="3.70915e-06" />
        <nuclide name="B11" wo="1.68974e-05" />
        <nuclide name="Fe54" wo="0.03855611055" />
        <nuclide name="Fe56" wo="0.627585036425" />
        <nuclide name="Fe57" wo="0.014750478" />
        <nuclide name="Fe58" wo="0.002000875025" />
        <nuclide name="Ni58" wo="0.059855207342" />
        <nuclide name="Mn55" wo="0.019794" />
        <nuclide name="Cr52" wo="0.157390026871" />
        <sab name="c_H_in_H2O" />
    </material>
    <material id="9" name="Bottom nozzle region">
        <density units="g/cm3" value="2.53" />
        <nuclide name="H1" wo="0.0245014" />
        <nuclide name="O16" wo="0.1944274" />
        <nuclide name="B10" wo="7.89917e-05" />
        <nuclide name="B11" wo="0.000359854" />
        <nuclide name="Fe54" wo="0.030411411144" />
        <nuclide name="Fe56" wo="0.495012237964" />
        <nuclide name="Fe57" wo="0.01163454624" />
        <nuclide name="Fe58" wo="0.001578204652" />
        <nuclide name="Ni58" wo="0.047211231662" />
        <nuclide name="Mn55" wo="0.0156126" />
        <nuclide name="Cr52" wo="0.124142524198" />
        <sab name="c_H_in_H2O" />
    </material>
</materials>

geometry.xml:

<?xml version='1.0' encoding='utf-8'?>
<geometry>
    <cell fill="200" id="1" region="-6 34 -35" universe="0" />
    <cell id="10" material="5" region="7 -8 31 -39" universe="0" />
    <cell id="11" material="6" region="5 -6 32 -34" universe="0" />
    <cell id="12" material="7" region="5 -6 36 -38" universe="0" />
    <cell fill="201" id="2" region="-6 35 -36" universe="0" />
    <cell id="21" material="1" region="-1" universe="1" />
    <cell id="22" material="2" region="1 -2" universe="1" />
    <cell id="23" material="3" region="2" universe="1" />
    <cell id="24" material="3" region="-3" universe="2" />
    <cell id="25" material="2" region="3 -4" universe="2" />
    <cell id="26" material="3" region="4" universe="2" />
    <cell id="27" material="1" region="-1" universe="3" />
    <cell id="28" material="2" region="1 -2" universe="3" />
    <cell id="29" material="4" region="2" universe="3" />
    <cell id="3" material="8" region="-7 31 -32" universe="0" />
    <cell id="30" material="4" region="-3" universe="4" />
    <cell id="31" material="2" region="3 -4" universe="4" />
    <cell id="32" material="4" region="4" universe="4" />
    <cell id="4" material="9" region="-5 32 -33" universe="0" />
    <cell id="5" material="12" region="-5 33 -34" universe="0" />
    <cell id="50" material="3" region="34 -35" universe="5" />
    <cell id="6" material="11" region="-5 36 -37" universe="0" />
    <cell fill="100" id="60" region="34 -35" universe="6" />
    <cell id="7" material="10" region="-5 37 -38" universe="0" />
    <cell id="70" material="4" region="35 -36" universe="7" />
    <cell id="8" material="7" region="-7 38 -39" universe="0" />
    <cell fill="101" id="80" region="35 -36" universe="8" />
    <cell id="9" material="9" region="6 -7 32 -38" universe="0" />
    <lattice id="100" name="Fuel assembly (lower half)">
        <pitch>1.26 1.26</pitch>
        <dimension>17 17</dimension>
        <lower_left>-10.71 -10.71</lower_left>
        <universes>
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 
1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 </universes>
    </lattice>
    <lattice id="101" name="Fuel assembly (upper half)">
        <pitch>1.26 1.26</pitch>
        <dimension>17 17</dimension>
        <lower_left>-10.71 -10.71</lower_left>
        <universes>
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 
3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 </universes>
    </lattice>
    <lattice id="200" name="Core lattice (lower half)">
        <pitch>21.42 21.42</pitch>
        <dimension>21 21</dimension>
        <lower_left>-224.91 -224.91</lower_left>
        <universes>
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 
5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 
5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 
5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 
5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 
5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 </universes>
    </lattice>
    <lattice id="201" name="Core lattice (upper half)">
        <pitch>21.42 21.42</pitch>
        <dimension>21 21</dimension>
        <lower_left>-224.91 -224.91</lower_left>
        <universes>
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 
7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 
7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 
7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 
7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 
7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 </universes>
    </lattice>
    <surface coeffs="0.0 0.0 0.41" id="1" type="z-cylinder" />
    <surface coeffs="0.0 0.0 0.475" id="2" type="z-cylinder" />
    <surface coeffs="0.0 0.0 0.56" id="3" type="z-cylinder" />
    <surface boundary="vacuum" coeffs="-229.0" id="31" type="z-plane" />
    <surface coeffs="-199.0" id="32" type="z-plane" />
    <surface coeffs="-193.0" id="33" type="z-plane" />
    <surface coeffs="-183.0" id="34" type="z-plane" />
    <surface coeffs="0.0" id="35" type="z-plane" />
    <surface coeffs="183.0" id="36" type="z-plane" />
    <surface coeffs="203.0" id="37" type="z-plane" />
    <surface coeffs="215.0" id="38" type="z-plane" />
    <surface boundary="vacuum" coeffs="223.0" id="39" type="z-plane" />
    <surface coeffs="0.0 0.0 0.62" id="4" type="z-cylinder" />
    <surface coeffs="0.0 0.0 187.6" id="5" type="z-cylinder" />
    <surface coeffs="0.0 0.0 209.0" id="6" type="z-cylinder" />
    <surface coeffs="0.0 0.0 229.0" id="7" type="z-cylinder" />
    <surface boundary="vacuum" coeffs="0.0 0.0 249.0" id="8" type="z-cylinder" />
</geometry>

pshriwise commented 5 years ago

Finally had a chance to dig into this today. Thanks for your patience @SterlingButters!

In version 0.9.0 there was an incorrect assumption about regular mesh index ordering in statepoint files (xyz vs. zyx). The problematic method can be seen here in filter.py's get_bin_index method.

After updating this method to

    def get_bin_index(self, filter_bin):                                                                                                                                                                                                      
        if len(self.mesh.dimension) == 3:                                                                                                                                                                                                      
            nx, ny, nz = self.mesh.dimension                                                                                                                                                                                                   
            i, j, k = filter_bin                                                                                                                                                                                                               
            val = ((k-1) * ny + (j-1)) * nx + (i-1)                                                                                                                                                                                            
        else:                                                                                                                                                                                                                                  
            nx, ny = self.mesh.dimension                                                                                                                                                                                                       
            i, j, k = filter_bin                                                                                                                                                                                                               
            val = (j-1) * nx + (i-1)                                                                                                                                                                                                           

        return val

the repeated structures problem seems to be resolved. If you're still on 0.9.0 and want to try this fix out to verify please feel free!

This issue has been resolved in develop by using a NumPy mask on the full list of indices the MeshFilter occupies and when re-running this problem using that branch I had similar success viewing the mesh tally data from this problem.

I'll apply this fix to a problem I have a better intuition for as well to make sure the fix isn't displaying misleading data and double-check it against develop.

SterlingButters commented 5 years ago

Just to clarify, this has since been fixed in 0.10.0?

pshriwise commented 5 years ago

It looks like this has been updated in our develop branch, but not in versions 0.9.0 or 0.10.0. If you're restricted to working with one of those versions, we can hopefully help you work with an older statepoint file.

SterlingButters commented 5 years ago

That's ok, I'll just edit filter.py - will post back in several hours with confirmation of functionality

SterlingButters commented 5 years ago

@pshriwise Was

    def get_bin_index(self, filter_bin):
        # Filter bins for a mesh are an (x,y,z) tuple. Convert (x,y,z) to a
        # single bin -- this is similar to subroutine mesh_indices_to_bin in
        # openmc/src/mesh.F90.
        if len(self.mesh.dimension) == 3:
            nx, ny, nz = self.mesh.dimension
            val = (filter_bin[0] - 1) * ny * nz + \
                  (filter_bin[1] - 1) * nz + \
                  (filter_bin[2] - 1)
        else:
            nx, ny = self.mesh.dimension
            val = (filter_bin[0] - 1) * ny + \
                  (filter_bin[1] - 1)

to

    def get_bin_index(self, filter_bin):                                                                                                                                                                                                      
        if len(self.mesh.dimension) == 3:                                                                                                                                                                                                      
            nx, ny, nz = self.mesh.dimension                                                                                                                                                                                                   
            i, j, k = filter_bin                                                                                                                                                                                                               
            val = ((k-1) * ny + (j-1)) * nx + (i-1)                                                                                                                                                                                            
        else:                                                                                                                                                                                                                                  
            nx, ny = self.mesh.dimension                                                                                                                                                                                                       
            i, j, k = filter_bin                                                                                                                                                                                                               
            val = (j-1) * nx + (i-1)                                                                                                                                                                                                           

        return val

the only change you made? I am still getting the repeated structures:

Is it possible develop contains other changes responsible for the fix you achieved?

pshriwise commented 5 years ago

Strange! Yes that was the only change I made. What was the location of the filter.py file you modified if you don't mind my asking?

SterlingButters commented 5 years ago

@pshriwise

(base) Sterlings-MBP:~ sterlingbutters$ cd /anaconda3/lib/python3.7/site-packages/openmc
(base) Sterlings-MBP:openmc sterlingbutters$ ls
__init__.py    examples.py  mixin.py           source.py
__pycache__    executor.py  model              statepoint.py
arithmetic.py  filter.py    nuclide.py         stats
capi           geometry.py  openmoc_compatible.py  summary.py
cell.py        lattice.py   particle_restart.py    surface.py
checkvalue.py  macroscopic.py   plots.py           tallies.py
clean_xml.py   material.py  plotter.py         tally_derivative.py
cmfd.py        mesh.py      region.py          trigger.py
data           mgxs     search.py          universe.py
element.py     mgxs_library.py  settings.py        volume.py

openmc-dev / openmc

Reopening of Issue #971: 3D Tally-Score Capability #1106