KUCA Modeling problems

Hello everybody, recently I have carried out the KUCA modeling(Case Ⅱ-1).


Here are my Python script


import openmc
import os
os.chdir('kuka')
# ==================================================
fuel=openmc.Material(name='93 fuel')
fuel.add_nuclide('U235',1.50694E-03)
fuel.add_nuclide('U238',1.08560E-04)
fuel.add_nuclide('Al27',5.56436E-02)
fuel.set_density('sum')#
# ==================================================
air=openmc.Material(name='air')
air.add_element('N' ,0.000039)
air.add_element('O' ,0.000011)
air.set_density('atom/b-cm', 5.000E-05)
# ==================================================
al=openmc.Material(name='Al')
al.add_nuclide('Al27',1)
al.set_density('atom/b-cm',6.00385E-02)
# ==================================================
# PE1/8
pe1_8=openmc.Material(name='PET1/8')
pe1_8.add_element('C',4.01084E-02)
pe1_8.add_element('H',8.02167E-02)
pe1_8.set_density('atom/cm3',(4.01084E-02+8.02167E-02)*1e24)
# ==================================================
# PE1/4
pe1_4=openmc.Material(name='PET1/4')
pe1_4.add_element('C',4.04356E-02)
pe1_4.add_element('H',8.08711E-02)
pe1_4.set_density('sum')
# ==================================================
# PE1/2
pe1_2=openmc.Material(name='PET1/2')
pe1_2.add_element('C',4.03280E-02)
pe1_2.add_element('H',8.06560E-02)
pe1_2.set_density('sum')
# ==================================================
# PE19
pe19=openmc.Material(name='PET19')
pe19.add_element('C',4.00042E-02)
pe19.add_element('H',8.00083E-02)
pe19.set_density('sum')
# ==================================================
# Pb-Bi
pb_bi=openmc.Material(name='pb_bi')
pb_bi.add_nuclide('Pb204',1.87461E-04)
pb_bi.add_nuclide('Pb206',3.25860E-03)
pb_bi.add_nuclide('Pb207',3.00266E-03)
pb_bi.add_nuclide('Pb208',7.15378E-03)
pb_bi.add_nuclide('Bi209',1.67670E-02)
pb_bi.set_density('sum')
# ==================================================
# Pb-Bi coating 1
pb_bi_c1=openmc.Material(name='pb_bi coating 1st')
pb_bi_c1.add_nuclide('H1',2.83301E-03)
pb_bi_c1.add_nuclide('H2',4.25015E-07)
pb_bi_c1.add_nuclide('O16',2.06885E-03+7.88039E-07+4.14757E-06)
pb_bi_c1.add_element('C',2.27058E-03)
pb_bi_c1.add_nuclide('Ti46',2.50941E-05)
pb_bi_c1.add_nuclide('Ti47',2.28983E-05)
pb_bi_c1.add_nuclide('Ti48',2.31493E-04)
pb_bi_c1.add_nuclide('Ti49',1.72522E-05)
pb_bi_c1.add_nuclide('Ti50',1.69385E-05)
pb_bi_c1.add_nuclide('S32',4.16818E-04+8.77326E-08)
pb_bi_c1.add_nuclide('S33',3.28997E-06)
pb_bi_c1.add_nuclide('S34',1.84677E-05)

