5.4.1. Switching From Steady-State Mode to Transient Mode in OpenSn
5.4.1.1. Using this Notebook
Before running this example, make sure the Python module of OpenSn is available in your environment.
5.4.1.1.1. Converting and Running this Notebook from the Terminal
To run this notebook from the terminal, use:
jupyter nbconvert --to python --execute keigen_to_transient_tutorial.ipynb
To run in parallel (for example, 4 ranks):
mpiexec -n 4 jupyter nbconvert --to python --execute keigen_to_transient_tutorial.ipynb
5.4.1.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.
5.4.1.2.1. Disable colorized output
For notebook rendering, plain (non-colorized) output is typically easier to read.
[38]:
import math
import os
import sys
import ctypes
from pathlib import Path
from mpi4py import MPI
size = MPI.COMM_WORLD.size
rank = MPI.COMM_WORLD.rank
# In notebooks, __file__ is typically undefined. Prefer cwd as the tutorial dir.
tutorial_dir = Path(__file__).resolve().parent if "__file__" in globals() else Path.cwd()
# If launched from repo root (or another directory), recover the tutorial directory.
if not (tutorial_dir / "xs1g_prompt_crit.cxs").exists():
candidate = Path.cwd() / "doc/source/tutorials/lbs/transient"
if candidate.exists():
tutorial_dir = candidate
repo_root = None
build_dir = None
# Prefer normal environment import (e.g., site-packages), then fall back to local build.
try:
import pyopensn
except ImportError:
# Find repo root robustly by walking upward.
for p in [tutorial_dir, *tutorial_dir.parents]:
has_src_pyopensn = (p / "pyopensn").exists() and (p / "CMakeLists.txt").exists()
has_bld_pyopensn = (p / "build" / "pyopensn").exists()
if has_bld_pyopensn or has_src_pyopensn:
repo_root = p
break
if repo_root is None:
repo_root = tutorial_dir
build_dir = repo_root / "build"
if (build_dir / "pyopensn").exists():
# Avoid stale module path order issues.
sys.path = [p for p in sys.path if p not in (str(build_dir), str(repo_root))]
sys.path.insert(0, str(build_dir))
# Preload local libopensn if available (Linux/macOS).
lib_candidates = (
build_dir / "libopensn.so",
build_dir / "libopensn.dylib",
build_dir / "lib" / "libopensn.so",
build_dir / "lib" / "libopensn.dylib",
)
for lib_path in lib_candidates:
if lib_path.exists():
ctypes.CDLL(str(lib_path), mode=ctypes.RTLD_GLOBAL)
break
import pyopensn
# If we used fallback import, ensure we loaded from local build.
pyopensn_path = Path(pyopensn.__file__).resolve()
expected_parent = (build_dir / "pyopensn").resolve()
if expected_parent.exists() and expected_parent not in pyopensn_path.parents:
raise RuntimeError(
f"Imported pyopensn from unexpected location: {pyopensn_path}. "
f"Expected under {expected_parent}"
)
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.xs import MultiGroupXS
from pyopensn.aquad import GLCProductQuadrature3DXYZ
from pyopensn.solver import (DiscreteOrdinatesProblem,
PowerIterationKEigenSolver,
TransientSolver,
)
5.4.1.3. Mesh
We construct a simple 3D Cartesian mesh with 64 cells using the built-in orthogonal mesh generator. The dimensions of the mesh are 8cm in each dimension.
[39]:
num_cells_1d = 4
L = 8.0
dx = L / num_cells_1d
nodes = [i * dx for i in range(num_cells_1d + 1)]
meshgen = OrthogonalMeshGenerator(node_sets=[nodes, nodes, nodes])
grid = meshgen.Execute()
[0] Done checking cell-center-to-face orientations
[0] 07:25:15.7 Establishing cell connectivity.
[0] 07:25:15.7 Vertex cell subscriptions complete.
[0] 07:25:15.7 Surpassing cell 7 of 64 (10%)
[0] 07:25:15.7 Surpassing cell 13 of 64 (20%)
[0] 07:25:15.7 Surpassing cell 20 of 64 (30%)
[0] 07:25:15.7 Surpassing cell 26 of 64 (40%)
[0] 07:25:15.7 Surpassing cell 32 of 64 (50%)
[0] 07:25:15.7 Surpassing cell 39 of 64 (60%)
[0] 07:25:15.7 Surpassing cell 45 of 64 (70%)
[0] 07:25:15.7 Surpassing cell 52 of 64 (80%)
[0] 07:25:15.7 Surpassing cell 58 of 64 (90%)
[0] 07:25:15.7 Surpassing cell 64 of 64 (100%)
[0] 07:25:15.7 Establishing cell boundary connectivity.
[0] 07:25:15.7 Done establishing cell connectivity.
[0] Number of cells per partition (max,min,avg) = 64,64,64
[0]
[0] Mesh statistics:
[0] Global cell count : 64
[0] Local cell count (avg,max,min): 64,64,64
[0] Ghost-to-local ratio (avg) : 0
[0]
5.4.1.4. Material IDs
When using the built-in mesh generator, no material IDs are assigned. We assign a material ID with value 0 for each cell in the spatial domain.
[40]:
grid.SetUniformBlockID(0)
[0] 07:25:15.7 Done setting block id 0 to all cells
5.4.1.5. Cross Sections
Load two one-group cross-section sets:
xs_critfor the steady-state k-eigen solvexs_superfor the transient reactivity insertion
This gives a clean steady-state-to-transient transition with a controlled cross section swap.
[41]:
xs_crit = MultiGroupXS()
xs_crit.LoadFromOpenSn(str(tutorial_dir / "xs1g_prompt_crit.cxs"))
xs_super = MultiGroupXS()
xs_super.LoadFromOpenSn(str(tutorial_dir / "xs1g_prompt_super.cxs"))
[0] Reading OpenSn cross-section file "/Users/dhawkins/opensn_dev/opensn/doc/source/tutorials/lbs/transient/xs1g_prompt_crit.cxs"
[0] *** WARNING *** Estimating absorption from the transfer matrices.
[0] Reading OpenSn cross-section file "/Users/dhawkins/opensn_dev/opensn/doc/source/tutorials/lbs/transient/xs1g_prompt_super.cxs"
[0] *** WARNING *** Estimating absorption from the transfer matrices.
5.4.1.6. 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). This creates a 3D angular quadrature for XYZ geometry.
[42]:
pquad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0)
5.4.1.7. Linear Boltzmann Solver
Create a single DiscreteOrdinatesProblem with reflecting boundaries, a single groupset, reasonable iterative options, and an initial xs_map that assigns our xs_crit material to all cells. This problem object will be reused across both the steady-state and transient solves.
Two options are key here:
save_angular_flux=Truefor storing the angular flux state needed by the transient solve.use_precursors=Falsefor prompt-only transient behavior in this tutorial.
Note: In this workflow, set save_angular_flux=True on the shared problem object before the steady-state solve so the same problem can transition cleanly into transient mode.
[43]:
phys = DiscreteOrdinatesProblem(
mesh=grid,
num_groups=1,
groupsets=[
{
"groups_from_to": (0, 0),
"angular_quadrature": pquad,
"inner_linear_method": "classic_richardson",
"l_abs_tol": 1.0e-8,
"l_max_its": 200,
},
],
xs_map=[{"block_ids": [0], "xs": xs_crit}],
boundary_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"},
],
options={
"save_angular_flux": True,
"use_precursors": False,
"verbose_inner_iterations": False,
"verbose_outer_iterations": True,
},
)
[0]
[0] Initializing DiscreteOrdinatesProblem
[0]
[0] Scattering order : 0
[0] Number of moments : 1
[0] Number of groups : 1
[0] Number of groupsets : 1
[0]
[0] ***** Groupset 0 *****
[0] Number of angles: 8
[0] Groups:
[0] 0
[0]
[0] Initializing spatial discretization metadata.
[0] Computing unit integrals.
[0] Ghost cell unit cell-matrix ratio: 0%
[0] Cell matrices computed.
[0] Initializing parallel arrays. G=1 M=1
[0] Done with parallel arrays.
[0] 07:25:15.7 Initializing sweep datastructures.
[0] 07:25:15.7 Done initializing sweep datastructures.
5.4.1.8. Steady-State Solve
Run a power-iteration k-eigenvalue solve to convergence. This converged state becomes the initial condition for the transient solve.
[44]:
keigen = PowerIterationKEigenSolver(problem=phys, max_iters=200, k_tol=1.0e-10)
keigen.Initialize()
keigen.Execute()
[0] 07:25:15.8 Iteration 1 k_eff 1 k_eff change 2.436051e-11 reactivity -2.436051e-06 CONVERGED
[0]
[0] Final k-eigenvalue : 1
[0] Final change : 2.43605e-11 (Total number of sweeps:200)
[0]
[0] 00:00:00.0 Finished solver execution PowerIterationKEigenSolver.
[0]
5.4.1.9. Transition to Transient Mode
Now, let’s switch the problem object to time-dependent mode and initialize the transient solver. We’ll pass the initial_state="existing"option to the transient solver so that it uses the converged steady-state solution as its initial condition.
[45]:
phys.SetTimeDependentMode()
transient_solver = TransientSolver(problem=phys, verbose=False, initial_state="existing")
transient_solver.Initialize()
[0] 00:00:00.0 Initializing solver TransientSolver.
[0] *** WARNING *** TransientSolver: fissionable material is present but use_precursors is disabled. Running prompt-only transient.
5.4.1.10. Update Cross Sections and Advance in Time (Transient Solve)
Swap from xs_crit to xs_super, then march with backward Euler (theta=1.0) using a Python time stepping loop.
Diagnostic printouts show per-step time advancement and fission production.
[46]:
phys.SetXSMap(xs_map=[{"block_ids": [0], "xs": xs_super}])
# Baseline fission production at transient start (post-swap, t=0+)
t0 = phys.GetTime()
fp0 = phys.ComputeFissionProduction("new")
if rank == 0:
print(f"Transient start time t0 = {t0:.6f}")
print(f"Initial fission production fp0 (post-swap) = {fp0:.8e}")
transient_solver.SetTimeStep(1.0e-2)
transient_solver.SetTheta(1.0)
t_end = 0.1
step = 0
while phys.GetTime() < t_end:
t_before = phys.GetTime()
transient_solver.Advance()
t_after = phys.GetTime()
fp = phys.ComputeFissionProduction("new")
if rank == 0:
print(
f"Time step = {step:03d}, "
f"start time = {t_before:.6f}, "
f"end time = {t_after:.6f}, "
f"fission production @end = {fp:.8e}"
)
step += 1
if step == 0:
raise RuntimeError("No transient steps executed. Reinitialize solver state or increase t_end.")
fp_end = phys.ComputeFissionProduction("new")
t1 = phys.GetTime()
Transient start time t0 = 0.000000
Initial fission production fp0 (post-swap) = 1.84320000e+02
Time step = 000, start time = 0.000000, end time = 0.010000, fission production @end = 1.84422441e+02
Time step = 001, start time = 0.010000, end time = 0.020000, fission production @end = 1.84533161e+02
Time step = 002, start time = 0.020000, end time = 0.030000, fission production @end = 1.84643948e+02
Time step = 003, start time = 0.030000, end time = 0.040000, fission production @end = 1.84754801e+02
Time step = 004, start time = 0.040000, end time = 0.050000, fission production @end = 1.84865720e+02
Time step = 005, start time = 0.050000, end time = 0.060000, fission production @end = 1.84976706e+02
Time step = 006, start time = 0.060000, end time = 0.070000, fission production @end = 1.85087759e+02
Time step = 007, start time = 0.070000, end time = 0.080000, fission production @end = 1.85198878e+02
Time step = 008, start time = 0.080000, end time = 0.090000, fission production @end = 1.85310064e+02
Time step = 009, start time = 0.090000, end time = 0.100000, fission production @end = 1.85421317e+02
Time step = 010, start time = 0.100000, end time = 0.110000, fission production @end = 1.85532636e+02
5.4.1.11. Post-processing: Prompt-Only Growth Check
For this simple one-group setup we compare the fission-production fraction to:
In the notebook, we evaluate this over the transient interval
:
Code equivalent:
ratio_analytic = math.exp(alpha * delta_t)
The numerical fraction is computed directly from fission production:
Code equivalent:
ratio_numeric = fp_end / fp0 if fp0 > 0.0 else 0.0
Then we compare
.
[47]:
sigma_t = 1.0
sigma_s = 0.7
sigma_a = sigma_t - sigma_s
nu = 2.0
sigma_f = 0.180000
alpha = nu * sigma_f - sigma_a
delta_t = t1 - t0
ratio_analytic = math.exp(alpha * delta_t)
ratio_numeric = fp_end / fp0 if fp0 > 0.0 else 0.0
rel_err = abs(ratio_numeric - ratio_analytic) / max(abs(ratio_analytic), 1.0e-16)
tol = 5.0e-3 # matches prompt-analytic regression tolerance
if rank == 0:
print(f"Elapsed transient time = {delta_t:.6f} s")
print(f"FP ratio numeric = {ratio_numeric:.8e}")
print(f"FP ratio analytic = {ratio_analytic:.8e}")
print(f"Relative error = {rel_err:.3e}")
assert rel_err < tol, (
f"FP ratio check failed: rel_err={rel_err:.3e} >= tol={tol:.3e} "
f"(numeric={ratio_numeric:.8e}, analytic={ratio_analytic:.8e}, delta_t={delta_t:.6e})"
)
if rank == 0:
print("KEIGEN_TO_TRANSIENT_OK 1")
Elapsed transient time = 0.110000 s
FP ratio numeric = 1.00657897e+00
FP ratio analytic = 1.00662183e+00
Relative error = 4.257e-05
5.4.1.12. Finalize (for Jupyter Notebook only)
In script/console mode, PyOpenSn finalization is handled automatically. In notebook mode, finalize explicitly only if your runtime requires it.
5.4.1.12.1. Summary
Reused the same
DiscreteOrdinatesProblemobject for both steady-state and transient solves,solved steady-state k-eigen to convergence,
switched mode with
SetTimeDependentMode(),initialized transient with
initial_state="existing",applied cross section change and advanced from the converged steady state.