Hi everyone,
I rotate a region around an axis, which is also the axis of a cylinder belonging to the region. Oddly enough, the cylinder (id=4) doesn’t merge with the new surfaces created and a new ‘quadric’ is created (id=8) that describes the same surface as the cylinder.
I get a strange behavior of the code : sometimes, a particle will get “trapped” on surface 4/8 and have as many events as possible (1 million in my case) before the code switch to the next particle.
An easy fix is to replace all 8 by 4 in the xml. I can also patch that behavior in the case of my project. Still, it might be worth further investigations. My way to patch this would be to fix the remove_redundant_surface by converting all surfaces coefficients in quadric coefficients to list redundant surfaces.
import openmc
import numpy as np
def rot(region, angle_in_deg):
return region.rotate((0.,0.,angle_in_deg), inplace=False)
mat0 = openmc.Material()
mat0.set_density("g/cm3", 10.)
mat0.add_element("O", 2.)
mat0.add_nuclide("U235", 0.93)
mat0.add_nuclide("U238", 0.07)
mat1 = openmc.Material()
mat1.set_density("g/cm3", 1.)
mat1.add_element("O", 2.)
mat1.add_element("H", 1)
materials = openmc.Materials([mat0, mat1])
materials.export_to_xml()
box = +openmc.XPlane(-500.) & -openmc.XPlane(-100.) & -openmc.XCylinder(r=60.)
drill = openmc.ZCylinder(r=165.)
drilled0 = box & +drill
drilled1 = rot(drilled0, 120.)
drilled2 = rot(drilled0, 240.)
cell1 = openmc.Cell(fill=mat0, region = drilled0)
cell2 = openmc.Cell(fill=mat0, region = drilled1)
cell3 = openmc.Cell(fill=mat0, region = drilled2)
cell4 = openmc.Cell(fill=mat1, region = ~drilled0 & ~drilled1 & ~drilled2)
universe = openmc.Universe(cells=[cell1, cell2, cell3, cell4])
xP0 = openmc.XPlane(-1000., boundary_type="reflective")
xP1 = openmc.XPlane(+1000., boundary_type="reflective")
yP0 = openmc.YPlane(-1000., boundary_type="reflective")
yP1 = openmc.YPlane(+1000., boundary_type="reflective")
zP0 = openmc.ZPlane(-1000., boundary_type="reflective")
zP1 = openmc.ZPlane(+1000., boundary_type="reflective")
world = +xP0 & -xP1 & +yP0 & -yP1 & +zP0 & -zP1
root = openmc.Cell(fill=universe, region = world)
geom = openmc.Geometry(root=[root])
geom.merge_surfaces = True
geom.export_to_xml()
point = openmc.stats.Point((-200.,0.,0.))
source = openmc.IndependentSource(space=point)
settings = openmc.Settings()
settings.particles = 1000
settings.inactive = 10
settings.batches = 1000
settings.source = source
settings.export_to_xml()
<geometry>
<cell id="1" material="1" region="1 -2 -3 4" universe="1"/>
<cell id="2" material="1" region="5 -6 -7 8" universe="1"/>
<cell id="3" material="1" region="9 -10 -11 8" universe="1"/>
<cell id="4" material="2" region="~(1 -2 -3 4) ~(5 -6 -7 8) ~(9 -10 -11 8)" universe="1"/>
<cell fill="1" id="5" region="13 -14 15 -16 17 -18" universe="2"/>
<surface coeffs="-500.0" id="1" type="x-plane"/>
<surface coeffs="-100.0" id="2" type="x-plane"/>
<surface coeffs="0.0 0.0 60.0" id="3" type="x-cylinder"/>
<surface coeffs="0.0 0.0 165.0" id="4" type="z-cylinder"/>
<surface coeffs="-0.49999999999999983 0.8660254037844387 0.0 -500.0" id="5" type="plane"/>
<surface coeffs="-0.49999999999999983 0.8660254037844387 0.0 -100.0" id="6" type="plane"/>
<surface coeffs="0.7500000000000003 0.2499999999999999 1.0000000000000002 0.8660254037844386 -0.0 0.0 -0.0 -0.0 0.0 -3600.000000000001" id="7" type="quadric"/>
<surface coeffs="1.0 1.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -27225.0" id="8" type="quadric"/>
<surface coeffs="-0.5000000000000004 -0.8660254037844384 0.0 -500.0" id="9" type="plane"/>
<surface coeffs="-0.5000000000000004 -0.8660254037844384 0.0 -100.0" id="10" type="plane"/>
<surface coeffs="0.7499999999999996 0.25000000000000044 1.0 -0.8660254037844392 0.0 0.0 0.0 -0.0 -0.0 -3600.0" id="11" type="quadric"/>
<surface boundary="reflective" coeffs="-1000.0" id="13" type="x-plane"/>
<surface boundary="reflective" coeffs="1000.0" id="14" type="x-plane"/>
<surface boundary="reflective" coeffs="-1000.0" id="15" type="y-plane"/>
<surface boundary="reflective" coeffs="1000.0" id="16" type="y-plane"/>
<surface boundary="reflective" coeffs="-1000.0" id="17" type="z-plane"/>
<surface boundary="reflective" coeffs="1000.0" id="18" type="z-plane"/>
</geometry>