pb_bi_c1.add_nuclide('Ba132',4.43050E-07)
pb_bi_c1.add_nuclide('Ba134',1.06025E-05)
pb_bi_c1.add_nuclide('Ba135',2.89167E-05)
pb_bi_c1.add_nuclide('Ba136',3.44526E-05)
pb_bi_c1.add_nuclide('Ba137',4.92619E-05)
pb_bi_c1.add_nuclide('Ba138',3.14521E-04)
pb_bi_c1.set_density('sum')
# ==================================================
# Pb-Bi coating 1
pb_bi_c2=openmc.Material(name='pb_bi coating 2nd')
pb_bi_c2.add_nuclide('H1',3.78991E-03)
pb_bi_c2.add_nuclide('H2',5.64072E-07)
pb_bi_c2.add_nuclide('O16',4.58782E-04+1.74753E-07+9.19753E-07)
pb_bi_c2.add_element('C',4.03671E-03)
pb_bi_c2.add_nuclide('Ti46',5.36896E-05)
pb_bi_c2.add_nuclide('Ti47',4.89918E-05)
pb_bi_c2.add_nuclide('Ti48',4.95287E-04)
pb_bi_c2.add_nuclide('Ti49',3.69116E-05)
pb_bi_c2.add_nuclide('Ti50',3.62405E-05)
pb_bi_c2.add_nuclide('Si28',1.86243E-05)
pb_bi_c2.add_nuclide('Si29',9.43026E-07)
pb_bi_c2.add_nuclide('Si30',6.25992E-07)
pb_bi_c2.set_density('sum')
# ==================================================
# control rod
rod=openmc.Material(name='B2O3')
rod.add_nuclide('B10',3.87448E-03)
rod.add_nuclide('B11',1.68447E-02)
rod.add_nuclide('O16',3.10787E-02)
rod.set_density('sum')
# ==================================================
# moderator PE1_8
m_pe1_8=openmc.Material(name='moderator PE1/8')
m_pe1_8.add_element('C',3.95860E-02)
m_pe1_8.add_element('H',7.77938E-02)
m_pe1_8.set_density('sum')
# ==================================================
# moderator PE10
m_pe10=openmc.Material(name='moderator PE10')
m_pe10.add_element('C',4.08960E-02)
m_pe10.add_element('H',7.97990E-02)
m_pe10.set_density('sum')
# ==================================================
# geometry
# region 
r_up1_l1_in=openmc.rectangular_prism(5.08,5.08,'z',    (-0.05,0.05))
r_up1_l1_gap=openmc.rectangular_prism(5.13,5.13,'z',  (-0.05,0.05) )
r_up1_l1_clad=openmc.rectangular_prism(5.43,5.43,'z',  (-0.05,0.05))

# region pb_bi
r_up1_l1_in_pi=openmc.rectangular_prism(5.00,5.00,'z',    (-0.05,0.05))
r_up1_l1_gap_pi=openmc.rectangular_prism(5.08,5.08,'z',  (-0.05,0.05) )
r_up1_l1_clad_pi=openmc.rectangular_prism(5.16,5.16,'z',  (-0.05,0.05))

#=======================================================

# control rod region
r=openmc.ZCylinder(r=1.69)
r1=openmc.ZCylinder(r=2)
r2=openmc.ZCylinder(r=2.25)
r3=openmc.ZCylinder(r=2.5)
r_0=openmc.ZPlane(z0=0)
r_1=openmc.ZPlane(z0=147)

rod_in=openmc.Cell(fill=rod,region=-r ) #& +r_0 & -r_1
rod_clad_in=openmc.Cell(fill=al,region=+r & -r1)
rod_gap=openmc.Cell(fill=None,region=+r1 & -r2)
rod_clad_out=openmc.Cell(fill=al,region=+r2 & -r3)
rod_out=openmc.Cell(fill=None,region=+r3)
u_rod=openmc.Universe(cells=[rod_in,rod_clad_in,rod_gap,rod_clad_out,rod_out])

#=======================================================
# Air
c_air=openmc.Cell(fill=air)
u1=openmc.Universe(name='air',cells=[c_air])

#  Al 
c_al=openmc.Cell(fill=al,region=r_up1_l1_in)
c_al_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_al_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_al_out=openmc.Cell(fill=air)
u_al=openmc.Universe(cells=[c_al,c_al_in,c_al_clad,c_al_out])

#  m_pe10
c_m_pe10=openmc.Cell(fill=m_pe10,region=r_up1_l1_in)
c_m_pe10_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_m_pe10_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_m_pe10_out=openmc.Cell(fill=air)
u_m_pe10=openmc.Universe(cells=[c_m_pe10,c_m_pe10_in,c_m_pe10_clad,c_m_pe10_out])

#  m_pe1_8
c_m_pe1_8=openmc.Cell(fill=m_pe1_8,region=r_up1_l1_in)
c_m_pe1_8_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_m_pe1_8_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_m_pe1_8_out=openmc.Cell(fill=air)
u_m_pe1_8=openmc.Universe(cells=[c_m_pe1_8,c_m_pe1_8_in,c_m_pe1_8_clad,c_m_pe1_8_out])

