Issue with DAGMC export in OpenMC (lost particles)

[generein]

I am using the following code to generate a paramak shell geometry, though the issue also happens with cylinders (paramak.CenterColumnShieldCylinder).

R = 5
t2 = 2
t3 = 3
t4 = 4
sphere_inner = paramak.SphericalShell(inner_radius=0,shell_thickness=R,name='sphere-inner',rotation_angle=360)
sphere_second = paramak.SphericalShell(inner_radius=R,shell_thickness=t2,name='sphere-second',rotation_angle=360)
sphere_third = paramak.SphericalShell(inner_radius=R+t2,shell_thickness=t3,name='sphere-third',rotation_angle=360)
sphere_fourth = paramak.SphericalShell(inner_radius=R+t2+t3,shell_thickness=t4,name='sphere-fourth',rotation_angle=360)
my_reactor = paramak.Reactor(shapes_and_components = [sphere_inner,sphere_second,sphere_third,sphere_fourth])
reactor_descr = 'spherical_shell_example'
my_reactor.export_stp(reactor_descr+'.stp',units='m')
my_reactor.export_dagmc_h5m(reactor_descr+'.h5m',min_mesh_size = 0.01, max_mesh_size = 20.)

I then load this .h5m file into an OpenMC environment and follow this example

import openmc
import openmc_plasma_source as ops  
import openmc_tally_unit_converter as otuc
import neutronics_material_maker as nmm

filedir = '/path/to/file/'
filename = 'spherical_shell_example.h5m'

bound_dag_univ = openmc.DAGMCUniverse(filename=filedir+filename).bounded_universe()
my_geometry = openmc.Geometry(root=bound_dag_univ)

mat1 = nmm.Material.from_library(name='DT_plasma').openmc_material
mat1.name = 'sphere-inner'
mat2 = nmm.Material.from_library(name='tungsten').openmc_material
mat2.name = 'sphere-second'
mat3 = nmm.Material.from_library(name='lithium-lead',temperature=900.0).openmc_material
mat3.name = 'sphere-third'
mat4 = nmm.Material.from_library(name='SS_316L_N_IG').openmc_material
mat4.name = 'sphere-fourth'

materials = openmc.Materials([mat1, mat2, mat3, mat4])

tally1 = openmc.Tally()
material_filter = openmc.MaterialFilter(mat3) 
tally1 = openmc.Tally(name='blanket_heating')
tally1.filters = [material_filter]
tally1.scores = ['heating']

my_tallies = openmc.Tallies([tally1])

my_settings = openmc.Settings(batches = 1, particles = 100, run_mode = 'fixed source')

my_settings.source = ops.FusionPointSource(fuel="DT", temperature=20000.0, coordinate=(0,0,0))

my_model = openmc.Model(
    materials=materials, geometry=my_geometry, settings=my_settings, tallies=my_tallies
)

statepoint_file = my_model.run()

sp = openmc.StatePoint(statepoint_file)

heating_tally = sp.get_tally(name='blanket_heating')
heating_tally.get_values().sum()

tally = otuc.process_tally(
    heating_tally,
)
print(tally)

This results in an error of 100% leakage, i.e. all particles lost, even at increased batches and/or particle numbers. DAGMC checks on the geometry show that it has no leaky surfaces or volumes. Material assignments work/load fine.

[Adambar24]

Hi, I dont know which version of paramak you are using but I do not have an attribute called SphericalShell and I cant find it on paramak API documents (https://paramak.readthedocs.io/en/main/API-Reference.html#parametric-components). Though this kind of problem is simple and robust using openmc's geometry classes. I have created an example below using your variables:

R = 5
t2 = 2
t3 = 3
t4 = 4

Sphere_inner = openmc.Sphere(r= R)
Sphere_second = openmc.Sphere(r= R + t2)
Sphere_third = openmc.Sphere(r = R + t2 + t3)
Sphere_fourth = openmc.Sphere(r = R + t2 + t3 + t4, boundary_type = 'vacuum') 
#Note that this surface is a vacuum since this is no longer a bound universe 

######## Create Regions and Cells ########

Inner_Region = -Sphere_inner
Second_Region = +Sphere_inner & -Sphere_second
Third_Region = +Sphere_second & -Sphere_third
Fourth_Region = +Sphere_third & -Sphere_fourth

Inner_cell = openmc.Cell(name='Inner_Sphere', fill= mat1, region= Inner_Region)
Second_cell = openmc.Cell(name='Second_Shell', fill= mat2, region= Second_Region)
Third_cell = openmc.Cell(name='Third_Shell', fill= mat3, region= Third_Region)
Fourth_cell = openmc.Cell(name='Fourth_Shell', fill= mat4, region= Fourth_Region)

######## Create Universe #######

root_universe = openmc.Universe(0, name='Root Universe')
root_universe.add_cells([Inner_cell, Second_cell, Third_cell, Fourth_cell])

geometry = openmc.Geometry(root_universe)

This should create a shell like in the image below. This is a 2d plane (xy) at the origin

[generein]

Thanks, @Adambar24 . I am now indeed using the OpenMC CSG geometries directly, but longer term I want to be able to use the STP export too. The SphericalShell is part of this PR.

[shimwell]

Thanks for reporting this, I have investigated a bit and found that setting the rotation angle to 180 appears to help.

I guess the surface senses are getting mixed up when a full rotation happens and the two cad surfaces that come into contact with each other are not getting cleaned up correctly.

So I'm thinking it can work if we make a sector model, would that be sufficient for now.

This is a bit trickier than it should be but here is an example of how to make a reflecting surface sector model, this one is 90 degrees so needs two surfaces, you would just need one

Alternatively we can make a sector model and get openmc to repeat/roated the geometry so it looks and behaves like a full 360 model. Not got an example but I can make one if needed

[shimwell]

I have done the sector model option in the code below, to note, the source is nudged inside the sector, the bounding region is now manually made, let me know if you want a 360 model by rotation the 180 degree version

import paramak
R = 5
t2 = 2
t3 = 3
t4 = 4
sphere_inner = paramak.SphericalShell(inner_radius=0,shell_thickness=R,name='sphere-inner',rotation_angle=180)
sphere_second = paramak.SphericalShell(inner_radius=R,shell_thickness=t2,name='sphere-second',rotation_angle=180)
sphere_third = paramak.SphericalShell(inner_radius=R+t2,shell_thickness=t3,name='sphere-third',rotation_angle=180)
sphere_fourth = paramak.SphericalShell(inner_radius=R+t2+t3,shell_thickness=t4,name='sphere-fourth',rotation_angle=180)
my_reactor = paramak.Reactor(shapes_and_components = [sphere_inner,sphere_second,sphere_third,sphere_fourth])
reactor_descr = 'spherical_shell_example'
my_reactor.export_stp(reactor_descr+'.stp',units='m')
my_reactor.export_dagmc_h5m(reactor_descr+'.h5m',min_mesh_size = 0.01, max_mesh_size = 20.)

import openmc
import openmc_plasma_source as ops  
import openmc_tally_unit_converter as otuc
import neutronics_material_maker as nmm
import math

filedir = ''
filename = 'spherical_shell_example.h5m'

# makes use of the dagmc geometry
dag_univ = openmc.DAGMCUniverse(filedir+filename)

# creates an edge of universe boundary surface
vac_surf = openmc.Sphere(r=10000, surface_id=9999, boundary_type="vacuum")

# adds reflective surface for the sector model at 0 degrees
# you could do this differently but I've done it like this so if you want to
# change the angle it is easy
reflective_1 = openmc.Plane(
    a=math.sin(0),
    b=-math.cos(0),
    c=0.0,
    d=0.0,
    surface_id=9991,
    boundary_type="reflective",
)

# specifies the region as below the universe boundary and inside the reflective surfaces
region = -vac_surf & -reflective_1

# creates a cell from the region and fills the cell with the dagmc geometry
containing_cell = openmc.Cell(cell_id=9999, region=region, fill=dag_univ)

my_geometry = openmc.Geometry(root=[containing_cell])

mat1 = nmm.Material.from_library(name='DT_plasma').openmc_material
mat1.name = 'sphere-inner'
mat2 = nmm.Material.from_library(name='tungsten').openmc_material
mat2.name = 'sphere-second'
mat3 = nmm.Material.from_library(name='lithium-lead',temperature=900.0).openmc_material
mat3.name = 'sphere-third'
mat4 = nmm.Material.from_library(name='SS_316L_N_IG').openmc_material
mat4.name = 'sphere-fourth'

my_materials = openmc.Materials([mat1, mat2, mat3, mat4]) # gets tally for all 4 materials in one go

tally1 = openmc.Tally()
material_filter = openmc.MaterialFilter([mat1, mat2, mat3, mat4]) 
tally1 = openmc.Tally(name='blanket_heating')
tally1.filters = [material_filter]
tally1.scores = ['heating']

my_tallies = openmc.Tallies([tally1])

my_settings = openmc.Settings(batches = 10, particles = 1000, run_mode = 'fixed source')

# note the source has been moved slightly to be inside the geometry
my_settings.source = ops.FusionPointSource(fuel="DT", temperature=20000.0, coordinate=(0,0.0001,0))

my_model = openmc.Model(
    materials=my_materials, geometry=my_geometry, settings=my_settings, tallies=my_tallies
)

statepoint_file = my_model.run()

sp = openmc.StatePoint(statepoint_file)

heating_tally = sp.get_tally(name='blanket_heating')
print(heating_tally.mean)  # units a eV per source particle