4.3. Reactor pin cell
[ ]:
import os
import sys
import numpy as np
4.3.1. Using this Notebook
Before running this example, make sure that the Python module of OpenSn was installed.
4.3.1.1. Converting and Running this Notebook from the Terminal
To run this notebook from the terminal, simply type:
jupyter nbconvert --to python --execute <notebook_name>.ipynb.
To run this notebook in parallel (for example, using 4 processes), simply type:
mpiexec -n 4 jupyter nbconvert --to python --execute <notebook_name>.ipynb.
[ ]:
from mpi4py import MPI
size = MPI.COMM_WORLD.size
rank = MPI.COMM_WORLD.rank
if rank == 0:
print(f"Running with {size} MPI processors.")
4.3.2. Import Requirements
Import required classes and functions from the Python interface of OpenSn. Make sure that the path to PyOpenSn is appended to Python’s PATH.
[ ]:
# assuming that the execute dir is the notebook dir
# this line is not necessary when PyOpenSn is installed using pip
sys.path.append("../../..")
from pyopensn.mesh import FromFileMeshGenerator
from pyopensn.xs import MultiGroupXS
from pyopensn.aquad import GLCProductQuadrature2DXY
from pyopensn.solver import DiscreteOrdinatesProblem, PowerIterationKEigenSolver
from pyopensn.fieldfunc import FieldFunctionGridBased, FieldFunctionInterpolationVolume
from pyopensn.logvol import RPPLogicalVolume
from pyopensn.context import UseColor, Finalize
[ ]:
UseColor(False)
4.3.3. Mesh
We load a 3x3 reactor pin lattice (.obj file). The mesh and its block IDs are shown below:
div>
[ ]:
meshgen = FromFileMeshGenerator(filename="pincell.obj")
grid = meshgen.Execute()
grid.ExportToPVTU("pincell_mesh")
4.3.4. Cross Sections
We load 361-group cross sections that were generated using OpenMC.
[ ]:
xs_filepath = "./mgxs_2B_one_eighth_SHEM-361.h5"
h5_mat_names = [
"fuel",
"fuel clad",
"fuel gap",
"moderator",
]
xs_dict = {}
xs_list = []
for name in h5_mat_names:
xs_dict[name] = MultiGroupXS()
xs_dict[name].LoadFromOpenMC(xs_filepath, name, 294.0)
xs_list = np.append(xs_list, xs_dict[name])
block_ids = [i for i in range(0, len(xs_list))]
scat_order = 3 # xs_list[0].scattering_order
num_groups = xs_list[0].num_groups
4.3.5. Angular Quadrature
We create a product Gauss-Legendre-Chebyshev angular quadrature and pass the total number of polar cosines (here n_polar = 2) and the number of azimuthal subdivisions in four quadrants (n_azimuthal = 4).
For more accurate results, we suggest using n_polar = 8 and n_azimuthal = 32
[ ]:
pquad = GLCProductQuadrature2DXY(n_polar=2, n_azimuthal=4)
4.3.6. Linear Boltzmann Solver
4.3.6.1. Options for the Linear Boltzmann Problem (LBS)
[ ]:
group_sets = [
{
"groups_from_to": (0, num_groups - 1),
"angular_quadrature": pquad,
"angle_aggregation_type": "polar",
"inner_linear_method": "classic_richardson",
"l_abs_tol": 1.0e-5,
"l_max_its": 300,
}
]
bound_conditions = [
{"name": "xmin", "type": "reflecting"},
{"name": "xmax", "type": "reflecting"},
{"name": "ymin", "type": "reflecting"},
{"name": "ymax", "type": "reflecting"},
{"name": "zmin", "type": "reflecting"},
{"name": "zmax", "type": "reflecting"},
]
xs_mapping = [
{"block_ids": [0], "xs": xs_list[0]},
{"block_ids": [1], "xs": xs_list[1]},
{"block_ids": [2], "xs": xs_list[2]},
{"block_ids": [3], "xs": xs_list[3]}
]
phys = DiscreteOrdinatesProblem(mesh=grid, num_groups=num_groups, groupsets=group_sets, xs_map=xs_mapping)
phys.SetOptions(
scattering_order=scat_order,
verbose_inner_iterations=True,
verbose_outer_iterations=True,
use_precursors=False,
power_default_kappa=1.0,
power_normalization=1.0,
save_angular_flux=False,
boundary_conditions=bound_conditions,
)
4.3.6.2. Putting the Linear Boltzmann Solver Together
We then create the physics solver, initialize it, and execute it.
[ ]:
k_solver = PowerIterationKEigenSolver(lbs_problem=phys, k_tol=1.0e-14)
k_solver.Initialize()
k_solver.Execute()
keff = k_solver.GetEigenvalue()
if rank ==0:
print(f"Eigenvalue = {keff}")
4.3.7. Post-Processing via Field Functions
[ ]:
fflist = phys.GetScalarFieldFunctionList()
vtk_basename = "pin_cell"
FieldFunctionGridBased.ExportMultipleToVTK([fflist[g] for g in range(num_groups)], vtk_basename)
4.3.8. Post-processing: Extract the average flux in a portion of the domain
We create an edit zone (logical volume) that is the entire domain.
We request the average (keyword "avg") of the scalar flux over the edit zone, for each group.
[ ]:
logvol_whole_domain = RPPLogicalVolume(infx=True, infy=True, infz=True)
[ ]:
flux = np.zeros(num_groups)
for g in range(0, num_groups):
ffi = FieldFunctionInterpolationVolume()
ffi.SetOperationType("sum")
ffi.SetLogicalVolume(logvol_whole_domain)
ffi.AddFieldFunction(fflist[g])
ffi.Initialize()
ffi.Execute()
flux[g] = ffi.GetValue()
flux /= np.sum(flux)
[ ]:
# 361-group structure (copied from openmc/openmc/mgxs/__init__.py)
group_edges = np.array([
0.00000e+00, 2.49990e-03, 4.55602e-03, 7.14526e-03, 1.04505e-02,
1.48300e-02, 2.00104e-02, 2.49394e-02, 2.92989e-02, 3.43998e-02,
4.02999e-02, 4.73019e-02, 5.54982e-02, 6.51999e-02, 7.64969e-02,
8.97968e-02, 1.04298e-01, 1.19995e-01, 1.37999e-01, 1.61895e-01,
1.90005e-01, 2.09610e-01, 2.31192e-01, 2.54997e-01, 2.79989e-01,
3.05012e-01, 3.25008e-01, 3.52994e-01, 3.90001e-01, 4.31579e-01,
4.75017e-01, 5.20011e-01, 5.54990e-01, 5.94993e-01, 6.24999e-01,
7.19999e-01, 8.00371e-01, 8.80024e-01, 9.19978e-01, 9.44022e-01,
9.63960e-01, 9.81959e-01, 9.96501e-01, 1.00904e+00, 1.02101e+00,
1.03499e+00, 1.07799e+00, 1.09198e+00, 1.10395e+00, 1.11605e+00,
1.12997e+00, 1.14797e+00, 1.16999e+00, 1.21397e+00, 1.25094e+00,
1.29304e+00, 1.33095e+00, 1.38098e+00, 1.41001e+00, 1.44397e+00,
1.51998e+00, 1.58803e+00, 1.66895e+00, 1.77997e+00, 1.90008e+00,
1.98992e+00, 2.07010e+00, 2.15695e+00, 2.21709e+00, 2.27299e+00,
2.33006e+00, 2.46994e+00, 2.55000e+00, 2.59009e+00, 2.62005e+00,
2.64004e+00, 2.70012e+00, 2.71990e+00, 2.74092e+00, 2.77512e+00,
2.88405e+00, 3.14211e+00, 3.54307e+00, 3.71209e+00, 3.88217e+00,
4.00000e+00, 4.21983e+00, 4.30981e+00, 4.41980e+00, 4.76785e+00,
4.93323e+00, 5.10997e+00, 5.21008e+00, 5.32011e+00, 5.38003e+00,
5.41025e+00, 5.48817e+00, 5.53004e+00, 5.61979e+00, 5.72015e+00,
5.80021e+00, 5.96014e+00, 6.05991e+00, 6.16011e+00, 6.28016e+00,
6.35978e+00, 6.43206e+00, 6.48178e+00, 6.51492e+00, 6.53907e+00,
6.55609e+00, 6.57184e+00, 6.58829e+00, 6.60611e+00, 6.63126e+00,
6.71668e+00, 6.74225e+00, 6.75981e+00, 6.77605e+00, 6.79165e+00,
6.81070e+00, 6.83526e+00, 6.87021e+00, 6.91778e+00, 6.99429e+00,
7.13987e+00, 7.38015e+00, 7.60035e+00, 7.73994e+00, 7.83965e+00,
7.97008e+00, 8.13027e+00, 8.30032e+00, 8.52407e+00, 8.67369e+00,
8.80038e+00, 8.97995e+00, 9.14031e+00, 9.50002e+00, 1.05793e+01,
1.08038e+01, 1.10529e+01, 1.12694e+01, 1.15894e+01, 1.17094e+01,
1.18153e+01, 1.19795e+01, 1.21302e+01, 1.23086e+01, 1.24721e+01,
1.26000e+01, 1.33297e+01, 1.35460e+01, 1.40496e+01, 1.42505e+01,
1.44702e+01, 1.45952e+01, 1.47301e+01, 1.48662e+01, 1.57792e+01,
1.60498e+01, 1.65501e+01, 1.68305e+01, 1.74457e+01, 1.75648e+01,
1.77590e+01, 1.79591e+01, 1.90848e+01, 1.91997e+01, 1.93927e+01,
1.95974e+01, 2.00734e+01, 2.02751e+01, 2.04175e+01, 2.05199e+01,
2.06021e+01, 2.06847e+01, 2.07676e+01, 2.09763e+01, 2.10604e+01,
2.11448e+01, 2.12296e+01, 2.13360e+01, 2.14859e+01, 2.17018e+01,
2.20011e+01, 2.21557e+01, 2.23788e+01, 2.25356e+01, 2.46578e+01,
2.78852e+01, 3.16930e+01, 3.30855e+01, 3.45392e+01, 3.56980e+01,
3.60568e+01, 3.64191e+01, 3.68588e+01, 3.73038e+01, 3.77919e+01,
3.87874e+01, 3.97295e+01, 4.12270e+01, 4.21441e+01, 4.31246e+01,
4.41721e+01, 4.52904e+01, 4.62053e+01, 4.75173e+01, 4.92591e+01,
5.17847e+01, 5.29895e+01, 5.40600e+01, 5.70595e+01, 5.99250e+01,
6.23083e+01, 6.36306e+01, 6.45923e+01, 6.50460e+01, 6.55029e+01,
6.58312e+01, 6.61612e+01, 6.64929e+01, 6.68261e+01, 6.90682e+01,
7.18869e+01, 7.35595e+01, 7.63322e+01, 7.93679e+01, 8.39393e+01,
8.87741e+01, 9.33256e+01, 9.73287e+01, 1.00594e+02, 1.01098e+02,
1.01605e+02, 1.02115e+02, 1.03038e+02, 1.05646e+02, 1.10288e+02,
1.12854e+02, 1.15480e+02, 1.16524e+02, 1.17577e+02, 1.20554e+02,
1.26229e+02, 1.32701e+02, 1.39504e+02, 1.46657e+02, 1.54176e+02,
1.63056e+02, 1.67519e+02, 1.75229e+02, 1.83295e+02, 1.84952e+02,
1.86251e+02, 1.87559e+02, 1.88877e+02, 1.90204e+02, 1.93078e+02,
1.95996e+02, 2.00958e+02, 2.12108e+02, 2.24325e+02, 2.35590e+02,
2.41796e+02, 2.56748e+02, 2.68297e+02, 2.76468e+02, 2.84888e+02,
2.88327e+02, 2.95922e+02, 3.19928e+02, 3.35323e+02, 3.53575e+02,
3.71703e+02, 3.90760e+02, 4.19094e+02, 4.53999e+02, 5.01746e+02,
5.39204e+02, 5.77146e+02, 5.92941e+02, 6.00099e+02, 6.12834e+02,
6.46837e+02, 6.77287e+02, 7.48517e+02, 8.32218e+02, 9.09681e+02,
9.82494e+02, 1.06432e+03, 1.13467e+03, 1.34358e+03, 1.58620e+03,
1.81183e+03, 2.08410e+03, 2.39729e+03, 2.70024e+03, 2.99618e+03,
3.48107e+03, 4.09735e+03, 5.00451e+03, 6.11252e+03, 7.46585e+03,
9.11881e+03, 1.11377e+04, 1.36037e+04, 1.48997e+04, 1.62005e+04,
1.85847e+04, 2.26994e+04, 2.49991e+04, 2.61001e+04, 2.73944e+04,
2.92810e+04, 3.34596e+04, 3.69786e+04, 4.08677e+04, 4.99159e+04,
5.51656e+04, 6.73794e+04, 8.22974e+04, 9.46645e+04, 1.15624e+05,
1.22773e+05, 1.40000e+05, 1.64999e+05, 1.95008e+05, 2.30014e+05,
2.67826e+05, 3.20646e+05, 3.83884e+05, 4.12501e+05, 4.56021e+05,
4.94002e+05, 5.78443e+05, 7.06511e+05, 8.60006e+05, 9.51119e+05,
1.05115e+06, 1.16205e+06, 1.28696e+06, 1.33694e+06, 1.40577e+06,
1.63654e+06, 1.90139e+06, 2.23130e+06, 2.72531e+06, 3.32871e+06,
4.06569e+06, 4.96585e+06, 6.06530e+06, 6.70319e+06, 7.40817e+06,
8.18730e+06, 9.04836e+06, 9.99999e+06, 1.16183e+07, 1.38403e+07,
1.49182e+07, 1.96403e+07])
# flip group edges to have highest energies first, since fastest groups have the lowest index
E = np.flip(group_edges)
# compute the group widths
dE = -np.diff(E)
# compute the group midpoints
Emid = E[:-1] + dE/2
[ ]:
import matplotlib.pyplot as plt
if rank ==0:
plt.figure()
y = Emid * flux / dE
y = np.insert(y, 0, y[0])
plt.semilogx(E, y, drawstyle='steps',label='flux')
plt.title('Lethargy Flux')
plt.legend()
plt.grid()
# plt.savefig("./images/pincell_lethargy_spectrum.png")
# plt.show()
plt.figure()
y = flux / dE
y = np.insert(y, 0, y[0])
plt.loglog(E, y, drawstyle='steps',label='flux')
plt.title('Flux')
plt.legend()
plt.grid()
# plt.savefig("./images/pincell_spectrum.png")
# plt.show()
The resulting spectra are shown below:
4.3.9. Finalize (for Jupyter Notebook only)
In Python script mode, PyOpenSn automatically handles environment termination. However, this automatic finalization does not occur when running in a Jupyter notebook, so explicit finalization of the environment at the end of the notebook is required. Do not call the finalization in Python script mode, or in console mode.
Note that PyOpenSn’s finalization must be called before MPI’s finalization.
[ ]:
MPI.COMM_WORLD.Barrier()
[ ]:
from IPython import get_ipython
def finalize_env():
Finalize()
MPI.Finalize()
ipython_instance = get_ipython()
if ipython_instance is not None:
ipython_instance.events.register("post_execute", finalize_env)