#  pe19
c_pe19=openmc.Cell(fill=pe19,region=r_up1_l1_in)
c_pe19_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_pe19_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_pe19_out=openmc.Cell(fill=air)
u_pe19=openmc.Universe(cells=[c_pe19,c_pe19_in,c_pe19_clad,c_pe19_out])


#  fuel
c_fuel=openmc.Cell(fill=fuel,region=r_up1_l1_in)
c_fuel_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_fuel_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_fuel_out=openmc.Cell(fill=air)
u_fuel=openmc.Universe(cells=[c_fuel,c_fuel_in,c_fuel_clad,c_fuel_out])

#  pb_bi
c_pi=openmc.Cell(fill=pb_bi,region=r_up1_l1_in_pi)
c_pi_in=openmc.Cell(fill=pb_bi_c1,region=~r_up1_l1_in_pi & r_up1_l1_gap_pi)
c_pi_clad=openmc.Cell(fill=pb_bi_c2,region=~ r_up1_l1_gap_pi & r_up1_l1_clad_pi)
c_pi_out=openmc.Cell(fill=air)
u_pi=openmc.Universe(cells=[c_pi,c_pi_in,c_pi_clad,c_pi_out])

#  F assembly
z_s_F=[0, 2, 24.5425 ,49.9425 ,89.9475 , 98.52 ,149.32, 152.40102]
z_list_F=[openmc.ZPlane(z0=i)for i in z_s_F]

lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_F=+z_list_F[0] & -z_list_F[1] & lat_box
r_low_F18=+z_list_F[1] & -z_list_F[2] & lat_box
r_low_F10=+z_list_F[2] & -z_list_F[3] & lat_box
r_fuel_F=+z_list_F[3] & -z_list_F[4] & lat_box
r_up_F18=+z_list_F[4] & -z_list_F[5] & lat_box
r_up_F10=+z_list_F[5] & -z_list_F[6] & lat_box
r_viod_F=+z_list_F[6] & -z_list_F[7] & lat_box

c_al_F=openmc.Cell(fill=al,region=r_al_F)

a=openmc.RectLattice()
a.pitch=(5.53,5.53,0.3175)
a.lower_left=(-5.53/2,-5.53/2,2)
a.universes=[[[u_m_pe1_8]]]*71

c_low_F18=openmc.Cell(fill=a,region=r_low_F18)

b=openmc.RectLattice()
b.pitch=(5.53,5.53,25.4)
b.universes=[[[u_m_pe10]]]
b.lower_left=(-5.53/2,-5.53/2,24.5425)

c_low_F10=openmc.Cell(fill=b,region=r_low_F10)

