5.4.2. Switching From Steady-State to Transient to Steady-State in OpenSn
5.4.2.1. Using this Notebook
Before running this example, make sure the Python module of OpenSn is available in your environment.
5.4.2.1.1. Converting and Running this Notebook from the Terminal
To run this notebook from the terminal, use:
jupyter nbconvert --to python --execute source_step_steady_to_transient_to_steady.ipynb
To run in parallel (for example, 4 ranks):
mpiexec -n 4 jupyter nbconvert --to python --execute source_step_steady_to_transient_to_steady.ipynb
5.4.2.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.2.2.1. Disable colorized output
For notebook rendering, plain (non-colorized) output is typically easier to read.
[1]:
import math
import os
import sys
import ctypes
from pathlib import Path
if "opensn_console" not in globals():
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.source import VolumetricSource
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.solver import (
DiscreteOrdinatesProblem,
SteadyStateSourceSolver,
TransientSolver,
)
from pyopensn.fieldfunc import FieldFunctionInterpolationVolume
from pyopensn.logvol import RPPLogicalVolume
else:
# In opensn_console mode, solver symbols are preloaded by the runtime.
rank = 0
def max_phi(problem):
fflist = problem.GetScalarFluxFieldFunction(only_scalar_flux=True)
monitor_volume = RPPLogicalVolume(infx=True, infy=True, infz=True)
field_interp = FieldFunctionInterpolationVolume()
field_interp.SetOperationType("max")
field_interp.SetLogicalVolume(monitor_volume)
field_interp.AddFieldFunction(fflist[0])
field_interp.Initialize()
field_interp.Execute()
return field_interp.GetValue()
5.4.2.3. Mesh
We construct a 1D Cartesian mesh (10 cells on ([0,1])) so this tutorial focuses on workflow behavior rather than geometry complexity.
[2]:
num_cells_1d = 10
L = 1.0
dx = L / num_cells_1d
nodes = [i * dx for i in range(num_cells_1d + 1)]
meshgen = OrthogonalMeshGenerator(node_sets=[nodes])
grid = meshgen.Execute()
OpenSn version 0.0.1
2026-03-03 09:59:46 Running OpenSn with 1 processes.
[0] Done checking cell-center-to-face orientations
[0] 00:00:00.0 Establishing cell connectivity.
[0] 00:00:00.0 Vertex cell subscriptions complete.
[0] 00:00:00.0 Surpassing cell 1 of 10 (10%)
[0] 00:00:00.0 Surpassing cell 2 of 10 (20%)
[0] 00:00:00.0 Surpassing cell 4 of 10 (30%)
[0] 00:00:00.0 Surpassing cell 5 of 10 (40%)
[0] 00:00:00.0 Surpassing cell 6 of 10 (50%)
[0] 00:00:00.0 Surpassing cell 7 of 10 (60%)
[0] 00:00:00.0 Surpassing cell 8 of 10 (70%)
[0] 00:00:00.0 Surpassing cell 9 of 10 (80%)
[0] 00:00:00.0 Surpassing cell 10 of 10 (90%)
[0] 00:00:00.0 Establishing cell boundary connectivity.
[0] 00:00:00.0 Done establishing cell connectivity.
[0] Number of cells per partition (max,min,avg) = 10,10,10
[0]
[0] Mesh statistics:
[0] Global cell count : 10
[0] Local cell count (avg,max,min): 10,10,10
[0] Ghost-to-local ratio (avg) : 0
[0]
5.4.2.4. Material IDs
Assign a single material ID (0) to all cells.
[3]:
grid.SetUniformBlockID(0)
[0] 00:00:00.0 Done setting block id 0 to all cells
5.4.2.5. Cross Sections
Use a one-group pure absorber model:
\[\Sigma_t = 1.0\]
\[\Sigma_s = 0.0\]
For this setup, the steady source-driven solution is
\[\phi_{ss} = Q / \Sigma_t\]
.
[4]:
sigma_t = 1.0
sigma_s = 0.0
q0 = 2.0 # initial steady-state source strength
q1 = 3.0 # post-step source strength
dt = 0.1
t_end = 0.5
xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t, sigma_s, 1.0)
src0 = VolumetricSource(block_ids=[0], group_strength=[q0], start_time=0.0, end_time=1.0e9)
src1 = VolumetricSource(block_ids=[0], group_strength=[q1], start_time=0.0, end_time=1.0e9)
5.4.2.6. Angular Quadrature
Use a low-order 1D slab quadrature.
[5]:
pquad = GLProductQuadrature1DSlab(n_polar=4, scattering_order=0)
5.4.2.7. Linear Boltzmann Solver
Create one DiscreteOrdinatesProblem object that is reused across all three phases:
steady-state source solve with
Q0,transient with source stepped to
Q1,final steady-state source solve at
Q1.
[6]:
phys = DiscreteOrdinatesProblem(
mesh=grid,
num_groups=1,
groupsets=[
{
"groups_from_to": (0, 0),
"angular_quadrature": pquad,
"inner_linear_method": "petsc_richardson",
"l_abs_tol": 1.0e-8,
"l_max_its": 200,
}
],
xs_map=[{"block_ids": [0], "xs": xs}],
volumetric_sources=[src0],
boundary_conditions=[
{"name": "zmin", "type": "reflecting"},
{"name": "zmax", "type": "reflecting"},
],
options={
"save_angular_flux": True,
"verbose_inner_iterations": False,
"verbose_outer_iterations": False,
},
)
[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: 4
[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] Volumetric source #0 has 10 total subscribing cells.
[0] 00:00:00.0 Initializing sweep datastructures.
[0] 00:00:00.0 Done initializing sweep datastructures.
5.4.2.8. Steady-State Solve (Initial Source (Q_0))
Solve a source-driven steady-state problem and compare against
\[\phi_{ss,0} = Q_0 / \Sigma_t\]
.
[7]:
steady0 = SteadyStateSourceSolver(problem=phys)
steady0.Initialize()
steady0.Execute()
phi_ss0_num = max_phi(phys)
phi_ss0_ana = q0 / sigma_t
if rank == 0:
print(f"Initial steady numeric phi_ss0 = {phi_ss0_num:.8e}")
print(f"Initial steady analytic phi_ss0 = {phi_ss0_ana:.8e}")
Initial steady numeric phi_ss0 = 2.00000000e+00
Initial steady analytic phi_ss0 = 2.00000000e+00
5.4.2.9. Transition to Transient Mode and Step Source to (Q_1)
Switch to time-dependent mode, initialize from the existing steady state, and replace the volumetric source with the stepped value Q1.
[12]:
phys.SetTimeDependentMode()
transient_solver = TransientSolver(problem=phys, dt=dt, theta=1.0, initial_state="existing")
transient_solver.Initialize()
phys.SetVolumetricSources(clear_volumetric_sources=True)
phys.SetVolumetricSources(volumetric_sources=[src1])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 2
1 phys.SetTimeDependentMode()
----> 2 transient_solver = TransientSolver(problem=phys, dt=dt, theta=1.0, verbose=false, initial_state="existing")
3 transient_solver.Initialize()
5 phys.SetVolumetricSources(clear_volumetric_sources=True)
NameError: name 'false' is not defined
5.4.2.10. Advance in Time (Transient Solve)
Advance from (t=0) to (t=t_{end}) using backward Euler.
[9]:
step = 0
while phys.GetTime() < t_end:
t_before = phys.GetTime()
transient_solver.Advance()
t_after = phys.GetTime()
phi_now = max_phi(phys)
if rank == 0:
print(
f"step={step:03d}, t={t_before:.4f}->{t_after:.4f}, phi_max={phi_now:.8e}"
)
step += 1
phi_td_num = max_phi(phys)
n_steps = step
step=000, t=0.0000->0.1000, phi_max=2.09090909e+00[0] TransientSolver dt = 1.0e-01 time = 0.1000
step=001, t=0.1000->0.2000, phi_max=2.17355372e+00
step=002, t=0.2000->0.3000, phi_max=2.24868520e+00
step=003, t=0.3000->0.4000, phi_max=2.31698654e+00
step=004, t=0.4000->0.5000, phi_max=2.37907868e+00
[0] TransientSolver dt = 1.0e-01 time = 0.2000
[0] TransientSolver dt = 1.0e-01 time = 0.3000
[0] TransientSolver dt = 1.0e-01 time = 0.4000
[0] TransientSolver dt = 1.0e-01 time = 0.5000
5.4.2.11. Return to Steady-State Mode (Final Source (Q_1))
Switch back to steady-state mode and solve again with the same final source Q1. The final steady solution should satisfy
\[\phi_{ss,1} = Q_1/\Sigma_t\]
.
[10]:
phys.SetSteadyStateMode()
steady1 = SteadyStateSourceSolver(problem=phys)
steady1.Initialize()
steady1.Execute()
phi_ss1_num = max_phi(phys)
phi_ss1_ana = q1 / sigma_t
5.4.2.12. Post-processing: Analytic Checks for the Full Workflow
For backward Euler with constant source (Q_1), the transient update is
\[\phi^{n+1} = \frac{\phi^n + \Delta t\,Q_1}{1 + \Sigma_t\,\Delta t}\]
.
With initial value
\[\phi^0 = Q_0/\Sigma_t\]
, after (N) steps:
\[\phi^N = \phi_{ss,1} + \left(\phi^0-\phi_{ss,1}ight)
\left(\frac{1}{1+\Sigma_t\Delta t}ight)^N,
\quad \phi_{ss,1}=\frac{Q_1}{\Sigma_t}.\]
We compare:
transient numerical (\phi`^N) vs analytic (:nbsphinx-math:phi`^N),
final steady numerical (\phi_{ss,1}) vs analytic (Q_1/\Sigma_t).
[11]:
r = 1.0 / (1.0 + sigma_t * dt)
phi_td_ana = phi_ss1_ana + (phi_ss0_ana - phi_ss1_ana) * (r ** n_steps)
err_td = abs(phi_td_num - phi_td_ana) / max(abs(phi_td_ana), 1.0e-16)
err_ss = abs(phi_ss1_num - phi_ss1_ana) / max(abs(phi_ss1_ana), 1.0e-16)
tol_td = 5.0e-3
tol_ss = 1.0e-3
ok = (err_td < tol_td) and (err_ss < tol_ss)
if rank == 0:
print(f"Transient-end numeric phi = {phi_td_num:.8e}")
print(f"Transient-end analytic phi = {phi_td_ana:.8e}")
print(f"Final steady numeric phi = {phi_ss1_num:.8e}")
print(f"Final steady analytic phi = {phi_ss1_ana:.8e}")
print(f"Relative error (transient) = {err_td:.3e}")
print(f"Relative error (steady) = {err_ss:.3e}")
print(f"SOURCE_STEP_WORKFLOW_OK {1 if ok else 0}")
assert ok, (
f"Workflow check failed: err_td={err_td:.3e}, err_ss={err_ss:.3e}, "
f"tol_td={tol_td:.3e}, tol_ss={tol_ss:.3e}"
)
Transient-end numeric phi = 2.37907868e+00
Transient-end analytic phi = 2.37907868e+00
Final steady numeric phi = 3.00000000e+00
Final steady analytic phi = 3.00000000e+00
Relative error (transient) = 2.326e-11
Relative error (steady) = 9.089e-11
SOURCE_STEP_WORKFLOW_OK 1
5.4.2.13. Finalize (for Jupyter Notebook only)
In Python script mode, PyOpenSn usually handles environment termination automatically. In Jupyter notebook mode, finalize explicitly only if your runtime requires it.