Skip to content

XDGMesh Object#3889

Open
pshriwise wants to merge 27 commits into
openmc-dev:developfrom
pshriwise:xdg-mesh
Open

XDGMesh Object#3889
pshriwise wants to merge 27 commits into
openmc-dev:developfrom
pshriwise:xdg-mesh

Conversation

@pshriwise
Copy link
Copy Markdown
Contributor

Description

This introduces an initial integration of the XDG libarary for support of CAD-based mesh geometry in OpenMC. This code change provides tally support for volumetric mesh.

Capabilities

Previously, a tracklength estimator was not available for libMesh meshes due to the lack of support for surface-to-surface raytracing natively. XDG adds this support for libMesh meshes to enable raytracing for meshes that partially overlap the domain or have gaps.

Tetrahedral Mesh MOAB LibMesh XDG (libMesh or MOAB)
Tracklength Estimator
Collision Estimator
Hexahedral Mesh (coming soon) MOAB LibMesh XDG (libMesh or MOAB)
Tracklength Estimator
Collision Estimator

Algorithmic Changes

Previously, the MOAB tracklength estimator employed an algorithm that first determined all of the segments lying within elements for a given track and then performed a point location call at the midpoint of each segment. This algorithm is robust but computationally inefficient.

image

In XDG, a more standard element adjacency walk has been implemented requiring only one point location call per track. In the case that a track exits the mesh as elements are traversed, surface-to-surface ray tracing is employed to determine the track's re-entry location, if any.

image

This algorithmic change brings us up to date with standard ways of approaching this problem -- to great effect!

Results

Simple problems

