get UnstructuredMesh centroids and volumes without loading from statepoint

[shimwell]

Description

It would be useful if we can get mesh information such as volumes and centroids of the unstructured mesh. Currently this is possible if one loads in an unstructured mesh from a statepoint file but not if an UnstructuredMesh is instantiated in python like this

 umesh = openmc.UnstructuredMesh(filename="umesh.h5m",library='moab')
 umesh.centroids
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/j/openmc/openmc/mesh.py", line 2066, in centroids
    return np.array([self.centroid(i) for i in range(self.n_elements)])
  File "/home/j/openmc/openmc/mesh.py", line 2071, in n_elements
    raise RuntimeError("No information about this mesh has "
RuntimeError: No information about this mesh has been loaded from a statepoint file.

This would be useful as knowing the centroids or volumes could help with defining a openmc.MeshSource. It would also bring the UnstructuredMesh behaviour more inline with the other meshes such as RegularMesh

Alternatives

We write the UnstructuredMesh to a statepoint file, then read in the statepoint file and get the UnstructuredMesh from hdf5 then it contains centroids and volume information

Compatibility

this would maintain the same api usage

[pshriwise]

A good thought. I'm sure you can imagine why this is difficult from the Python API as the libraries we rely on for unstructured mesh (MOAB and libMesh) interact with OpennMC through a C++ interface. The mesh volumes are already exposed through the CAPI in the openmc.lib module. If the centroids were also to be added would that suit your needs?

A valid OpenMC problem would still need to be defined and loaded before determining your mesh source strengths, but maybe it's still useful?

[shimwell]

latest attempt

import openmc
import openmc.lib

openmc.config['cross_sections'] = '/home/j/endf-b8.0-hdf5/endfb-viii.0-hdf5/cross_sections.xml'

umesh = openmc.UnstructuredMesh(filename="umesh.h5m",library='moab')

surf1 = openmc.Sphere(r=50000, boundary_type="vacuum")
region1 = -surf1

cell1 = openmc.Cell(region=region1)

my_geometry = openmc.Geometry([cell1])

my_source = openmc.IndependentSource()
my_source.space = openmc.stats.Point((0, 0, 0))
my_source.angle = openmc.stats.Isotropic()
my_source.energy = openmc.stats.Discrete([14e6], [1])

all_sources = [my_source] * 308
mesh_source = openmc.MeshSource(
    mesh=umesh,
    sources=all_sources,
)

my_settings = openmc.Settings()
my_settings.batches = 1
my_settings.particles = 1
my_settings.run_mode = "fixed source"
my_settings.source = mesh_source

model = openmc.model.Model(my_geometry, None, my_settings )
model.export_to_model_xml()

openmc.lib.init()
openmc.lib.UnstructuredMesh()

[pshriwise]

At the end, perhaps try: unstructured_mesh = openmc.lib.meshes[umesh.id]

[shimwell]

Thanks Patrick

It is a bit tricky with the need to include the mesh in the model so that openmc.lib can find it.

Originally I tried to include it in the model via a MeshSource but that didn't work and raised the error

"/home/j/openmc/openmc/source.py", line 503, in populate_xml_element
    for idx in self.mesh.indices:
AttributeError: 'UnstructuredMesh' object has no attribute 'indices'

umesh = openmc.UnstructuredMesh(filename="umesh.h5m",library='moab')
mesh_source = openmc.MeshSource(
    mesh=umesh,
    sources=all_sources,
)
settings.source = mesh_source # this line raises the error

So I shall include the mesh via a tally and then I should be able to get the fully populated unstrucutred mesh object and then I should be able to make a meshSource from it.

Adding centroids would be great. I think longer term we might want to add some extra code to meshSource or settings.source to simplify the process of having a meshsource with and unstrucutred mesh

[shimwell]

update. I tried reading the unstructured mesh in from the statepoint and did get populated centroids and volumes. However it appears that an unstrucutred mesh can't be used as a mesh for an openmc.MeshSource. When trying to use a umesh read in from the statepoint we get the AttributeError

model_gamma.run(cwd="photons")
  File "/home/j/openmc/openmc/model/model.py", line 712, in run
    self.export_to_model_xml(**export_kwargs)
  File "/home/j/openmc/openmc/model/model.py", line 509, in export_to_model_xml
    settings_element = self.settings.to_xml_element(mesh_memo)
  File "/home/j/openmc/openmc/settings.py", line 1791, in to_xml_element
    self._create_source_subelement(element, mesh_memo)
  File "/home/j/openmc/openmc/settings.py", line 1076, in _create_source_subelement
    root.append(source.to_xml_element())
  File "/home/j/openmc/openmc/source.py", line 77, in to_xml_element
    self.populate_xml_element(element)
  File "/home/j/openmc/openmc/source.py", line 503, in populate_xml_element
    for idx in self.mesh.indices:
AttributeError: 'UnstructuredMesh' object has no attribute 'indices'

I tried to add that attribute with
umesh_from_sp.indices = [i+1 for i in range(len(mesh_vols))]

but this just raised a type error when trying to write the xml

model_gamma.run(cwd="photons")
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/j/openmc/openmc/model/model.py", line 712, in run
    self.export_to_model_xml(**export_kwargs)
  File "/home/j/openmc/openmc/model/model.py", line 509, in export_to_model_xml
    settings_element = self.settings.to_xml_element(mesh_memo)
  File "/home/j/openmc/openmc/settings.py", line 1791, in to_xml_element
    self._create_source_subelement(element, mesh_memo)
  File "/home/j/openmc/openmc/settings.py", line 1076, in _create_source_subelement
    root.append(source.to_xml_element())
  File "/home/j/openmc/openmc/source.py", line 77, in to_xml_element
    self.populate_xml_element(element)
  File "/home/j/openmc/openmc/source.py", line 504, in populate_xml_element
    idx = tuple(i - 1 for i in idx)
TypeError: 'int' object is not iterable

I have a not so minimal example over here if it helps

[pshriwise]

Glad to hear that reading the volumes and centroids went alright @shimwell!

It seems we are in fact missing the indices property for the UnstructuredMesh class. This should really be present on the MeshBase class as an abstract method to protect from this as we expect all mesh implementations to have it. Should be an easy fix and I can update it shortly.

The reason that the population of the indices attribute you have above isn't working is because that property is expected to be a tuple of indices. See https://docs.openmc.org/en/latest/pythonapi/generated/openmc.RegularMesh.html?highlight=meshbase#openmc-regularmesh

Hopefully that helps you move forward in the meantime.

[shimwell]

Thanks Patrick and Paul (on slack), I appreciate your help going through some of these less trodden workflows.

I changed

umesh_from_sp.indices = [i+1 for i in range(len(mesh_vols))]

to

umesh_from_sp.indices = [(i+1, i+1, i+1) for i in range(len(mesh_vols))]

in my example and it now writes the model.xml file, so it does get past that indices check. The code still fails on the run command as now later in the code it recieves a tuple with three values when unstructured meshes have 1D indexing.

    gamma_sp_filename = model_gamma.run(cwd="photons")
  File "/home/j/openmc/openmc/model/model.py", line 712, in run
    self.export_to_model_xml(**export_kwargs)
  File "/home/j/openmc/openmc/model/model.py", line 509, in export_to_model_xml
    settings_element = self.settings.to_xml_element(mesh_memo)
  File "/home/j/openmc/openmc/settings.py", line 1791, in to_xml_element
    self._create_source_subelement(element, mesh_memo)
  File "/home/j/openmc/openmc/settings.py", line 1076, in _create_source_subelement
    root.append(source.to_xml_element())
  File "/home/j/openmc/openmc/source.py", line 77, in to_xml_element
    self.populate_xml_element(element)
  File "/home/j/openmc/openmc/source.py", line 505, in populate_xml_element
    elem.append(self.sources[idx].to_xml_element())
IndexError: too many indices for array: array is 1-dimensional, but 3 were indexed
>>> 

So I guess the indices could be moved to MeshBase and used by regular mesh, spherical mesh etc but it the UnstrucutreMesh class might need its own implementation as it's indices appear to be 1D in places

[shimwell]

Would it be possible to keep this issue open till that add_library_data method that loads in the centroids and volumes using openmc.lb is possible

[pshriwise]

Yep yep. Hopefully I'll be able to close that out soon!