DAGMC CAD-based geometry as universes

[pshriwise]

We currently only support DAGMC models where the entire geometry is represented as CAD-based tesselations. This PR adds a DAGMC Universe derived from the standard CSG universe. This allows one to construct geometries from both CAD and CSG. This enables repeated CAD geometries in OpenMC like this lattice of DAGMC pincells from our regression tests

Or, as a more interesting case, one can rotate a section of a tokamak into a set of CSG cells to model the entire device with a smaller DAGMC memory footprint. (Thanks @shimwell for providing this quarter-tokamak model as a demonstration.)

Here is the entire tokamak model in Cubit:

We can add this model as a fill to the appropriate cells and rotate the model 90 degrees to get this X-Y view:

Adding two more cells and rotations gives use the entire device:

Input Changes

These updates do require some changes to the current format when the DAGMC model represents the entire geometry. A geometry.xml file is now always required. To represent a geometry composed of a single DAGMC model as we do now, a geometry.xml file that gives the same behavior should contain

<?xml version='1.0' encoding='utf-8'?>
<geometry>
<dagmc id="1" filename="dagmc.h5m"/>
</geometry>

where the id is the universe ID and the filename is the name of the DAGMC model file.

Note that the filename parameter remoes the restriction that DAGMC files have the name dagmc.h5m.

Notable code changes

• Addition of the DAGUniverse class
• Addition of Cell::geom_type and Surface::geom_type
• Removal of the settings::dagmc parameter from the CAPI and Python API
• Added a Universe::find_cell method
• Fewer #ifdef DAGMC's throughout the code thanks to the addition of the GeometryType enum

ID management

One difficulty of mixing CSG and CAD geometry is the unknown ID spaces of the CAD geometry before simulation initialization. The DAGUniverse class on the Python side has an option, auto_geom_ids, to automatically set cell and surface IDs outside of the CSG cell & surface space to avoid conflicts. If this option is not enabled, then an error will occur if the ID spaces overlap.

A similar option, auto_mat_ids, is also available for materials to address the case where some of the DAGMC universes are using UWUW materials (material definitions embedded in the .h5m along with the geometry.

Particle tracking considerations

Coincident DAGMC/CSG surfaces

Thanks to our preference for hitting surface at higher levels in the geometry, this is hasn't been a problem in the models I've used for testing.

Managing a particle's ray history

The ray history is a record of facets hit by DAGMC rays to avoid duplicate hits when coincident with a surface or reflecting from a surface. This history is reset when

• a particle changes direction or undergoes a collision,
• a particle is not on a surface,
• and a particle is terminated.

Now that a particle can move directly from one DAGMC universe to another, another case is added for

• a particle exiting the current DAGUniverse

For more details on this implementation, please see this google doc w/ the proposal for this work.

[makeclean]

Ooh very exciting :)

…

________________________________ From: Patrick Shriwise ***@***.***> Sent: Tuesday, April 27, 2021 10:12:21 PM To: openmc-dev/openmc ***@***.***> Cc: Davis, Andrew ***@***.***>; Review requested ***@***.***> Subject: Re: [openmc-dev/openmc] DAGMC CAD-based geometry as universes (#1825) @pshriwise<https://github.com/pshriwise> requested your review on: #1825<#1825> DAGMC CAD-based geometry as universes. — You are receiving this because your review was requested. Reply to this email directly, view it on GitHub<#1825 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AASTUSV3N6KP3ENVKEFXRCLTK4SDLANCNFSM43VYSAZA>.

[gridley]

In case anyone else is wondering, this seems to add pretty much zero overhead to non-DAGMC models. Seems like a really useful upgrade.

[pshriwise]

Sorry for the initial CI failures here. Needed to make a correction to the properties of the DAGUniverse Python class. This should be good to try out now!

[shimwell]

@pshriwise is it worth me adding a rotation to this example #1800 . I assume that I just remove the reflecting surface (perhaps using the new dagmc_remove_tags package) and rotate 90, 180 and 270 degrees to make the complete geometry

[pshriwise]

@pshriwise is it worth me adding a rotation to this example #1800 . I assume that I just remove the reflecting surface (perhaps using the new dagmc_remove_tags package) and rotate 90, 180 and 270 degrees to make the complete geometry

I'd welcome it! It would make #1800 dependent on this PR though. I'll leave it up to you as to whether or not you're ok with waiting on these changes to get that example in.

[makeclean]

Looks mostly good here, I think a few extra new lines, some extra comments and its good :)

