Merge branch '1k_openmc_model' into 'master'

1k openmc model See merge request !3

Merge branch '1k_openmc_model' into 'master'
1k openmc model See merge request !3
fb643297 · Elia Merzari · 5466364b · c71cea29 · fb643297 · fb643297
Commit fb643297 authored Oct 28, 2020 by Elia Merzari
--- a/146_pebble/geometry.xml
+++ b/146_pebble/geometry.xml
--- a/146_pebble/master.i
+++ b/146_pebble/master.i
@@ -121,7 +121,7 @@
 [Outputs]
  exodus = true
  interval = 100
-  file_base = ../../../../../../../../scratch/gastdr/coe/146_pebbles/out
+  file_base = ./
 []

 [MultiApps]

--- a/146_pebble/openmc.i
+++ b/146_pebble/openmc.i
@@ -27,7 +27,9 @@

 [Problem]
  type = OpenMCProblem
-  power = 146000 # 1000 per pebble
+  power = 1460 # 1000 per pebble
+  tally_type = 'cell'
+  mesh_template = 'sphere.e'
  centers = '-8.065973074285714972e+00 2.236820423458646623e+00 1.168503152842105308e+01
 -4.108761094736841812e+00 7.264663405714286704e+00 1.148013507969924873e+01
 -8.190250614135338836e+00 -1.669100062556391206e+00 1.174674951699248204e+01

--- a/146_pebble/openmc_master.i
+++ b/146_pebble/openmc_master.i
@@ -237,10 +237,10 @@

 [Transfers]
  [./heat_source_from_openmc]
-    type = MultiAppUserObjectTransfer
+    type = MultiAppCopyTransfer
    direction = from_multiapp
    multi_app = openmc
-    user_object = heat_source
+    source_variable = heat_source
    variable = heat_source
  []
  [./average_temp_to_openmc]

--- a/146_pebble/plots.xml
+++ b/146_pebble/plots.xml
 <?xml version='1.0' encoding='utf-8'?>
 <plots>
-  <plot basis="xy" color_by="material" filename="plot1" id="1" type="slice">
+  <plot color_by="material" filename="plotv1" id="5" type="voxel">
    <origin>0 0 0</origin>
-    <width>30.25 30.25</width>
-    <pixels>1000 1000</pixels>
-    <color id="1" rgb="165 42 42" />
-    <color id="2" rgb="176 196 222" />
-    <color id="3" rgb="0 0 255" />
-    <color id="4" rgb="255 165 0" />
-    <color id="5" rgb="0 255 255" />
-    <color id="6" rgb="0 191 255" />
-    <color id="7" rgb="0 0 128" />
-    <color id="8" rgb="255 255 0" />
-  </plot>
-  <plot basis="xy" color_by="material" filename="plot2" id="2" type="slice">
-    <origin>0 0 1</origin>
-    <width>30.25 30.25</width>
-    <pixels>1000 1000</pixels>
-    <color id="1" rgb="165 42 42" />
-    <color id="2" rgb="176 196 222" />
-    <color id="3" rgb="0 0 255" />
-    <color id="4" rgb="255 165 0" />
-    <color id="5" rgb="0 255 255" />
-    <color id="6" rgb="0 191 255" />
-    <color id="7" rgb="0 0 128" />
-    <color id="8" rgb="255 255 0" />
-  </plot>
-  <plot basis="xy" color_by="material" filename="plot3" id="3" type="slice">
-    <origin>0 0 1</origin>
-    <width>30.25 30.25</width>
-    <pixels>1000 1000</pixels>
-    <color id="1" rgb="165 42 42" />
-    <color id="2" rgb="176 196 222" />
-    <color id="3" rgb="0 0 255" />
-    <color id="4" rgb="255 165 0" />
-    <color id="5" rgb="0 255 255" />
-    <color id="6" rgb="0 191 255" />
-    <color id="7" rgb="0 0 128" />
-    <color id="8" rgb="255 255 0" />
-  </plot>
-  <plot basis="xz" color_by="material" filename="plot4" id="4" type="slice">
-    <origin>0 0 5.08</origin>
-    <width>31.49 31.49</width>
-    <pixels>1000 1000</pixels>
+    <width>20 20 20</width>
+    <pixels>400 400 400</pixels>
    <color id="1" rgb="165 42 42" />
    <color id="2" rgb="176 196 222" />
    <color id="3" rgb="0 0 255" />

