Nuclear data dependent zero result for heating-local tally

User Support
1

I have a rather strange issue, which is that the tally for ‘heating-local’ using a mesh filter is giving either non-zero or identically zero results depending on nuclear data I use. When I use the nuclear data cross sections required for the python test suite, nndc_hdf5, I get the zero results. If however I use endfb71_hdf5 I get non-trivial results. Below is a short quick script that I find I can reproduce these effects. I’d really appreciate knowing if anyone else can reproduce this issue, and then what the reason is.

#!/usr/bin/env python3

import openmc
import os
import subprocess
import matplotlib.pyplot as plt

def createMaterials():

    copper = openmc.Material(1, "copper")
    copper.add_element('Cu', 1.0)
    copper.set_density('g/cm3', 8.96)
    copper.temperature = 300
    
    air = openmc.Material(name='air')
    air.set_density('g/cc', 0.001205)
    air.add_element('N', 0.784431)
    air.add_element('O', 0.210748)
    air.add_element('Ar',0.0046)
    
    mats = openmc.Materials([copper, air])
    mats.export_to_xml()
    print("Created materials.xml")

    materialsDict={}
    for mat in mats:
        materialsDict[mat.name]=mat
        print(mat)

    return materialsDict

def createPlots( materials ):

    copper = materials['copper']
    air = materials['air']

    p_xy = openmc.Plot()
    p_xy.width = (25.0, 25.0)
    p_xy.pixels = (400, 400)
    p_xy.origin = (0.,0.,0.)
    p_xy.color_by = 'material'
    p_xy.colors = {copper: 'orange', air: 'blue'}
    p_xy.basis = 'xy'

    p_xz = openmc.Plot()
    p_xz.width = (25.0, 25.0)
    p_xz.pixels = (400, 400)
    p_xz.origin = (0.,0.,0.)
    p_xz.color_by = 'material'
    p_xz.colors = {copper: 'orange', air: 'blue'}
    p_xz.basis = 'xz'

    p_yz = openmc.Plot()
    p_yz.width = (25.0, 25.0)
    p_yz.pixels = (400, 400)
    p_yz.origin = (5.,0.,0.)
    p_yz.color_by = 'material'
    p_yz.colors = {copper: 'orange', air: 'blue'}
    p_yz.basis = 'yz'
    
    plots = openmc.Plots([p_xy,p_xz,p_yz])
    plots.export_to_xml()

def createGeometry( materials):

    len=20
    left = openmc.XPlane(x0=-len/2, boundary_type='vacuum')
    right = openmc.XPlane(x0=len/2, boundary_type='vacuum')
    bottom = openmc.YPlane(y0=-len/2, boundary_type='vacuum')
    top = openmc.YPlane(y0=len/2, boundary_type='vacuum')
    back = openmc.ZPlane(z0=-len/2, boundary_type='vacuum')
    front = openmc.ZPlane(z0=len/2, boundary_type='vacuum')

    len=5
    offset=5
    left2 = openmc.XPlane(x0=-len/2+offset)
    right2 = openmc.XPlane(x0=len/2+offset)
    bottom2 = openmc.YPlane(y0=-len/2)
    top2 = openmc.YPlane(y0=len/2)
    back2 = openmc.ZPlane(z0=-len/2)
    front2 = openmc.ZPlane(z0=len/2)
    

    big_box = +left & -right & +bottom & -top & +back & -front
    copper_region = +left2 & -right2 & +bottom2 & -top2 & +back2 & -front2
    air_region = big_box | ~copper_region

    copper = openmc.Cell(1, 'copper')
    copper.fill = materials["copper"]
    copper.region = copper_region

    air = openmc.Cell(2, 'air')
    air.fill = materials["air"]
    air.region = air_region
    
    root = openmc.Universe(cells=(copper,air))
    geom = openmc.Geometry(root)
    geom.export_to_xml()
    print("Created geometry.xml")
    
def createSettings():

    # Create neutron source
    point = openmc.stats.Point((0, 0, 0))
    src = openmc.Source(space=point)
    src.particle='neutron'

    # Create settings object
    settings = openmc.Settings()
    settings.source = src
    settings.run_mode = 'fixed source'
    settings.photon_transport = False
    
    # Turn on dagmc
    settings.batches = 10
    settings.inactive = 2
    settings.particles = 5000
    settings.export_to_xml()
    print("Created settings.xml")

def createTallies():

    # Filters
    
    # Create mesh filter for tally
    mesh = openmc.RegularMesh()
    mesh.dimension = [50, 50]
    mesh.lower_left = [-10,-10]
    mesh.upper_right = [10,10]
    mesh_filter = openmc.MeshFilter(mesh)

    # Create mesh tally to score flux and fission rate
    tally = openmc.Tally()
    tally.filters = [mesh_filter]
    tally.scores = ['heating-local', 'flux']
    tally.estimator = 'tracklength'    
    tallies = openmc.Tallies([tally])
    tallies.export_to_xml()
    print("Created tallies.xml")

def setupInput():    
    mats = createMaterials()
    createGeometry(mats)
    createPlots(mats)
    createSettings()
    createTallies()

setupInput()
#openmc.plot_geometry()
openmc.run()

with openmc.StatePoint("statepoint.10.h5") as sp:
    tally = sp.tallies[1]    
    flux = tally.get_slice(scores=['flux'])
    heating = tally.get_slice(scores=['heating-local'])

    flux.mean.shape = (50, 50)
    heating.mean.shape = (50, 50)

    fig = plt.subplot(111)
    fig.imshow(flux.mean)
    plt.show()
    
    fig2 = plt.subplot(111)
    fig2.imshow(heating.mean)
    plt.show()
2

Welcome to the forum @helen-brooks. heating-local tallies rely on a specially-processed set of KERMA cross sections. Namely, the KERMAs are calculated assuming that photon energy is deposited locally. To obtain these KERMAs requires running NJOY/HEATR with a special flag set. When we are generating libraries based on ACE files that have already been generated by someone else, NJOY has already been run and thus we have no ability to add the special KERMA needed for heating-local. However, when we are generating libraries starting from ENDF files (using openmc.data.IncidentNeutron.from_njoy), we do have the opportunity to run NJOY with whatever run sequence is needed.

What this means is that the “official” libraries that we distribute on openmc.org have been processed by our team and do have the data needed for heating-local. All the libraries that have been converted from other sources, e.g. the libraries based on MCNP data, do not have the data needed for heating-local and will thus produce a zero value in the tally. We ought to have a warning/error about this so that it doesn’t come as a surprise to a user.