c=openmc.RectLattice()
c.pitch=(5.53,5.53,0.15875)
c.lower_left=(-5.53/2,-5.53/2,49.9425)
c.universes=[[[u_fuel]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]
            ,[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*36

c_fuel_F=openmc.Cell(fill=c,region=r_fuel_F)

d=openmc.RectLattice()
d.pitch=(5.53,5.53,0.3175)
d.lower_left=(-5.53/2,-5.53/2,89.9475)
d.universes=[[[u_m_pe1_8]]]*27

c_up_F18=openmc.Cell(fill=d,region=r_up_F18)

e=openmc.RectLattice()
e.pitch=(5.53,5.53,25.4)
e.lower_left=(-5.53/2,-5.53/2,98.52)
e.universes=[[[u_m_pe10]]]*2

c_up_F10=openmc.Cell(fill=e,region=r_up_F10)

c_void_F=openmc.Cell(fill=None,region=r_viod_F)
u_F=openmc.Universe(cells=[c_al_F,c_low_F18,c_low_F10,c_fuel_F,c_up_F18,c_up_F10,c_void_F])

# ================================================================
#  f assembly
z_s_f=[0, 2, 25.1775 ,50.5775, 60.1025 , 79.1525 , 88.6775 , 97.25 , 148.05 , 151.13102]
z_list_f=[openmc.ZPlane(z0=i)for i in z_s_f]
lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_f=+z_list_f[0] & -z_list_f[1] & lat_box
r_low_f18=+z_list_f[1] & -z_list_f[2] & lat_box
r_low_f10=+z_list_f[2] & -z_list_f[3] & lat_box
r_low_p  =+z_list_f[3] & -z_list_f[4] & lat_box
r_fuel_f=+z_list_f[4] & -z_list_f[5] & lat_box
r_up_p  =+z_list_f[5] & -z_list_f[6] & lat_box
r_up_f18=+z_list_f[6] & -z_list_f[7] & lat_box
r_up_f10=+z_list_f[7] & -z_list_f[8] & lat_box
r_viod_f=+z_list_f[8] & -z_list_f[9] & lat_box

c_al_f=openmc.Cell(fill=al,region=r_al_f)

aa=openmc.RectLattice()
aa.pitch=(5.53,5.53,0.3175)
aa.lower_left=(-5.53/2,-5.53/2,2)
aa.universes=[[[u_m_pe1_8]]]*73

c_low_f18=openmc.Cell(fill=aa,region=r_low_f18)

bb=openmc.RectLattice()
bb.pitch=(5.53,5.53,25.4)
bb.universes=[[[u_m_pe10]]]
bb.lower_left=(-5.53/2,-5.53/2,z_s_f[2])

c_low_f10=openmc.Cell(fill=bb,region=r_low_f10)

cc=openmc.RectLattice()
cc.pitch=(5.53,5.53,9.525)
cc.universes=[[[u_m_pe10]]]
cc.lower_left=(-5.53/2,-5.53/2,z_s_f[3])

c_low_p=openmc.Cell(fill=cc,region=r_low_p)


dd=openmc.RectLattice()
dd.pitch=(5.53,5.53,19.05)
dd.lower_left=(-5.53/2,-5.53/2,z_s_f[4])
dd.universes=[[[u_pi]]]

c_fuel_f=openmc.Cell(fill=dd,region=r_fuel_f)

ee=openmc.RectLattice()
ee.pitch=(5.53,5.53,9.525)
ee.lower_left=(-5.53/2,-5.53/2,z_s_f[5])
ee.universes=[[[u_m_pe10]]]

c_up_p=openmc.Cell(fill=ee,region=r_up_p)

ff=openmc.RectLattice()
ff.pitch=(5.53,5.53,0.3175)
ff.universes=[[[u_m_pe1_8]]]*27
ff.lower_left=(-5.53/2,-5.53/2,z_s_f[6])

c_up_f18=openmc.Cell(fill=ff,region=r_up_f18)

gg=openmc.RectLattice()
gg.pitch=(5.53,5.53,25.4)
gg.universes=[[[u_m_pe10]]]*2
gg.lower_left=(-5.53/2,-5.53/2,z_s_f[7])

c_up_f10=openmc.Cell(fill=gg,region=r_up_f10)

c_void_f=openmc.Cell(fill=None,region=r_viod_f)

u_f=openmc.Universe(cells=[c_al_f,c_low_f18,c_low_f10,c_low_p,c_fuel_f,c_up_p,c_up_f18,c_up_f10,c_void_f])


# ================================================================
#  16 assembly
z_s_16=[0, 2, 13.1125 , 35.655 , 61.055 , 78.835 , 87.4075 , 138.2075 , 149.32 , 152.40102]
z_list_16=[openmc.ZPlane(z0=i)for i in z_s_16]
lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_16    =+z_list_16[0] & -z_list_16[1] & lat_box
r_low_al_16=+z_list_16[1] & -z_list_16[2] & lat_box
r_low_1618 =+z_list_16[2] & -z_list_16[3] & lat_box
r_low_1610 =+z_list_16[3] & -z_list_16[4] & lat_box
r_fuel_16  =+z_list_16[4] & -z_list_16[5] & lat_box
r_up_1618  =+z_list_16[5] & -z_list_16[6] & lat_box
r_up_1610  =+z_list_16[6] & -z_list_16[7] & lat_box
r_up_al_16 =+z_list_16[7] & -z_list_16[8] & lat_box
r_viod_16  =+z_list_16[8] & -z_list_16[9] & lat_box

c_al_16=openmc.Cell(fill=al,region=r_al_16)

aaa=openmc.RectLattice()
aaa.pitch=(5.53,5.53,0.15875)
aaa.lower_left=(-5.53/2,-5.53/2,2)
aaa.universes=[[[u_al]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*10

c_low_al_16=openmc.Cell(fill=aaa,region=r_low_al_16)

bbb=openmc.RectLattice()
bbb.pitch=(5.53,5.53,0.3175)
bbb.universes=[[[u_m_pe1_8]]]*71
bbb.lower_left=(-5.53/2,-5.53/2,z_s_16[2])

c_low_1618=openmc.Cell(fill=bbb,region=r_low_1618)

ccc=openmc.RectLattice()
ccc.pitch=(5.53,5.53,25.4)
ccc.universes=[[[u_m_pe10]]]
ccc.lower_left=(-5.53/2,-5.53/2,z_s_16[3])

c_low_1610=openmc.Cell(fill=ccc,region=r_low_1610)

ddd=openmc.RectLattice()
ddd.pitch=(5.53,5.53,0.15875)
ddd.lower_left=(-5.53/2,-5.53/2,z_s_16[4])
ddd.universes=[[[u_fuel]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*16

c_fuel_16=openmc.Cell(fill=ddd,region=r_fuel_16)

eee=openmc.RectLattice()
eee.pitch=(5.53,5.53,0.3175)
eee.lower_left=(-5.53/2,-5.53/2,z_s_16[5])
eee.universes=[[[u_m_pe1_8]]]*27

c_up_1618=openmc.Cell(fill=eee,region=r_up_1618)

fff=openmc.RectLattice()
fff.pitch=(5.53,5.53,25.4)
fff.universes=[[[u_m_pe10]]]*2
fff.lower_left=(-5.53/2,-5.53/2,z_s_16[6])

c_up_1610=openmc.Cell(fill=fff,region=r_up_1610)

ggg=openmc.RectLattice()
ggg.pitch=(5.53,5.53,0.15875)
ggg.universes=[[[u_al]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*10
ggg.lower_left=(-5.53/2,-5.53/2,z_s_16[7])

c_up_al_16=openmc.Cell(fill=ggg,region=r_up_al_16)

c_void_16=openmc.Cell(fill=None,region=r_viod_16)

u_16=openmc.Universe(cells=[c_al_16,c_low_al_16,c_low_1618,c_low_1610,c_fuel_16,c_up_1618,c_up_1610,c_up_al_16,c_void_16])
# ============================================================================================================================
# A’-I left
l_ai_lat=openmc.RectLattice()
l_ai_lat.pitch=(5.53,5.53)
l_ai_lat.lower_left=(0,0)
l_ai_lat.universes=[[u_al]*4]*7

l_ai_box=openmc.rectangular_prism(5.53*4,5.53*7,'z',(5.53*4/2,5.53*7/2))
l_ai_cell=openmc.Cell(fill=l_ai_lat,region=l_ai_box)
# ============================================================================================================================
# J-P left
l_jp_lat=openmc.RectLattice()
l_jp_lat.pitch=(5.53,5.53)
l_jp_lat.lower_left=(0,39.01)
l_jp_lat.universes=[[u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_al,u_pe19,u_al],
                    [u_al,u_al,u_pe19,u_pe19],
                    [u_al]*4,
                    [u_al]*4,]

l_jp_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2,5.53*6/2+39.01))
l_jp_cell=openmc.Cell(fill=l_jp_lat,region=l_jp_box)
# ============================================================================================================================
# Q-T left
l_qt_lat=openmc.RectLattice()
l_qt_lat.pitch=(5.53,5.53)
l_qt_lat.lower_left=(0,72.49)
l_qt_lat.universes=[[u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_pe19,u_al  ,u_al  ],
                    [u_al,u_pe19,u_al  ,u_al  ]]

l_qt_box=openmc.rectangular_prism(5.53*4,5.53*3,'z',(5.53*4/2,5.53*3/2+72.49))
l_qt_cell=openmc.Cell(fill=l_qt_lat,region=l_qt_box)
# ============================================================================================================================
# U-Z left
l_uz_lat=openmc.RectLattice()
l_uz_lat.pitch=(5.53,5.53)
l_uz_lat.lower_left=(0,89.38)
l_uz_lat.universes=[[u_al]*4]*6

l_uz_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2,5.53*6/2+89.38))
l_uz_cell=openmc.Cell(fill=l_uz_lat,region=l_uz_box)
# ============================================================================================================================
# U-Z left1
l_uz1_lat=openmc.RectLattice()
l_uz1_lat.pitch=(5.53,5.53)
l_uz1_lat.lower_left=(22.42,89.38)
l_uz1_lat.universes=[[u_al]*6,
                     [u_pe19]*6,
                     [u_pe19,u_pe19,u_al,u_al,u_pe19,u_pe19],
                     [u_pe19,u_pe19,u_al,u_al,u_pe19,u_pe19],
                     [u_pe19]*6,
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10]]

l_uz1_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+22.42,5.53*6/2+89.38))
l_uz1_cell=openmc.Cell(fill=l_uz1_lat,region=l_uz1_box)
# ============================================================================================================================
# Q-T left1
l_qt1_lat=openmc.RectLattice()
l_qt1_lat.pitch=(5.53,5.53)
l_qt1_lat.lower_left=(22.42,72.49)
l_qt1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_rod   ,u_F]]

l_qt1_box=openmc.rectangular_prism(5.53*6,5.53*3,'z',(5.53*6/2+22.42,5.53*3/2+72.49))
l_qt1_cell=openmc.Cell(fill=l_qt1_lat,region=l_qt1_box)
# ============================================================================================================================
# J-P left1
l_jp1_lat=openmc.RectLattice()
l_jp1_lat.pitch=(5.53,5.53)
l_jp1_lat.lower_left=(22.42,39.01)
l_jp1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_rod   ,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_rod   ,u_F]]

l_jp1_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+22.42,5.53*6/2+39.01))
l_jp1_cell=openmc.Cell(fill=l_jp1_lat,region=l_jp1_box)
# ============================================================================================================================
# A’-I left1
l_ai1_lat=openmc.RectLattice()
l_ai1_lat.pitch=(5.53,5.53)
l_ai1_lat.lower_left=(22.42,0)
l_ai1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_al]*6]

l_ai1_box=openmc.rectangular_prism(5.53*6,5.53*7,'z',(5.53*6/2+22.42,5.53*7/2))
l_ai1_cell=openmc.Cell(fill=l_ai1_lat,region=l_ai1_box)
# ============================================================================================================================
# A’-I left2
u2=openmc.Universe(name='empty uni',cells=[openmc.Cell()])
l_ai2_lat=openmc.RectLattice()
l_ai2_lat.pitch=(5.53,5.53)
l_ai2_lat.lower_left=(55.9,0)
l_ai2_lat.universes=[[u_F,u_m_pe10,u_F],
                     [u_F,u_m_pe10,u_F],
                     [u_F,u_m_pe10,u_F],
                     [u_16,u_m_pe10,u_16],
                     [u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_al,u_m_pe10,u_al]]