The tally capabilities of this library have been tested on several simple problems (benchmark problems as well, but we'll save those for a later date). Below is a set of results using XDG for a superimposed mesh tally (CSG geometry) in OpenMC. A significant improvement is seen. The improved threaded scalability comes from the transition to a modern raytracing kernel (Embree) and, in the case of MOAB, bypassing most of the standard interface calls for element connectivity, coordinates, and adjacency. Serial improvements are the outcome of algorithmic changes to the tally routine (walking elements instead of many point location calls). The combination of these result in massive speedups in a problem like this where most of the cost is in the tally kernel computing element crossings (up to 300x). A per-element result comparison of the XDG tally to the MOAB tally using a tracklength estimator is shown below. The distribution of the differences is a little unusual, possibly indicating something systematic, but the absolute scale of the differences is quite small.

image image image image

Production problem: WISTEL-D

@connormoreno was kind enough to provide us with a WISTEL-D model to test this out on.

This model uses meshes in a number of ways:

  • The geometry is a DAGMC surface mesh (here we're using double-down as well for the raytracer)
  • A mesh source is used to sample neutron locations within the plasma volume. This is supported by XDG's find_element point locatoin capability.
  • An unstructured mesh heating tally is included for heating in the first wall.

It's noted in the figure, but the maximum difference between the tally results of a MOAB tracklength estimator (current capability) and the XDG mesh (also using MOAB but with the improvements noted above) is extremely small and the performance improvement is significant. This simulation ran 300M particles in total.

image image

Testing

The existing unstructured mesh tests for the direct MOAB/libMesh implementations have been adapted for use in XDG with the same measures of tally equivalence using tetrahedral mesh containing groups of elements that conform to the boundaries of a regular mesh.

Checklist

  • I have performed a self-review of my own code
  • I have run clang-format (version 18) on any C++ source files (if applicable)
  • I have followed the style guidelines for Python source files (if applicable)
  • I have made corresponding changes to the documentation (if applicable)
  • I have added tests that prove my fix is effective or that my feature works (if applicable)

@pshriwise pshriwise requested a review from paulromano as a code owner March 18, 2026 20:41
@nuclearkevin
Copy link
Copy Markdown
Member

Just stopping by to say this is awesome! Looking forward to integrating this into Cardinal and testing it out on some tally AMR problems

Comment thread src/xdg.cpp Outdated
Comment thread src/initialize.cpp Outdated
Co-authored-by: Kevin Sawatzky <66632997+nuclearkevin@users.noreply.github.com>
@shimwell
Copy link
Copy Markdown
Member

I had a question about the Python API design. Currently users create unstructured meshes like this:

umesh = openmc.UnstructuredMesh("mesh.vtk", library="moab")
umesh = openmc.UnstructuredMesh("mesh.exo", library="libmesh")                                              

Since XDGMesh inherits from UnstructuredMesh on the C++ side, would it make sense to follow the same pattern rather than introducing a separate class? For example:

umesh = openmc.UnstructuredMesh("mesh.e", library="xdg-moab")
umesh = openmc.UnstructuredMesh("mesh.e", library="xdg-libmesh")

Instead of:

from openmc.xdg import XDGMesh         
xdg_mesh = XDGMesh("mesh.e", library="moab")

From a user perspective XDG feels like another backend for unstructured meshes rather than a fundamentally different concept, so keeping it under the same class might make the API easier to discover and reduce some of the duplication. Happy to hear your thoughts if there's a reason to keep them separate though am I missing something something

Comment thread openmc/xdg.py
Comment on lines +26 to +28
Location of the unstructured mesh file. Supported files for 'moab'
library are .h5 and .vtk. Supported files for 'libmesh' library are
exodus mesh files .exo.
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a couple of questions here.

Looking at the tests it looks like supported files for 'moab' library are .h5 and .e (which is also exodus format?). I was happy when I saw .vtk in the docstring as a supported file type can we please have a vtk file in the test suit as well.

Would it be good to add .e to the moab supported files, does moab allow .exo? if not then I guess we want a comment to say .e is exodus format.

Suggested change
Location of the unstructured mesh file. Supported files for 'moab'
library are .h5 and .vtk. Supported files for 'libmesh' library are
exodus mesh files .exo.
Location of the unstructured mesh file. Supported files for 'moab'
library are .h5, .vtk and .e (exodus). Supported files for 'libmesh'
library are exodus mesh files .exo.

jon-proximafusion added a commit to shimwell/openmc that referenced this pull request Apr 13, 2026
Combines the wheel-building infrastructure from making-wheel-3 with XDG
mesh support from PR openmc-dev#3889. Adds XDG (with MOAB and libMesh) build steps
to manylinux.Dockerfile and enables OPENMC_USE_XDG in the wheel build.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
jon-proximafusion added a commit to shimwell/openmc that referenced this pull request Apr 13, 2026
Combines the wheel-building infrastructure from making-wheel-3 with XDG
mesh support from PR openmc-dev#3889. Adds XDG (with MOAB and libMesh) build steps
to manylinux.Dockerfile and enables OPENMC_USE_XDG in the wheel build.
@pshriwise
Copy link
Copy Markdown
Contributor Author

From a user perspective XDG feels like another backend for unstructured meshes rather than a fundamentally different concept, so keeping it under the same class might make the API easier to discover and reduce some of the duplication. Happy to hear your thoughts if there's a reason to keep them separate though am I missing something something

You make a good point here. XDG really does serve as another backend in this context. There are places where I'd like to differentiate XDG-based functionality into new objects (i.e. universes) , but this doesn't need to be one of them.

@jon-proximafusion
Copy link
Copy Markdown
Contributor

If anyone else out there like me wants to test this feature then I have built some handy wheels with dagmc (also with embree and dd) and xdg and they are downloadable as artifacts from this action run.

linux
https://github.com/shimwell/openmc/actions/runs/24345620993

mac
https://github.com/shimwell/openmc/actions/runs/24345621022

@eepeterson
Copy link
Copy Markdown
Contributor

From a user perspective XDG feels like another backend for unstructured meshes rather than a fundamentally different concept, so keeping it under the same class might make the API easier to discover and reduce some of the duplication. Happy to hear your thoughts if there's a reason to keep them separate though am I missing something something

You make a good point here. XDG really does serve as another backend in this context. There are places where I'd like to differentiate XDG-based functionality into new objects (i.e. universes) , but this doesn't need to be one of them.

I'd disagree with this interpretation since XDG is not really a mesh backend it's an API that OpenMC can leverage to have consistent mesh-based operations across many actual backends. If XDG is the touch point between OpenMC and any external mesh library (MOAB, LibMesh, MFEM, etc) it keeps the OpenMC code and build options much cleaner.

@pshriwise
Copy link
Copy Markdown
Contributor Author

From a user perspective XDG feels like another backend for unstructured meshes rather than a fundamentally different concept, so keeping it under the same class might make the API easier to discover and reduce some of the duplication. Happy to hear your thoughts if there's a reason to keep them separate though am I missing something something

You make a good point here. XDG really does serve as another backend in this context. There are places where I'd like to differentiate XDG-based functionality into new objects (i.e. universes) , but this doesn't need to be one of them.

I'd disagree with this interpretation since XDG is not really a mesh backend it's an API that OpenMC can leverage to have consistent mesh-based operations across many actual backends. If XDG is the touch point between OpenMC and any external mesh library (MOAB, LibMesh, MFEM, etc) it keeps the OpenMC code and build options much cleaner.

Also a good point @eepeterson. And I tend to agree for reasons in addition to multiple mesh backends. This PR introduces XDG as a capability purely for unstructured mesh tallies, but in subsequent PRs I'll be introducing capabilities that allow for the use of the XDGMesh object for treatment as geometry -- including both surface and volumetric meshes. The case of handling a surface mesh in particular isn't something the current UnstructuredMesh class is intended to handle and in that light there's some merit to having a separate mesh class for handling XDG meshes.

The idea once the aditional XDGUniverse class is introduced is for an XDGMesh object to be passed to the universe, allowing spatial data structures to be built only once for that particular mesh and leveraged for both geometry and tally operations as needed. For example, one might want to use an XDG volumetric mesh as a surface geometry in OpenMC but perform mesh tallies using only a volume or two.

It's still possible it will become clear that we really should wrap this class into the UnstructuredMesh class, but I think given the distinct usecases in upcoming PRs we can leave it as-is. Also, it will be easier to fold in than take out in the future if we need to so I think it's in a good place for now.

Additional thoughts are welcome as usual!

Comment thread openmc/mesh.py
@library.setter
def library(self, lib: str):
cv.check_value('Unstructured mesh library', lib, ('moab', 'libmesh'))
cv.check_value('Unstructured mesh library', lib, ('xdg','moab', 'libmesh'))
Copy link
Copy Markdown
Contributor

@jon-proximafusion jon-proximafusion May 20, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If xdgmesh is going to be separate to Unstructuredmesh then I guess 'xdg' should not be accepted here as a valid input but I guess this is needed when reading it in from the statepoint.

Suggested change
cv.check_value('Unstructured mesh library', lib, ('xdg','moab', 'libmesh'))
cv.check_value('Unstructured mesh library', lib, ('moab', 'libmesh'))

@jon-proximafusion
Copy link
Copy Markdown
Contributor

Thanks both, you two know the future of XDG so are best to decide if XDG should be its own class.

Just when I was viewing this as a self contained PR it looks like an UnstructuredMesh in three ways:

  • The Python validator on UnstructuredMesh.library accepts "xdg" (added in this PR at openmc/mesh.py:2657), implying XDGMesh is just another backend.

  • On statepoint readback, the mesh comes back as an openmc.mesh.UnstructuredMesh with library="xdg", not as an XDGMesh. So downstream code using statepoint.meshes[...], .volumes, .write_data_to_vtk(...) keeps working unchanged.

  • XDGMesh.from_hdf5 is a NotImplementedError stub, so the readback path relies entirely on UnstructuredMesh.from_hdf5.

As XDGMesh is meant to stay separate from UnstructuredMesh I would be keen to know what you have planned for reading the mesh results from the statepoint, will it be read back out as an xdgmesh in the future?

Copy link
Copy Markdown
Contributor

@paulromano paulromano left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very excited for this new feature @pshriwise! Here's a first round of comments:

Comment thread .github/workflows/ci.yml
Comment on lines 49 to 51
dagmc: [n]
xdg: [n]
libmesh: [n]
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wondering if maybe we should just have a single configuration that covers DAGMC, libMesh, and XDG. Do you see any downside of that?

Comment thread include/openmc/xdg.h
Comment on lines +89 to +99
//! Get the distance to the nearest boundary for a given position and
//! direction \param[in] g GeometryState object containing position and
//! direction \return NextMeshCell struct containing distance, face index, and
//! next indices
NextMeshCell distance_to_bin_boundary(GeometryState& g) const;

//! Get the distance to the nearest boundary for a given position and
//! direction \param[in] r Position to check \param[in] u Direction to check
//! \return Distance to the nearest boundary
NextMeshCell distance_to_bin_boundary(
int bin, const Position& r, const Direction& u) const;
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not defined/used

Comment thread include/openmc/mesh.h
Comment on lines +50 to +54
struct NextMeshCell {
double distance {INFTY};
int face_idx {-1};
std::array<int, 3> next_ijk;
};
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is only used by distance_to_bin_boundary, which itself is never used


if(@OPENMC_USE_OPENMP@)
find_package(OpenMP REQUIRED)
if(NOT TARGET OpenMP::OpenMP_CXX)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now that #3653 is merged, you'll need to update this file. One difference between #3653 and what you have here is that you also have the calls to if(NOT TARGET ...). Not sure if those are strictly necessary? @ahnaf-tahmid-chowdhury do you have an opinion on the matter?

Comment thread openmc/xdg.py
Comment on lines +102 to +105
@classmethod
def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
raise NotImplementedError("XDGMesh.from_hdf5 is not implemented")
# TODO: add length_multiplier
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you planning on adding this in this PR?

Comment thread src/xdg.cpp
// TODO: Make more robust (including mesh entrance/re-entrance)
xdg::Position p0 {r0.x, r0.y, r0.z};
xdg::Position p1 {r1.x, r1.y, r1.z};
double length_rcp = 1 / (p1 - p0).length();
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider length_inv, inverse_length, one_over_length, etc. I figured out "rcp" = reciprocal but it took me a few seconds of head scratching

Comment thread src/xdg.cpp
Comment on lines +119 to +127
auto track_segments = xdg_->segments(p0, p1);
// remove elements with lengths of zero
track_segments.erase(
std::remove_if(track_segments.begin(), track_segments.end(),
[](const std::pair<xdg::MeshID, double>& p) { return p.second == 0.0; }),
track_segments.end());
for (const auto& track_segment : track_segments) {
bins.push_back(mesh_id_to_bin(track_segment.first));
lengths.push_back(track_segment.second * length_rcp);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like segments returns a vector, which means this would end up allocating memory. Have you considered any design alternatives that would avoid this?

Comment on lines +317 to +330
# {
# "mesh_filename": "test_mesh_hexes.e",
# "mesh_kind": "hex",
# "holes": None,
# "elem_per_voxel": 1,
# "libraries": ("libmesh",),
# },
# {
# "mesh_filename": "test_mesh_hexes.exo",
# "mesh_kind": "hex",
# "holes": None,
# "elem_per_voxel": 1,
# "libraries": ("libmesh",),
# },
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are these commented out?

boundary_type='vacuum')
water_max_z = openmc.ZPlane(z0=1.0,
name="maximum z",
boundary_type='vacuum')
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All of these box-like regions could be greatly simplified with openmc.model.RectangularParallelepiped

model.settings.batches = 10

# source setup
r = openmc.stats.Uniform(a=0.0, b=0.0)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would effectively make it a point source, which doesn't seem to be what was intended.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants