[BUG] - lost particles on openmc xml geometries

[shimwell]

Describe the bug

I've been running the geometry xml files produced with openmc and have noticed a few files fail in transport with lost particles. I converted all 52 of the files in the testings/inputStep folder and only found 8 files with leaky geometry, which I think is quite good as the other 44 work nicely.

I've made some CI to test geometry produced work in transport without leakage and while testing locally noticed not all the geometry files pass the new test

To Reproduce

I've made a script that will convert all the stp files and then run them in openmc. for this to work you must have openmc and freecad installed. conda install -c conda-forge openmc freecad doe the trick for me.

try:
    import freecad  # importing conda package if present
except:
    pass
import Part
from FreeCAD import Import
import openmc
from pathlib import Path
import GEOUNED

input_step_files = [
    Path('testing/inputSTEP/large/SCDR.stp'),
    Path('testing/inputSTEP/placa.stp'),
    Path('testing/inputSTEP/placa2.stp'),
    Path('testing/inputSTEP/Misc/sphereBarCyl2.stp'),
    Path('testing/inputSTEP/Misc/sphereBarCyl1.stp'),
    Path('testing/inputSTEP/DoubleCylinder/placa3.step'),
    Path('testing/inputSTEP/DoubleCylinder/placa.stp'),
    Path('testing/inputSTEP/Torus/face2.stp'),
]

for leaky_stp_file in input_step_files:

    template = (
        "[Files]\n"
        "title = 'Input Test'\n"
        f"stepFile = {leaky_stp_file.resolve()}\n"
        f"geometryName = {leaky_stp_file.stem}\n"
        "outFormat = ('openMC_XML')\n"
        "[Parameters]\n"
        "compSolids = False\n"
        "volCARD = False\n"
        "volSDEF = True\n"
        "voidGen = True\n"
        "dummyMat = True\n"
        "minVoidSize = 100\n"
        "cellSummaryFile = False\n"
        "cellCommentFile = False\n"
        "debug = False\n"
        "simplify = full\n"
        "[Options]\n"
        "forceCylinder = False\n"
        "splitTolerance = 0\n"
        "newSplitPlane = True\n"
        "nPlaneReverse = 0\n"
    )

    with open("config.ini", mode="w") as file:
        file.write(template)

    GEO = GEOUNED.GEOUNED('config.ini')
    GEO.SetOptions()
    GEO.outFormat = ("openMC_XML")
    GEO.Start()

    Import.insert(str(leaky_stp_file), "converted-cad")
    result = Part.Shape()
    result.read(str(leaky_stp_file))

    openmc_xml_file=Path(leaky_stp_file.stem).with_suffix('.xml')
    materials = openmc.Materials()
    geometry = openmc.Geometry.from_xml(path=openmc_xml_file, materials=materials)
    cells = geometry.get_all_cells()

    bb = result.BoundBox
    llx, lly, llz, urx, ury, urz = bb.XMin, bb.YMin, bb.ZMin, bb.XMax, bb.YMax, bb.ZMax
    llx, lly, llz, urx, ury, urz = llx/10, lly/10, llz/10, urx/10, ury/10, urz/10

    source = openmc.IndependentSource()
    source.space = openmc.stats.Box(
        lower_left=(llx, lly, llz),
        upper_right=(urx, ury, urz)
    )
    source.energy = openmc.stats.Discrete([14e6], [1])

    materials = openmc.Materials()
    geometry = openmc.Geometry.from_xml(openmc_xml_file, materials)

    settings = openmc.Settings()
    settings.batches = 10
    settings.max_lost_particles = 1
    # could be increased to be even more sure there are no leaky particles
    settings.particles = 100_000
    settings.run_mode = "fixed source"
    settings.source = source
    model = openmc.Model(geometry, materials, settings)

    model.run()

IMPORTANT I've identified the bug surfaces for each model in the code snippet

Expected behavior

These geometries should all run without losing any particles

Error message

After particle 23372 crossed surface 5 it could not be located in any
          cell and it did not leak.

Screenshots

If applicable, add screenshots to help explain your problem.

Please complete the following information):

• OS: ubuntu
• GEOUNED dev
• Python version 3.11

[shimwell]

setting simplify = 'no' gets a few of the geometries working. :tada:

However a few still don't work (placa, placa2 placa3) I have updated the example code to use simplify =no and removed the step files that work in this mode

try:
    import freecad  # importing conda package if present
except:
    pass
import Part
from FreeCAD import Import
import openmc
from pathlib import Path
import GEOUNED

# uncomment the stp file you want to check and run the script

input_step_file = Path('testing/inputSTEP/large/SCDR.stp')
# input_step_file = Path('testing/inputSTEP/placa.stp') # only has a single plane as a vacuum
# input_step_file = Path('testing/inputSTEP/placa2.stp') # crossing surface 33
# input_step_file = Path('testing/inputSTEP/DoubleCylinder/placa3.step') # crossing surface 51

# fails with a different error, broken geometry.xml?
# input_step_file = Path('testing/inputSTEP/DoubleCylinder/placa.stp') 
#  geometry = openmc.Geometry.from_xml(path=openmc_xml_file, materials=materials)

template = (
    "[Files]\n"
    "title = 'Input Test'\n"
    f"stepFile = {input_step_file.resolve()}\n"
    f"geometryName = {input_step_file.stem}\n"
    "outFormat = ('openMC_XML')\n"
    "[Parameters]\n"
    "compSolids = False\n"
    "volCARD = False\n"
    "volSDEF = True\n"
    "voidGen = True\n"
    "dummyMat = True\n"
    "minVoidSize = 100\n"
    "cellSummaryFile = False\n"
    "cellCommentFile = False\n"
    "debug = False\n"
    "simplify = no\n"
    "[Options]\n"
    "forceCylinder = False\n"
    "splitTolerance = 0\n"
    "newSplitPlane = True\n"
    "nPlaneReverse = 0\n"
)

with open("config.ini", mode="w") as file:
    file.write(template)

GEO = GEOUNED.GEOUNED('config.ini')
GEO.SetOptions()
GEO.outFormat = ("openMC_XML")
GEO.Start()

Import.insert(str(input_step_file), "converted-cad")
result = Part.Shape()
result.read(str(input_step_file))

openmc_xml_file=Path(input_step_file.stem).with_suffix('.xml')
materials = openmc.Materials()
geometry = openmc.Geometry.from_xml(path=openmc_xml_file, materials=materials)
cells = geometry.get_all_cells()

bb = result.BoundBox
llx, lly, llz, urx, ury, urz = bb.XMin, bb.YMin, bb.ZMin, bb.XMax, bb.YMax, bb.ZMax
llx, lly, llz, urx, ury, urz = llx/10, lly/10, llz/10, urx/10, ury/10, urz/10

source = openmc.IndependentSource()
source.space = openmc.stats.Box(
    lower_left=(llx, lly, llz),
    upper_right=(urx, ury, urz)
)
source.energy = openmc.stats.Discrete([14e6], [1])

materials = openmc.Materials()
geometry = openmc.Geometry.from_xml(openmc_xml_file, materials)

settings = openmc.Settings()
settings.batches = 10
settings.max_lost_particles = 1
# could be increased to be even more sure there are no leaky particles
settings.particles = 100_000_00
settings.run_mode = "fixed source"
settings.source = source
model = openmc.Model(geometry, materials, settings)

model.run()

[psauvan]

I have found two geometries producing lost particles (placa3.stp and Triangle.stp). The first geometry (placa3.stp) lost particles because two planes that should be identical have been considered as differents. This has been introduced by the FreeCAD spliting method. As you can in the pictures below the plate is cut by several planes. The plane identified in the first picture is reproduced as two planes in the second picture. To avoid lost particle in this case plane tolerance has to be set to 0.05 planeDistance = 0.05

For the second model, the lost particle is due to numerical precision leading to overlaping cell in a small region as show in the picture. To resolve this lost particle problem the overlapping should be fix manually or using the GEOUNED option forceNoOverlap = True.

If the lost particles you found was in that two models, the issue can be considered as closed.

[shimwell]

Super I can try these changes locally and then I can make a PR to change the CI slightly and we can see them working on the CI.

Thanks for looking into this.

I guess let us keep the issue open a bit while I follow up on these useful suggestions.