l_ai2_box=openmc.rectangular_prism(5.53*3,5.53*7,'z',(5.53*3/2+55.9,5.53*7/2))
l_ai2_cell=openmc.Cell(fill=l_ai2_lat,region=l_ai2_box)
# ============================================================================================================================
# J-P left2
l_jp2_lat=openmc.RectLattice()
l_jp2_lat.pitch=(5.53,5.53)
l_jp2_lat.lower_left=(55.9,39.01)
l_jp2_lat.universes=[[u_F,u_F,u_F],
                     [u_f,u_f,u_f],
                     [u_f,u_f,u_f],
                     [u_f,u_m_pe10 ,u_f],
                     [u_F,u_m_pe10 ,u_F],
                     [u_F,u_m_pe10 ,u_F]]

l_jp2_box=openmc.rectangular_prism(5.53*3,5.53*6,'z',(5.53*3/2+55.9,5.53*6/2+39.01))
l_jp2_cell=openmc.Cell(fill=l_jp2_lat,region=l_jp2_box)
# ============================================================================================================================
# Q-T left2
l_qt2_lat=openmc.RectLattice()
l_qt2_lat.pitch=(5.53,5.53)
l_qt2_lat.lower_left=(55.9,72.49)
l_qt2_lat.universes=[[u_pe19,u_pe19,u_pe19],
                     [u_pe19,u_pe19,u_pe19],
                     [u_F   ,u_F   ,u_F   ]]

l_qt2_box=openmc.rectangular_prism(5.53*3,5.53*3,'z',(5.53*3/2+55.9,5.53*3/2+72.49))
l_qt2_cell=openmc.Cell(fill=l_qt2_lat,region=l_qt2_box)
# ===========================================================================================================================
# U-Z left2
l_uz2_lat=openmc.RectLattice()
l_uz2_lat.pitch=(5.53,5.53)
l_uz2_lat.lower_left=(55.9,89.38)
l_uz2_lat.universes=[[u_al]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_m_pe10]*3]

l_uz2_box=openmc.rectangular_prism(5.53*3,5.53*6,'z',(5.53*3/2+55.9,5.53*6/2+89.38))
l_uz2_cell=openmc.Cell(fill=l_uz2_lat,region=l_uz2_box)
# ===========================================================================================================================
# U-Z left3
l_uz3_lat=openmc.RectLattice()
l_uz3_lat.pitch=(5.53,5.53)
l_uz3_lat.lower_left=(72.79,89.38)
l_uz3_lat.universes=[[u_al]*6,
                     [u_pe19,u_pe19,u_pe19,u_pe19,u_pe19,u_al],
                     [u_pe19,u_pe19,u_al  , u_al ,u_pe19,u_al],
                     [u_pe19,u_pe19,u_al  , u_al ,u_pe19,u_al],
                     [u_pe19,u_pe19,u_pe19,u_pe19,u_pe19,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_al]]

l_uz3_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+72.79,5.53*6/2+89.38))
l_uz3_cell=openmc.Cell(fill=l_uz3_lat,region=l_uz3_box)
# ============================================================================================================================
# Q-T left3
l_qt3_lat=openmc.RectLattice()
l_qt3_lat.pitch=(5.53,5.53)
l_qt3_lat.lower_left=(72.79,72.49)
l_qt3_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F     ,u_rod   ,u_m_pe10,u_m_pe10,u_pe19,u_pe19]]

l_qt3_box=openmc.rectangular_prism(5.53*6,5.53*3,'z',(5.53*6/2+72.79,5.53*3/2+72.49))
l_qt3_cell=openmc.Cell(fill=l_qt3_lat,region=l_qt3_box)
# ============================================================================================================================
# J-P left3
l_jp3_lat=openmc.RectLattice()
l_jp3_lat.pitch=(5.53,5.53)
l_jp3_lat.lower_left=(72.79,39.01)
l_jp3_lat.universes=[[u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_rod,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_rod,u_m_pe10,u_m_pe10,u_pe19,u_pe19]]

l_jp3_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+72.79,5.53*6/2+39.01))
l_jp3_cell=openmc.Cell(fill=l_jp3_lat,region=l_jp3_box)
# ============================================================================================================================
# A’-I left3
l_ai3_lat=openmc.RectLattice()
l_ai3_lat.pitch=(5.53,5.53)
l_ai3_lat.lower_left=(72.79,0)
l_ai3_lat.universes=[[u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_al,u_al,u_al,u_al,u_pe19,u_pe19]]

l_ai3_box=openmc.rectangular_prism(5.53*6,5.53*7,'z',(5.53*6/2+72.79,5.53*7/2))
l_ai3_cell=openmc.Cell(fill=l_ai3_lat,region=l_ai3_box)
# ============================================================================================================================
# A’-I left3
l_ai4_lat=openmc.RectLattice()
l_ai4_lat.pitch=(5.53,5.53)
l_ai4_lat.lower_left=(106.27,0)
l_ai4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al]]

l_ai4_box=openmc.rectangular_prism(5.53*4,5.53*7,'z',(5.53*4/2+106.27,5.53*7/2))
l_ai4_cell=openmc.Cell(fill=l_ai4_lat,region=l_ai4_box)
# ============================================================================================================================
# J-P left4
l_jp4_lat=openmc.RectLattice()
l_jp4_lat.pitch=(5.53,5.53)
l_jp4_lat.lower_left=(106.27,39.01)
l_jp4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al]]

l_jp4_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2+106.27,5.53*6/2+39.01))
l_jp4_cell=openmc.Cell(fill=l_jp4_lat,region=l_jp4_box)
# ============================================================================================================================
# Q-T left4
l_qt4_lat=openmc.RectLattice()
l_qt4_lat.pitch=(5.53,5.53)
l_qt4_lat.lower_left=(106.27,72.49)
l_qt4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al]]