[makeclean]

will u->id_ always be sorted smallest to largest?

[pshriwise]

It isn't a requirement for model::universes, no. Are you concerned the next ID isn't correct?

[makeclean]

just wondering if it leads to any issues, like does it matter if the order isnt monotonically increasing?

[pshriwise]

The IDs can be in any order here. Whenever we find one that's bigger, we update next_univ_id. If we happened to find the largest ID first, we'd only update that variable once.

[pshriwise]

Are you happy with this implementation @makeclean?

[makeclean]

Do you need a graveyard even in a universe based model? Do OpenMC cells not bound this?

[pshriwise]

The only case where you now need one is if the DAGMC model is the root universe of the Geometry. I'll add some detection for that so this warning doesn't appear when it shouldn't. Thanks!

[pshriwise]

Reorganized the code a bit so we can perform this check in 648054f. Hopefully this is what you had in mind!

[shimwell]

Looks like code block 17 from the cad-based-geometry.ipynb has the old settings.dagmc = True method

https://github.com/pshriwise/openmc/blob/5b20fbbf6a377b16d201a34c5cf65471662ac7cc/examples/jupyter/cad-based-geometry.ipynb#L524

[shimwell]

I've been using the branch on an ARC reactor model and have reproduced the example pictures with a different geometry (still 90 degree rotated 3 times).

I've found the ability to turn off the auto assign id numbers useful
It is nice to have a geometry.xml file when using DAGMC
It is great that this branch allows flexibility with the dagmc h5m file (previously only dagmc.h5m was accepted).

Actually this feature really saved me and was exactly what I needed to get some simulation results.

Perhaps this is unrelated to this PR but one thing I noticed was that the terminal prints out a lot of (perhaps too much) information when plotting an DAGMC geometry using this method

import os
# makes the 3d "cube" style geometry
vox_plot = openmc.Plot()
vox_plot.type = 'voxel'
vox_plot.width = (1700., 1700., 1800.)
vox_plot.pixels = (300, 300, 300)
vox_plot.filename = 'plot_3d_tokamak'
vox_plot.color_by = 'material'
plots = openmc.Plots([vox_plot])
plots.export_to_xml()
openmc.plot_geometry()
os.system('openmc-voxel-to-vtk plot_3d_tokamak.h5 -o plot_3d_tokamak.vti')

[pshriwise]

import os
# makes the 3d "cube" style geometry
vox_plot = openmc.Plot()
vox_plot.type = 'voxel'
vox_plot.width = (1700., 1700., 1800.)
vox_plot.pixels = (300, 300, 300)
vox_plot.filename = 'plot_3d_tokamak'
vox_plot.color_by = 'material'
plots = openmc.Plots([vox_plot])
plots.export_to_xml()
openmc.plot_geometry()
os.system('openmc-voxel-to-vtk plot_3d_tokamak.h5 -o plot_3d_tokamak.vti')

Yeah, I'm seeing that too. Thanks for sharing this. Appears to be a warning about a corner case in MOAB's triangle intersection routine. I'll have to look into that separately.

[pshriwise]

Looks like code block 17 from the cad-based-geometry.ipynb has the old settings.dagmc = True method

Thanks! Just took care of that one.

[pshriwise]

I believe I've addressed all of your comments @paulromano, @makeclean. Let me know if I've missed anything!

I did some further testing using the CSG and DAGMC ATR models originally used for comparison in our DAGMC implementation.

Model	k-eff	std. dev. (pcm)
CSG	1.02115	9
DAGMC	1.02123	10
DAGMC/ w CSG bounds	1.02123	10

I ran the CSG and DAGMC models w/ 200000 particles per batch, 500 batches, 20 inactive. I then replaced the external boundary volume of the DAGMC model (a.k.a. the graveyard volume in DAGMC terms) with CSG planes and a vacuum BC.
A 200 x 200 x 1 mesh tally ranging from +/- 100 cm in X and Y and +/- 312 in Z (full axial height) was applied to produce the comparisons in flux shown below. (The flux units are unadjusted for reactor power as that isn't necessary for the purposes of comparison.)

I'd expect that the flux comparison between the DAGMC model and the DAGMC model bounded by CSG planes would be exactly the same, which is what we see here. I don't expect the comparison between the CSG and DAGMC models to be the same due to small differences in the geometry caused by the CAD tessellation, however. The results looked reasonable to me, with the largest relative differences in the flux occurring near the edge of the core where the variance in either flux tally is highest. Overall, very few of the flux voxels were statistically different.

I then inserted CSG X and Y planes at the origin w/ reflective BC's to test individual quadrants of the model.

Model	k-eff	std. dev.(pcm)
CSG Q1	1.00750	9
DAGMC Q1	1.00757	10
----------	---------	----------------
CSG Q2	1.01705	10
DAGMC Q2	1.01695	10
----------	---------	----------------
CSG Q3	1.03430	10
DAGMC Q3	1.03425	9
----------	---------	----------------
CSG Q4	1.02592	9
DAGMC Q4	1.02603	9

Again, the results looked reasonable when comparing the the CSG model for each quadrant to the DAGMC model. I've included the results of this comparison for quadrant 4:

There were a few statistically different results in the mesh, but, as I mentioned above, this isn't unexpected due to the fact that the geometries are in fact different.

[paulromano]

@pshriwise Thanks for incorporating all the feedback. I have a few more requests:

1. The new <dagmc_universe> element needs to be documented in docs/source/io_formats/geometry.rst
2. The geom_type attribute for surfaces needs to be documented in docs/source/io_formats/summary.rst
3. It looks like a geometry that has a DAGMC universe does not "round trip" correctly to and from XML. That is, if you geom.export_to_xml() and then new_geom = openmc.Geometry.from_xml(...); new_geom.export_to_xml('new_geom.xml'), the new_geom.xml file is not the same as the original.

[pshriwise]

@pshriwise Thanks for incorporating all the feedback. I have a few more requests:

1. The new <dagmc_universe> element needs to be documented in docs/source/io_formats/geometry.rst
2. The geom_type attribute for surfaces needs to be documented in docs/source/io_formats/summary.rst
3. It looks like a geometry that has a DAGMC universe does not "round trip" correctly to and from XML. That is, if you geom.export_to_xml() and then new_geom = openmc.Geometry.from_xml(...); new_geom.export_to_xml('new_geom.xml'), the new_geom.xml file is not the same as the original.

Thanks for @paulromano! Docs have all been updated.

I've updated the export/import code for XML so that it handles DAGMC universes properly and can reproduce the same geometry.xml reliably. I also decided that we should really be able to reproduce the geometry.xml from a summary.h5 file as well, so I've added a DAGUniverse::to_hdf5 method that writes out all of the necessary information to reproduce the Python object (and corresponding XML element if desired). I've updated the summary.h5 documentation to reflect these changes as well.

I noticed along the way that surface names weren't being picked up in openmc.Surface.from_xml_element, so I took care of that in 25cf7a3.

[paulromano]

Thanks for the latest round of updates. A few more requests based on those changes... we're getting close!

[paulromano]

@makeclean Do you have any further comments on this PR?

[makeclean]

LGTM thanks @pshriwise