--- a/146_pebble/settings.xml
+++ b/146_pebble/settings.xml
@@ -5,7 +5,7 @@
  <batches>150</batches>
  <inactive>50</inactive>
  <source strength="1.0">
-    <space type="fission">
+    <space type="box">
      <parameters>-15.125 -15.125 -10.665 15.125 15.125 20.825</parameters>
    </space>
  </source>
@@ -18,7 +18,7 @@
    <domain_type>material</domain_type>
    <domain_ids>1</domain_ids>
    <samples>10000000</samples>
-    <lower_left>-1.4 -1.4 -1.4</lower_left>
-    <upper_right>1.4 1.4 1.4</upper_right>
+    <lower_left>-10 -10 -10</lower_left>
+    <upper_right>10 10 10</upper_right>
  </volume_calc>
 </settings>
--- a/146_pebble/sphere.e
+++ b/146_pebble/sphere.e
--- a/1568_pebble_yuhsiang_E524k/NEKMOOSE
+++ b/1568_pebble_yuhsiang_E524k/NEKMOOSE
+      parameter (lsurf_m=170000)
+      common /point_cloudx/ pc_x(lsurf_m*4)
+      common /point_cloudy/ pc_y(lsurf_m*4)
+      common /point_cloudz/ pc_z(lsurf_m*4)
+      common /point_cloudt/ pc_t(lsurf_m*4)
+      common /point_cloudf/ pc_f(lsurf_m*4)
+      common /ref_element/  pc_flag(lsurf_m*4)
+      common /tot_surf/ nw_bdt       
+      common /data_reference/ plmask(lelt), pflmask(lelt), 
+     & pplist(lelt), offset, nw_bd1
+      common /test_passing/ flux_moose, temp_nek
+      real pc_x, pc_y, pc_z, pc_flag, pc_t, pc_f
+      integer plmask, pflmask, pllist, offset, nw_bd1, nw_bdt 
+
--- a/1568_pebble_yuhsiang_E524k/cardinal.i
+++ b/1568_pebble_yuhsiang_E524k/cardinal.i
--- a/1568_pebble_yuhsiang_E524k/core_1568_pebbles.py
+++ b/1568_pebble_yuhsiang_E524k/core_1568_pebbles.py
+# -------------- Legend ---------------------
+# b_  = box
+# c_  = cell
+# cz_ = cylinder z
+# l_  = lattice
+# m_  = material
+# px_ = plane x
+# py_ = plane y
+# pz_ = plane z
+# r_  = region
+# s_  = sphere
+# u_  = universe
+# -------------- Modules --------------------
+import os
+import math
+import numpy as np
+import openmc
+import openmc.model
+import scipy.stats
+import pandas as pd
+from   math   import pi
+from   IPython.display import Image
+import matplotlib.pyplot as plt
+# -------------- Functions ------------------
+def density_flibe(p, T):
+    # assumes the default value for drho_dp of 1.7324e-7
+    rho = -0.4884*T + 1.7324e-7*(p - 101325.0) + 2413.0 # kg/m3
+    return ('kg/m3', rho)
+# -------------- Arbitrary Parameters -------
+debug                 = 0              # Set 1 to print debug info
+voxel                 = 0              # Set 1 to generate 3d voxel plots              and 0 for 2d ppm plots
+random_distribution   = 0              # Set 1 for TRISO particles random distribution and 0 for regular lattice
+# -------------- Units Parameters -----------
+mm                    = 0.1   # default cm
+μm                    = 1e-4  # default cm
+# -------------- Geometry Parameters --------
+# TRISO particle
+radius_fuel           = 400.0*μm/2     # UCBTH-14-002, Table 2-1
+thickness_c_buffer    = 100.0*μm       # UCBTH-14-002, Table 2-1
+thickness_pyc_inner   =  35.0*μm       # UCBTH-14-002, Table 2-1
+thickness_sic         =  35.0*μm       # UCBTH-14-002, Table 2-1
+thickness_pyc_outer   =  35.0*μm       # UCBTH-14-002, Table 2-1
+radius_c_buffer       = radius_fuel      + thickness_c_buffer
+radius_pyc_inner      = radius_c_buffer  + thickness_pyc_inner
+radius_sic            = radius_pyc_inner + thickness_sic
+radius_pyc_outer      = radius_sic       + thickness_pyc_outer
+packing_fraction      = 0.40           # UCBTH-14-002, Table 2-1
+if random_distribution==0:
+   pitch_triso_lattice   = (radius_fuel+thickness_c_buffer+thickness_pyc_inner+thickness_sic+thickness_pyc_outer)*(4*pi/3/packing_fraction)**0.33333333333333333333333
+# Pebble - inner zone graphite; central active zone triso; outer zone graphite
+# hyuan@anl.gov email April 2, 2019 8:40 PM
+# In the actual TAMU experiment, the pebble diameter is 22.225mm.
+# But in CAD, we made it 22.225mm*0.95 = 21.11375 mm To avoid any pebble contacts.
+# So the pebble radius is 21.11375 mm/2 = 10.556875 mm = 1.0556875 cm
+# To scale up tp 3 cm pebble radius ,  I should scale Nek mesh and CAD file by a factor
+# 3/1.0556875 = 2.8417500444
+tamu_exp_factor       = 1/1.0556875               # Factor to scale geometry to TAMU experiment with pebble radius 1.0556875 cm
+radius_pebble_inner   = 2.5/2                     # UCBTH-14-002, Table 2-1 (differs slightly from Cisneros, Table 5-2) ; scaled by TAMU experiment factor
+radius_pebble_outer   = 3.0/2                     # UCBTH-14-002, Table 2-1; Cisneros, Table 5-2                        ; scaled by TAMU experiment factor
+radius_pebble_central = radius_pebble_outer - 0.1 # UCBTH-14-002, Table 2-1; Cisneros, Table 5-2      ; scaled by TAMU experiment factor
+tolerance             = 0.00                      # Tolerance for pebbles universe filling
+radius_pebble_flibe   = radius_pebble_outer+tolerance
+# Pebble centers coordinates (x,y,z)
+print("Reading pseudo-random pebble centers from file pebble_centers.txt")
+pebble_centers        = np.loadtxt('list_pebbles_1568.csv', delimiter=',', skiprows=1)[:, 1:]
+print("File pebble_centers.txt reading completed")
+rescaled_file = 'pebble_centers.txt'
+np.savetxt(rescaled_file, pebble_centers.reshape((-1,3)))
+print("Saved rescaled centers to " + rescaled_file)
+# Vessel
+extra_thickness       =  12.125  # change this value to set keff=1
+x_vessel              =  0.0  # x position vessel
+y_vessel              =  0.0  # y position vessel
+z1_vessel = np.min(pebble_centers[:, 2]) - radius_pebble_outer
+z2_vessel = np.max(pebble_centers[:, 2]) + radius_pebble_outer
+radius_vessel = np.max(np.linalg.norm(pebble_centers[:, :-1], axis=1)) + radius_pebble_outer
+z1_vessel             -= extra_thickness # bottom plane vessel ; scaled by TAMU experiment factor
+z2_vessel             += extra_thickness # top plane    vessel ; scaled by TAMU experiment factor
+radius_vessel         += extra_thickness #                       scaled by TAMU experiment factor
+height_vessel         =  z2_vessel - z1_vessel
+# -------------- Printing Parameters ---------
+print ("GEOMETRY PARAMETERS")
+print ("TRISO particles radius_fuel           [cm] = ", radius_fuel)
+print ("TRISO particles radius_c_buffer       [cm] = ", radius_c_buffer)
+print ("TRISO particles radius_pyc_inner      [cm] = ", radius_pyc_inner)
+print ("TRISO particles radius_sic            [cm] = ", radius_sic)
+print ("TRISO particles radius_pyc_outer      [cm] = ", radius_pyc_outer)
+print ("TRISO particles packing fraction           = ", packing_fraction)
+if random_distribution==0:
+   print ("TRISO particles regular lattice pitch_triso_lattice [cm] = ", pitch_triso_lattice)
+print ("Pebble radius_pebble_inner   [cm] = ", radius_pebble_inner)
+print ("Pebble radius_pebble_outer   [cm] = ", radius_pebble_outer)
+print ("Pebble radius_pebble_central [cm] = ", radius_pebble_central)
+print ("Pebble radius_pebble_flibe   [cm] = ", radius_pebble_flibe)
+print ("Vessel radius_vessel         [cm] = ", radius_vessel)
+print ("Vessel z1_vessel (min z)     [cm] = ", z1_vessel)
+print ("Vessel height_vessel         [cm] = ", height_vessel)
+print ("Pebbles centers [cm] :")
+print (pebble_centers)
+# -------------- Materials Parameters --------
+# TRISO particle
+enrichment_uranium   = 0.199
+fraction_u234        = 2e-3                 # U234/(U234+U235) mass from Giacint facility fuel
+density_fuel         = ('g/cm3', 10.5)      # UCBTH-14-002, Table 2-1
+density_c_buffer     = ('g/cm3', 1.0)       # Cisneros    , Table 3-1
+density_pyc          = ('g/cm3', 1.87)      # Cisneros    , Table 3-1
+density_sic          = ('g/cm3', 3.2)       # Cisneros    , Table 3-1
+# Pebble
+density_graphite_inner = ('g/cm3', 1.59368) # Cisneros    , Table 3-1
+density_graphite_outer = ('g/cm3', 1.74)    # Cisneros    , Table 3-1
+# FLiBe coolant
+enrichment_li7       = 0.99995
+temperature_inlet    = 273.15 + 600.0       # UCBTH-14-002, Table 1-1
+temperature_outlet   = 273.15 + 700.0       # UCBTH-14-002, Table 1-1
+temperature_flibe    = (temperature_inlet+temperature_outlet)/2
+# --------------------------------------------------
+# NO PARAMETERS DEFINITION BEYOND THIS LINE
+# -------------- Materials Definition --------------
+# TRISO particle
+# Fuel from Nagley et al. Fabrication of Uranium Oxycarbide Kernels for HTR Fuel https://inldigitallibrary.inl.gov/sites/sti/sti/4886646.pdf Table 2
+m_fuel = openmc.Material(name='m_fuel - uranium oxycarbide - triso partciles')
+m_fuel.add_nuclide('U234', 89.58*   enrichment_uranium*   fraction_u234 , percent_type='wo')
+m_fuel.add_nuclide('U235', 89.58*   enrichment_uranium*(1-fraction_u234), percent_type='wo')
+m_fuel.add_nuclide('U238', 89.58*(1-enrichment_uranium)                 , percent_type='wo')
+m_fuel.add_nuclide('C0'   ,  1.80                                        , percent_type='wo')
+m_fuel.add_element('O'   ,  8.62                                        , percent_type='wo')
+m_fuel.set_density(*density_fuel)
+#
+m_graphite_c_buffer = openmc.Material(name='m_graphite_c_buffer - triso partciles')
+m_graphite_c_buffer.set_density(*density_c_buffer)
+m_graphite_c_buffer.add_nuclide('C0', 1.0)
+m_graphite_c_buffer.add_s_alpha_beta('c_Graphite')
+#
+m_graphite_pyc = openmc.Material(name='m_graphite_pyc - triso particles')
+m_graphite_pyc.set_density(*density_pyc)
+m_graphite_pyc.add_nuclide('C0', 1.0)
+m_graphite_pyc.add_s_alpha_beta('c_Graphite')
+#
+m_sic = openmc.Material(name='sic - triso partciles')
+m_sic.add_nuclide('C0' , 1.0)
+m_sic.add_element('Si', 1.0)
+m_sic.set_density(*density_sic)
+# TODO: Add SiC S(alpha, beta) ; ENDF/B-VIIIb4 has data for both Si and C in SiC
+#
+m_graphite_matrix = openmc.Material(name='m_graphite_matrix - triso particles')
+m_graphite_matrix.set_density('g/cm3', 1.6)
+m_graphite_matrix.add_nuclide('C0', 1.0)
+m_graphite_matrix.add_s_alpha_beta('c_Graphite')
+# TODO: Need real specifications for the carbon filler
+# Pebble
+m_graphite_inner = openmc.Material(name='m_graphite_inner - pebble inner zone')
+m_graphite_inner.set_density(*density_graphite_inner)
+m_graphite_inner.add_nuclide('C0', 1.0)
+m_graphite_inner.add_s_alpha_beta('c_Graphite')
+#
+m_graphite_outer = openmc.Material(name='m_graphite_outer - pebble outer zone')
+m_graphite_outer.set_density(*density_graphite_outer)
+m_graphite_outer.add_nuclide('C0', 1.0)
+m_graphite_outer.add_s_alpha_beta('c_Graphite')
+# FLiBe coolant - From Cisneros, appendix B, material 24
+m_flibe = openmc.Material(name='m_flibe - 2LiF-BeF2')
+m_flibe.set_density(*density_flibe(101.325e3, temperature_flibe))
+m_flibe.add_nuclide('Li7', 2.0*     enrichment_li7)
+m_flibe.add_nuclide('Li6', 2.0*(1 - enrichment_li7))
+m_flibe.add_element('Be' , 1.0)
+m_flibe.add_element('F'  , 4.0)
+# TODO: FLiBe coolant - no S(alpha, beta) data available up to ENDF/B-VIIIb4
+# -------------- Geometry Definition --------------
+# TRISO particle
+s_fuel             = openmc.Sphere(R=radius_fuel)
+s_c_buffer         = openmc.Sphere(R=radius_c_buffer)
+s_pyc_inner        = openmc.Sphere(R=radius_pyc_inner)
+s_sic              = openmc.Sphere(R=radius_sic)
+s_pyc_outer        = openmc.Sphere(R=radius_pyc_outer)
+c_triso_fuel       = openmc.Cell(name='c_triso_fuel'     , fill=m_fuel,              region=-s_fuel)
+c_triso_c_buffer   = openmc.Cell(name='c_triso_c_buffer' , fill=m_graphite_c_buffer, region=+s_fuel      & -s_c_buffer)
+c_triso_pyc_inner  = openmc.Cell(name='c_triso_pyc_inner', fill=m_graphite_pyc,      region=+s_c_buffer  & -s_pyc_inner)
+c_triso_sic        = openmc.Cell(name='c_triso_sic'      , fill=m_sic,               region=+s_pyc_inner & -s_sic)
+c_triso_pyc_outer  = openmc.Cell(name='c_triso_pyc_outer', fill=m_graphite_pyc,      region=+s_sic       & -s_pyc_outer)
+c_triso_matrix     = openmc.Cell(name='c_triso_matrix'   , fill=m_graphite_matrix,   region=+s_pyc_outer)
+u_triso            = openmc.Universe(cells=[c_triso_fuel, c_triso_c_buffer, c_triso_pyc_inner, c_triso_sic, c_triso_pyc_outer, c_triso_matrix])
+# Pebble
+s_pebble_inner     = openmc.Sphere(R=radius_pebble_inner)
+s_pebble_central   = openmc.Sphere(R=radius_pebble_central)
+s_pebble_outer     = openmc.Sphere(R=radius_pebble_outer)
+c_pebble_inner     = openmc.Cell(name='c_pebble_inner'  , fill=m_graphite_inner, region=-s_pebble_inner)
+c_pebble_central   = openmc.Cell(name='c_pebble_central',                        region=+s_pebble_inner   & -s_pebble_central)
+c_pebble_outer     = openmc.Cell(name='c_pebble_inner'  , fill=m_graphite_outer, region=+s_pebble_central & -s_pebble_outer)
+c_pebble_flibe     = openmc.Cell(name='c_flibe'         , fill=m_flibe         , region=+s_pebble_outer)
+# Fill c_pebble_central with TRISO particles
+if random_distribution==0:
+# TRISO particles regular lattice using 'pitch_triso_lattice'
+      l_triso                    = openmc.RectLattice(name='l_triso')
+      l_triso.lower_left         = (-pitch_triso_lattice/2, -pitch_triso_lattice/2, -pitch_triso_lattice/2)
+      l_triso.pitch              = ( pitch_triso_lattice  ,  pitch_triso_lattice  ,  pitch_triso_lattice)
+      l_triso.outer              = u_triso
+      l_triso.universes          = np.tile(u_triso, (1,1,1))
+#     l_triso.universes[0, 0, 0] = universe_triso
+else:
+# TRISO particles random distribution using 'packing_fraction'
+      c_triso_pyc_outer_random   = openmc.Cell(name='c_triso_pyc_outer_random', fill=m_graphite_pyc, region=+s_sic)
+      u_triso_random             = openmc.Universe(cells=[c_triso_fuel, c_triso_c_buffer, c_triso_pyc_inner, c_triso_sic, c_triso_pyc_outer_random])
+      if not os.path.exists('triso_centers.npy'):
+            print("Writing random TRISO centers to file triso_centers.npy")
+            r_triso_random             = +s_pebble_inner & -s_pebble_central
+            spheres_random             = openmc.model.pack_spheres(radius=radius_pyc_outer, region=r_triso_random, pf=packing_fraction, initial_pf=0.15)
+            triso_random               = [openmc.model.TRISO(radius_pyc_outer, u_triso_random, i) for i in spheres_random]
+            centers = np.vstack([i.center for i in triso_random])
+            if debug==1:
+               print(triso_random[0])
+               print(triso_random[1])
+               print(centers.min(axis=0))
+               print(centers.max(axis=0))
+            print("Calculated packing fraction of the TRISO particles random distribution = ",len(triso_random)*radius_pyc_outer**3/(radius_pebble_central**3-radius_pebble_inner**3))
+            np.save('triso_centers.npy', triso_random)
+      print("Reading random TRISO centers from file triso_centers.npy")
+      triso_random = np.load('triso_centers.npy')
+      lower_left, upper_right    = c_pebble_central.region.bounding_box
+      shape                      = (5, 5, 5)
+      pitch                      = (upper_right - lower_left)/shape
+      l_triso                    = openmc.model.create_triso_lattice(triso_random, lower_left, pitch, shape, m_graphite_matrix)
+      print("File triso_centers.npy reading completed")
+c_pebble_central.fill   = l_triso
+u_pebble                = openmc.Universe(cells=[c_pebble_inner, c_pebble_central, c_pebble_outer, c_pebble_flibe])
+# Vessel
+cz_vessel               = openmc.ZCylinder(x0=x_vessel, y0=y_vessel, R=radius_vessel, boundary_type='reflective') # transmission | reflective | vacuum
+pz_vessel_bottom        = openmc.ZPlane(z0=z1_vessel, boundary_type='reflective')
+pz_vessel_top           = openmc.ZPlane(z0=z2_vessel, boundary_type='reflective')
+cell_name = ['cell_pebble_' + str(i) for i in range (len(pebble_centers))]
+s_pebble = [openmc.Sphere(x0=pebble_centers[i,0], y0=pebble_centers[i,1], z0=pebble_centers[i,2], R=radius_pebble_flibe, boundary_type='transmission') for i in range (len(pebble_centers))]
+c_pebble = [openmc.Cell(name=cell_name[i], fill=u_pebble, region=-s_pebble[i]) for i in range (len(pebble_centers))]
+r_vessel                = -cz_vessel & +pz_vessel_bottom & -pz_vessel_top
+for i in range(len(pebble_centers)):
+    c_pebble[i].translation = pebble_centers[i]
+    r_vessel                = r_vessel & +s_pebble[i]
+c_vessel                = openmc.Cell(name='c_vessel', fill=m_flibe, region = r_vessel)
+#
+#  Pebbles pseudo-random distribution using openmc.model.TRISO function - IT DOES NOT WORK WITH TALLY CARDS
+#  c_pebble_outer_random   = openmc.Cell(name='c_pebble_outer_random', fill=m_graphite_outer, region=+s_pebble_outer)
+#  u_pebble_random         = openmc.Universe(cells=[c_pebble_inner, c_pebble_central, c_pebble_outer])
+#  print("Reading random pebble centers from file pebble_centers.txt")
+#  spheres_random          = np.loadtxt('pebble_centers.txt')
+#  spheres_random          = spheres_random*radius_pebble_outer # pebbles centers have cm unit and must be scaled by TAMU experiment factor
+#  print("File pebble_centers.txt reading completed")
+#  print(spheres_random)
+#  pebble_random           = [openmc.model.TRISO(radius_pebble_outer, u_pebble_random, i) for i in spheres_random]
+#  lower_left, upper_right = c_vessel.region.bounding_box
+#  shape                   = (5, 5, 5)
+#  pitch                   = (upper_right - lower_left)/shape
+#  l_pebble                = openmc.model.create_triso_lattice(pebble_random, lower_left, pitch, shape, m_flibe)
+#  c_vessel.fill           = l_pebble
+#
+# Global geometry
+# c_all = [[c_vessel, c_pebble[i]]  for i in range (len(pebble_centers))] # The for loop does not work
+u_zero                  = openmc.Geometry([c_vessel] + list(c_pebble))
+# -------------- Printing Cells --------------
+print ("Cell of the vessel:")
+print (c_vessel)
+print ("Cells of the pebbles:")
+print (c_pebble)
+# -------------- Settings --------------------
+settings           = openmc.Settings()
+settings.source    = openmc.Source(space=openmc.stats.Box(*c_vessel.bounding_box, only_fissionable=True))
+settings.particles = 10000
+settings.inactive  = 50
+settings.batches   = 150
+settings.temperature = dict(default=573, method='interpolation',
+                            multipole=True, range=(300.0, 1500.0), tolerance=1000.0)
+# Fuel volume calculation
+volume_fuel                  = openmc.VolumeCalculation([m_fuel], 10000000, *c_pebble_central.bounding_box)
+settings.volume_calculations = [volume_fuel]
+model              = openmc.model.Model(geometry=u_zero, settings=settings)
+model.export_to_xml()
+# -------------- Tally -----------------------
+# print(l_pebble)
+# print(c_pebble_central)
+# print(u_pebble_random)
+# print(pebble_random)
+# # filter1            = openmc.DistribcellFilter([13])
+# filter1            = openmc.UniverseFilter([u_pebble_random])
+# filter2            = openmc.CellFilter([13])
+# tally              = openmc.Tally(name='tally kappa-fission')
+# tally.filters      = [filter1, filter2]
+# tally.scores       = ['kappa-fission']
+#
+# tally              = openmc.Tally(name='tally kappa-fission')
+# tally.filters      = [openmc.CellFilter([c_pebble_1, c_pebble_2])]
+# tally.scores       = ['kappa-fission']
+# tallies            = openmc.Tallies([tally])
+# tallies.export_to_xml()
+# -------------- Plots --------------
+# Plot parameters
+plot_width   = max(2*radius_vessel,height_vessel)
+plot_zcenter = z1_vessel+height_vessel/2
+m_colors = {m_fuel:              'brown',
+            m_graphite_c_buffer: 'LightSteelBlue',
+            m_graphite_pyc:      'blue',
+            m_sic:               'orange',
+            m_graphite_matrix:   'cyan',
+            m_graphite_inner:    'DeepSkyBlue',
+            m_graphite_outer:    'Navy',
+            m_flibe:             'yellow'}
+#
+plot1          = openmc.Plot()
+plot1.filename = 'plot1'
+plot1.width    = (2*radius_vessel, 2*radius_vessel)
+plot1.basis    = 'xy'
+plot1.origin   = (0, 0, 0)
+plot1.pixels   = (1000, 1000)
+plot1.color_by = 'material'
+plot1.colors   = m_colors
+#
+plot2          = openmc.Plot()
+plot2.filename = 'plot2'
+plot2.width    = (2*radius_vessel, 2*radius_vessel)
+plot2.basis    = 'xy'
+plot2.origin   = (0, 0, 1)
+plot2.pixels   = (1000, 1000)
+plot2.color_by = 'material'
+plot2.colors   = m_colors
+#
+plot3          = openmc.Plot()
+plot3.filename = 'plot3'
+plot3.width    = (2*radius_vessel, 2*radius_vessel)
+plot3.basis    = 'xy'
+plot3.origin   = (0, 0, 1)
+plot3.pixels   = (1000, 1000)
+plot3.color_by = 'material'
+plot3.colors   = m_colors
+#
+plot4          = openmc.Plot()
+plot4.filename = 'plot4'
+plot4.width    = (plot_width, plot_width)
+plot4.basis    = 'xz'
+plot4.origin   = (0, 0 ,plot_zcenter)
+plot4.pixels   = (1000, 1000)
+plot4.color_by = 'material'
+plot4.colors   = m_colors
+#
+plots          = openmc.Plots([plot1, plot2, plot3, plot4])
+if voxel==1:
+   plotv1 = openmc.Plot()
+   plotv1.filename = 'plotv1'
+   plotv1.width    = (2*radius_pebble_outer, 2*radius_pebble_outer,  2*radius_pebble_outer)
+   plotv1.origin   = (0, 0 ,0)
+   plotv1.pixels   = (300, 300, 300)
+   plotv1.color_by = 'material'
+   plotv1.colors   = m_colors
+   plotv1.type     = 'voxel'
+   plots           = openmc.Plots([plotv1])
+plots.export_to_xml()
--- a/1568_pebble_yuhsiang_E524k/geometry.xml
+++ b/1568_pebble_yuhsiang_E524k/geometry.xml
--- a/1568_pebble_yuhsiang_E524k/materials.xml
+++ b/1568_pebble_yuhsiang_E524k/materials.xml
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material depletable="true" id="1" name="m_fuel - uranium oxycarbide - triso partciles">
+    <density units="g/cm3" value="10.5" />
+    <nuclide name="U234" wo="0.03565284" />
+    <nuclide name="U235" wo="17.790767159999998" />
+    <nuclide name="U238" wo="71.75358" />
+    <nuclide name="C0" wo="1.8" />
+    <nuclide name="O16" wo="8.616527990092724" />
+    <nuclide name="O17" wo="0.003472009907274255" />
+  </material>
+  <material id="2" name="m_graphite_c_buffer - triso partciles">
+    <density units="g/cm3" value="1.0" />
+    <nuclide ao="1.0" name="C0" />
+    <sab name="c_Graphite" />
+  </material>
+  <material id="3" name="m_graphite_pyc - triso particles">
+    <density units="g/cm3" value="1.87" />
+    <nuclide ao="1.0" name="C0" />
+    <sab name="c_Graphite" />
+  </material>
+  <material id="4" name="sic - triso partciles">
+    <density units="g/cm3" value="3.2" />
+    <nuclide ao="1.0" name="C0" />
+    <nuclide ao="0.9222968" name="Si28" />
+    <nuclide ao="0.0468316" name="Si29" />
+    <nuclide ao="0.0308716" name="Si30" />
+  </material>
+  <material id="5" name="m_graphite_matrix - triso particles">
+    <density units="g/cm3" value="1.6" />
+    <nuclide ao="1.0" name="C0" />
+    <sab name="c_Graphite" />
+  </material>
+  <material id="6" name="m_graphite_inner - pebble inner zone">
+    <density units="g/cm3" value="1.59368" />
+    <nuclide ao="1.0" name="C0" />
+    <sab name="c_Graphite" />
+  </material>
+  <material id="7" name="m_graphite_outer - pebble outer zone">
+    <density units="g/cm3" value="1.74" />
+    <nuclide ao="1.0" name="C0" />
+    <sab name="c_Graphite" />
+  </material>
+  <material id="8" name="m_flibe - 2LiF-BeF2">
+    <density units="kg/m3" value="1962.13354" />
+    <nuclide ao="1.9999" name="Li7" />
+    <nuclide ao="9.999999999998899e-05" name="Li6" />
+    <nuclide ao="1.0" name="Be9" />
+    <nuclide ao="4.0" name="F19" />
+  </material>
+</materials>
--- a/1568_pebble_yuhsiang_E524k/nek.i
+++ b/1568_pebble_yuhsiang_E524k/nek.i
+[Mesh]
+  type = NekMesh
+[]
+
+[Variables]
+  [dummy]
+  []
+[]
+
+[AuxVariables]
+  [avg_flux]
+  []
+[]
+
+[Problem]
+  type = NekProblem
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 50000
+[]
+
+[Outputs]
+  [exo]
+    type = Exodus
+    output_dimension = 3
+    interval = 100
+    execute_on = 'timestep_end'
+  []
+[]
+
+[Postprocessors]
+  [total_flux]
+    type = Receiver
+    default = 0
+  []
+[]
--- a/1568_pebble_yuhsiang_E524k/nek_master.i
+++ b/1568_pebble_yuhsiang_E524k/nek_master.i
+[Mesh]
+  type = FileMesh
+  file = spheres_generator_in.e
+[]
+
+[Kernels]