l_qt4_box=openmc.rectangular_prism(5.53*4,5.53*3,'z',(5.53*4/2+106.27,5.53*3/2+72.49))
l_qt4_cell=openmc.Cell(fill=l_qt4_lat,region=l_qt4_box)
# ===========================================================================================================================
# U-Z left4
l_uz4_lat=openmc.RectLattice()
l_uz4_lat.pitch=(5.53,5.53)
l_uz4_lat.lower_left=(106.27,89.38)
l_uz4_lat.universes=[[u_al]*4]*6

l_uz4_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2+106.27,5.53*6/2+89.38))
l_uz4_cell=openmc.Cell(fill=l_uz4_lat,region=l_uz4_box)

# ============================================================================================================================
box=openmc.rectangular_prism(150,150,'z',(150/2,150/2),'vacuum')
z1=openmc.ZPlane(z0=0)
z1.boundary_type='vacuum'
z2=openmc.ZPlane(z0=150)
z2.boundary_type='vacuum'
box_region= +z1& -z2 & box
all_cell=openmc.Cell(fill=air)
cells=[l_ai_cell,l_jp_cell,l_qt_cell,l_uz_cell,l_uz1_cell,l_qt1_cell,l_jp1_cell,l_ai1_cell,l_ai2_cell,l_jp2_cell,l_qt2_cell,l_uz2_cell
        ,l_uz3_cell,l_qt3_cell,l_jp3_cell,l_ai3_cell,l_ai4_cell,l_jp4_cell,l_qt4_cell,l_uz4_cell,all_cell]
u_all=openmc.Universe(cells=cells)
cell_all=openmc.Cell(fill=u_all,region=box_region)
geom = openmc.Geometry([cell_all])
geom.export_to_xml()

mats = list(geom.get_all_materials().values())
openmc.Materials(mats).export_to_xml()
point = openmc.stats.Point((61, 27, 51))
src = openmc.Source(space=point)
settings = openmc.Settings()
settings.source = src
settings.batches = 200
settings.inactive = 10
settings.particles = 1000
# batches = 100
# inactive = 50
# particles = 1000
# Instantiate a Settings object
# settings_file = openmc.Settings()
# settings_file.batches = batches
# settings_file.inactive = inactive
# settings_file.particles = particles
# # Create an initial uniform spatial source distribution over fissionable zones
# bounds = [55.3, 22.12, 50, 72, 77, 88]
# uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
# settings_file.source = openmc.Source(space=uniform_dist)
settings.export_to_xml()

openmc.run(geometry_debug=1)
p = openmc.Plot()
p.filename = 'pinplot'
p.width = (10, 10)
p.pixels = (2000, 2000)
p.color_by = 'cell'
# p.colors={fuel:'blue',m_pe10:'green',m_pe10:'red',air:'black',al:'yellow'}
# p.colors={fuel:'blue'}
p.basis='xy'
p.origin=(55.9,0,100)
# p.origin=(10*5.53+0.3,4*5.53+0.3,50.0425)
plots = openmc.Plots([p])
plots.export_to_xml()
openmc.plot_geometry()

And the k result only have 0.24544,leakage fraction have 0.83594,which is far from the benckmark.


I tried to find out if there was a geometric error, compared the Plot with the benchmark I have not find any differences.But when I use the debug mode ,the system alerts me that there have detected overlapping cells.

I’ve been checking for a long time and can’t find my mistake :frowning: , so if anyone would like to help me find my mistake or suggest where I might have made it I’d really appreciate it!Thanks in advance!

Hi @world3,

Neat benchmark model! Unfortunately, I don’t have time to look through the Python code anytime soon, but I wanted to point you to the interactive plotter (GitHub - openmc-dev/plotter: Native plotting GUI for model design and verification), which might help you to debug some of the problems you’re having with the geometry. It should make it easier to determine the regions in the model with overlaps by turning on the “overlaps” option in the color dialog box or using the Ctrl + P keyboard shortcut.

I hope this helps you out!

Best,

Patrick

Thank you very much! It’s help me alot :smiley: