5.3.2. 3D Uncollided+Collided Calculation on an Unstructured Mesh
This notebook demonstrates a complete uncollided-plus-collided fixed-source calculation on a small tetrahedral mesh. The problem is intentionally simple:
single energy group,
single material,
single point source,
serial uncollided generation followed by a collided solve.
The mesh is a coarse tetrahedral cube (test/assets/mesh/uncollided_cube_coarse.msh).
5.3.2.1. Prerequisites
Run this notebook with the OpenSn Python module available. Because the uncollided stage must run with exactly one MPI rank, this end-to-end notebook should be executed in serial:
jupyter nbconvert --to notebook --execute uncollided.ipynb
If you later want to run the collided stage in parallel, use the companion solve_collided.py script after generating the uncollided file.
[1]:
import csv
import math
import os
import sys
from pathlib import Path
from mpi4py import MPI
size = MPI.COMM_WORLD.size
rank = MPI.COMM_WORLD.rank
if size != 1:
raise RuntimeError("This notebook must be run with exactly one MPI rank.")
if rank == 0:
print(f"Running with {size} MPI process.")
Running with 1 MPI process.
5.3.2.2. Import OpenSn
The notebook lives deep under doc/source/tutorials, so we append the repository root before importing pyopensn.
[2]:
repo_root = Path.cwd()
for candidate in (repo_root, *repo_root.parents):
if (candidate / "build/python").exists() and (candidate / "test/assets/mesh").exists():
repo_root = candidate
break
else:
raise RuntimeError("Could not locate the OpenSn repository root.")
sys.path.insert(0, str(repo_root / "build"))
from pyopensn.aquad import GLCProductQuadrature3DXYZ
from pyopensn.fieldfunc import (
FieldFunctionInterpolationLine,
FieldFunctionInterpolationPoint,
FieldFunctionInterpolationVolume,
)
from pyopensn.logvol import RPPLogicalVolume
from pyopensn.math import Vector3
from pyopensn.mesh import FromFileMeshGenerator
from pyopensn.solver import (
DiscreteOrdinatesProblem,
SteadyStateSourceSolver,
UncollidedProblem,
UncollidedSolver,
)
from pyopensn.source import PointSource
from pyopensn.xs import MultiGroupXS
5.3.2.3. Problem Setup
We use a single-material cube, a one-group total cross section of 0.8, and an isotropic scattering ratio of 0.45. The point source is placed safely inside the domain so the uncollided algorithm can localize it to one cell.
[3]:
mesh_file = repo_root / "test/assets/mesh/uncollided_cube_coarse.msh"
uncollided_file = Path("uncollided.h5").resolve()
source_location = (0.010, 0.012, 0.014)
sample_point = (0.024, 0.016, 0.008)
def make_mesh():
grid = FromFileMeshGenerator(filename=str(mesh_file)).Execute()
grid.SetUniformBlockID(0)
return grid
def make_xs():
xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=0.8, c=0.45)
return xs
def make_whole_domain():
return RPPLogicalVolume(
xmin=-0.001,
xmax=0.033,
ymin=-0.001,
ymax=0.033,
zmin=-0.001,
zmax=0.033,
)
def point_value(field_function, point):
interpolation = FieldFunctionInterpolationPoint()
interpolation.SetPointOfInterest(Vector3(*point))
interpolation.AddFieldFunction(field_function)
interpolation.Execute()
return interpolation.GetPointValue()
def volume_integral(field_function, logical_volume):
interpolation = FieldFunctionInterpolationVolume()
interpolation.SetOperationType("sum")
interpolation.SetLogicalVolume(logical_volume)
interpolation.AddFieldFunction(field_function)
interpolation.Execute()
return interpolation.GetValue()
line_start = Vector3(0.0, source_location[1], source_location[2])
line_end = Vector3(0.032, source_location[1], source_location[2])
plot_dir = Path("tutorial_output")
plot_dir.mkdir(exist_ok=True)
def csv_to_dict(filename):
data = {}
with open(filename, newline="") as csvfile:
reader = csv.DictReader(csvfile)
for row in reader:
for key, value in row.items():
data.setdefault(key, []).append(float(value))
return data
def export_line_data(field_function, csv_stem, start_point, end_point, num_points=200):
for csv_file in plot_dir.glob(f"{csv_stem}_*.csv"):
csv_file.unlink()
ffline = FieldFunctionInterpolationLine()
ffline.SetInitialPoint(start_point)
ffline.SetFinalPoint(end_point)
ffline.SetNumberOfPoints(num_points)
ffline.AddFieldFunction(field_function)
ffline.Execute()
ffline.ExportToCSV(str(plot_dir / csv_stem))
csv_file = next(plot_dir.glob(f"{csv_stem}_*.csv"))
return csv_to_dict(csv_file)
def analytic_uncollided_scalar_flux(distance, sigma_t=0.8, strength=1.0):
return strength * math.exp(-sigma_t * distance) / (4.0 * math.pi * distance * distance)
def relative_error(reference, value):
return abs(value - reference) / abs(reference)
5.3.2.4. Step 1: Generate the Uncollided File
The uncollided problem takes a list of point sources and writes the flux moments needed by the first-collision solve.
[4]:
try:
uncollided_file.unlink()
except FileNotFoundError:
pass
grid = make_mesh()
xs = make_xs()
whole_domain = make_whole_domain()
uncollided_problem = UncollidedProblem(
mesh=grid,
num_groups=1,
groupsets=[{"groups_from_to": [0, 0]}],
xs_map=[{"block_ids": [0], "xs": xs}],
point_sources=[PointSource(location=list(source_location), strength=[1.0])],
near_source=[whole_domain],
scattering_order=0,
)
uncollided_solver = UncollidedSolver(
problem=uncollided_problem,
file_name=str(uncollided_file),
progress_interval=25,
)
uncollided_solver.Initialize()
uncollided_solver.Execute()
if rank == 0:
print("Uncollided file:", uncollided_file)
OpenSn version 1.0.1
2026-06-09 18:03:54 Running OpenSn with 1 processes.
[0] FromFileMeshGenerator: Generating UnpartitionedMesh
Uncollided file:[0] Making unpartitioned mesh from Gmsh file /Users/dhawkins/opensn_dev/opensn/test/assets/mesh/uncollided_cube_coarse.msh (format v4.1 binary)
[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 74 of 735 (10%)
[0] 00:00:00.0 Surpassing cell 147 of 735 (20%)
[0] 00:00:00.0 Surpassing cell 221 of 735 (30%)
[0] 00:00:00.0 Surpassing cell 294 of 735 (40%)
[0] 00:00:00.0 Surpassing cell 368 of 735 (50%)
[0] 00:00:00.0 Surpassing cell 442 of 735 (60%)
[0] 00:00:00.0 Surpassing cell 515 of 735 (70%)
[0] 00:00:00.0 Surpassing cell 588 of 735 (80%)
[0] 00:00:00.0 Surpassing cell 662 of 735 (90%)
[0] 00:00:00.0 Surpassing cell 735 of 735 (100%)
[0] 00:00:00.0 Establishing cell boundary connectivity.
[0] 00:00:00.0 Done establishing cell connectivity.
[0] Done processing /Users/dhawkins/opensn_dev/opensn/test/assets/mesh/uncollided_cube_coarse.msh.
[0] Number of nodes read: 237
[0] Number of cells read: 735
[0] Number of cells per partition (max,min,avg) = 735,735,735
[0]
[0] Mesh statistics:
[0] Global cell count : 735
[0] Local cell count (avg,max,min): 735,735,735
[0] Ghost-to-local ratio (avg) : 0
[0]
[0] 00:00:00.0 Done setting block id 0 to all cells
[0]
[0] Initializing UncollidedProblem
[0]
[0] Initializing spatial discretization.
[0] 00:00:00.0 Computing unit integrals.
[0] 00:00:00.1 Ghost cell unit cell-matrix ratio: 0%
[0] 00:00:00.1 Cell matrices computed.
[0] 00:00:00.1 Initializing parallel arrays. G=1 M=1
[0] 00:00:00.1 Done with parallel arrays.
[0] Point source at [0.01 0.012 0.014] assigned to cell 52 with shape values [ 0.189239 0.420696 0.327311 0.0627551 ] and volume weight 1
[0] Point source has 1 subscribing cells on processor 0
[0] Point source has 1 global subscribing cells.
[0] 00:00:00.1 Initializing solver UncollidedSolver.
[0] 00:00:00.1 Starting solver execution UncollidedSolver.
[0] Starting uncollided transport for 1 physical source points. Progress will be reported every 25%.
[0] Uncollided progress: processing source point 1 / 1.
[0]
[0] Ray-tracing near-source region.
[0] Near-source current consistency diagnostic:
[0] Mismatched cells = 11 / 735
[0] Mismatched cell-group pairs = 11 / 735
[0] Mismatched cell fraction = 1.496599e-02
[0] Aggregate relative mismatch = 1.927983e-01
[0] Maximum relative mismatch = 5.058453e-01
[0] NearSourceCurrentMismatchedCellFraction=1.496599e-02
[0] NearSourceCurrentAggregateRelativeMismatch=1.927983e-01
[0] Uncollided progress: 1 / 1 source points (100.0%), elapsed 00:00:00, ETA 00:00:00.
[0] Global outflow consistency:
[0] Integrated vacuum outflow = 9.852756e-01
[0] Conservative outflow = 9.856949e-01
[0] Relative difference = 4.192985e-04
[0]
[0] Balance table:
[0] Removal rate = 1.430510e-02
[0] Production rate = 1.000000e+00
[0] Out-flow rate = 9.856949e-01
[0] Balance (Production - Loss) = 0.000000e+00
[0] Conservation error = 0.000000e+00
[0]
[0] 00:00:00.6 Finished solver execution UncollidedSolver.
/Users/dhawkins/opensn_dev/opensn/doc/source/tutorials/lbs/src_driven/uncollided/uncollided.h5
5.3.2.5. Compare Against the Analytic Uncollided Solution
For a unit isotropic point source in a homogeneous medium, the continuous uncollided scalar flux is exp(-\sigma_t r) / (4 \pi r^2) away from the source. We sample the OpenSn uncollided solution along the source-centered line and compare it against that reference. Because the continuous solution is singular at the source and the PWLD field is represented on a coarse tetrahedral mesh, we exclude a small near-source neighborhood from the line-based error summary and report the tutorial sample
point separately.
[5]:
uncollided_scalar_flux = uncollided_problem.GetScalarFluxFieldFunction()[0]
uncollided_line = export_line_data(
uncollided_scalar_flux,
"uncollided_line",
line_start,
line_end,
)
analytic_exclusion_radius = 8.0e-3
analytic_distances = []
analytic_values = []
open_sn_values = []
relative_errors = []
for x_value, phi_value in zip(uncollided_line["x"], uncollided_line["phi_g000_m00"]):
distance = abs(x_value - source_location[0])
if distance <= analytic_exclusion_radius:
continue
analytic_value = analytic_uncollided_scalar_flux(distance)
analytic_distances.append(distance)
analytic_values.append(analytic_value)
open_sn_values.append(phi_value)
relative_errors.append(relative_error(analytic_value, phi_value))
sample_distance = math.dist(sample_point, source_location)
analytic_sample_flux = analytic_uncollided_scalar_flux(sample_distance)
open_sn_sample_flux = point_value(uncollided_scalar_flux, sample_point)
if rank == 0:
print(f"AnalyticExclusionRadius={analytic_exclusion_radius:.12e}")
print(f"AnalyticComparisonPoints={len(relative_errors)}")
print(f"UncollidedLineMeanRelativeError={sum(relative_errors)/len(relative_errors):.12e}")
print(f"UncollidedLineMaxRelativeError={max(relative_errors):.12e}")
print(f"UncollidedSampleDistance={sample_distance:.12e}")
print(f"UncollidedSampleOpenSn={open_sn_sample_flux:.12e}")
print(f"UncollidedSampleAnalytic={analytic_sample_flux:.12e}")
print(
f"UncollidedSampleRelativeError="
f"{relative_error(analytic_sample_flux, open_sn_sample_flux):.12e}"
)
[0] Exported CSV file for field func "phi_g000_m00" to "tutorial_output/uncollided_line_phi_g000_m00.csv"
AnalyticExclusionRadius=8.000000000000e-03
AnalyticComparisonPoints=101
UncollidedLineMeanRelativeError=2.892774652525e-02
UncollidedLineMaxRelativeError=1.231411685501e-01
UncollidedSampleDistance=1.574801574802e-02
UncollidedSampleOpenSn=2.708264556451e+02
UncollidedSampleAnalytic=3.168597199861e+02
UncollidedSampleRelativeError=1.452796346063e-01
5.3.2.6. Optional: Visualize the Uncollided Flux
The notebook does not require matplotlib to run. If you want a quick line plot of the uncollided scalar flux, uncomment the next cell and run it in an environment where matplotlib is available.
[6]:
#import matplotlib.pyplot as plt
#
# uncollided_scalar_flux = uncollided_problem.GetScalarFluxFieldFunction()[0]
# uncollided_line = export_line_data(
# uncollided_scalar_flux,
# "uncollided_line",
# line_start,
# line_end,
# )
#
# plt.figure()
# plt.semilogy(uncollided_line["x"], uncollided_line["phi_g000_m00"])
# plt.xlabel("x (m)")
# plt.ylabel(r"$\phi_u(x)$")
# plt.title("Uncollided scalar flux along the source-centered line")
# plt.grid()
# plt.show()
5.3.2.7. Step 2: Solve the Collided Problem
The collided problem does not define the point source again. It reads the uncollided file through uncollided_flux and constructs the first-collision source internally.
[ ]:
grid = make_mesh()
xs = make_xs()
whole_domain = make_whole_domain()
quadrature = GLCProductQuadrature3DXYZ(
n_polar=2,
n_azimuthal=16,
scattering_order=0,
)
problem = DiscreteOrdinatesProblem(
mesh=grid,
num_groups=1,
groupsets=[
{
"groups_from_to": [0, 0],
"angular_quadrature": quadrature,
"angle_aggregation_type": "single",
"inner_linear_method": "petsc_gmres",
"gmres_restart_interval": 30,
"l_abs_tol": 1.0e-8,
"l_max_its": 200,
}
],
xs_map=[{"block_ids": [0], "xs": xs}],
uncollided_flux=str(uncollided_file),
)
solver = SteadyStateSourceSolver(problem=problem, compute_balance=True)
solver.Initialize()
solver.Execute()
5.3.2.8. Post-Process the Total Flux
After the steady-state solve converges, the reported field functions contain the total scalar flux. We sample one point and compute the domain integral.
[8]:
scalar_flux = problem.GetScalarFluxFieldFunction()[0]
point_flux = point_value(scalar_flux, sample_point)
total_integral = volume_integral(scalar_flux, whole_domain)
balance = solver.ComputeBalanceTable()["balance"]
if rank == 0:
print(f"Total scalar flux at {sample_point}: {point_flux:.12e}")
print(f"Total scalar-flux integral: {total_integral:.12e}")
print(f"Balance residual: {balance:.12e}")
Total scalar flux at (0.024, 0.016, 0.008): 2.735814472835e+02
Total scalar-flux integral: 1.798893012550e-02
Balance residual: -2.220446049250e-15
5.3.2.9. Optional: Visualize the Total Flux
The notebook does not require matplotlib to run. If you want a quick line plot of the total scalar flux, uncomment the next cell and run it in an environment where matplotlib is available.
[9]:
# import matplotlib.pyplot as plt
#
# total_line = export_line_data(
# scalar_flux,
# "total_line",
# line_start,
# line_end,
# )
#
# plt.figure()
# plt.semilogy(total_line["x"], total_line["phi_g000_m00"])
# plt.xlabel("x (m)")
# plt.ylabel(r"$\phi(x)$")
# plt.title("Total scalar flux along the source-centered line")
# plt.grid()
# plt.show()