11. Discrete-Ordinates Problems

The discrete-ordinates problem object is the center of OpenSn inputs. It gathers the mesh, materials, groupsets, sources, boundary conditions, problem options, and transport state into one object that can then be driven by a solver.

The two main problem classes exposed in Python are:

• pyopensn.solver.DiscreteOrdinatesProblem

• pyopensn.solver.DiscreteOrdinatesCurvilinearProblem

For first-collision calculations, OpenSn also provides:

• pyopensn.solver.UncollidedProblem

• pyopensn.solver.UncollidedSolver

11.1. Overview

Use pyopensn.solver.DiscreteOrdinatesProblem for Cartesian discrete-ordinates problems. This is the standard problem class for most steady-state, transient, and k-eigenvalue workflows.

The usual Cartesian construction pattern is:

phys = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=num_groups,
    groupsets=groupsets,
    xs_map=xs_map,
    boundary_conditions=boundary_conditions,
    point_sources=point_sources,
    volumetric_sources=volumetric_sources,
    options={...},
    sweep_type="AAH",
    use_gpus=False,
    time_dependent=False,
)

The problem object owns:

• the mesh and material mapping,

• the scalar and angular transport state,

• the groupset definitions,

• the current forward or adjoint mode,

• the current steady-state or time-dependent mode,

• the field-function interface associated with the problem.

In practice, that means the problem is where the transport model is assembled:

• the spatial domain and geometry ids come from the mesh,

• the material assignment comes from xs_map,

• the energy groupings and inner iteration setup come from groupsets,

• driving terms come from sources and boundary conditions,

• global transport behavior comes from problem-level options.

Note

A solver does not define the transport model. The problem does. The solver only decides how that already-defined model is advanced or converged.

11.2. Constructor Summary

The main constructor inputs are:

• mesh: a pyopensn.mesh.MeshContinuum

• num_groups: total number of energy groups

• groupsets: one or more groupset dictionaries

• xs_map: block-id to cross-section mapping

• boundary_conditions: optional boundary specifications

• point_sources: optional point sources

• volumetric_sources: optional volumetric sources

• options: problem-level settings

• sweep_type: "AAH" or "CBC"

• use_gpus: request GPU sweep support where available

• time_dependent: whether the problem starts in time-dependent mode

A complete problem definition often looks like:

phys = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=2,
    groupsets=groupsets,
    xs_map=[
        {"block_ids": [0], "xs": fuel_xs},
        {"block_ids": [1], "xs": moderator_xs},
    ],
    boundary_conditions=boundary_conditions,
    volumetric_sources=volumetric_sources,
    options={
        "verbose_inner_iterations": False,
        "verbose_outer_iterations": True,
    },
    sweep_type="AAH",
    time_dependent=False,
)

Each of these pieces is discussed elsewhere in the guide, but it is useful to see them in one place.

11.3. Constructor Inputs

11.3.1. mesh

This is the transport mesh.

It must already exist and already carry the block ids and boundary labels that the rest of the problem will use.

11.3.2. num_groups

This is the total number of energy groups in the problem.

It must be consistent with:

• the loaded cross sections,

• the groupset definitions,

• any group-wise source data.

11.3.3. groupsets

This is the list of groupset dictionaries that define angular quadratures, inner iteration, angle aggregation, and optional DSA configuration.

Groupsets are detailed separately in Groupsets, but from the problem perspective they define how the energy range is partitioned for transport iterations.

11.3.4. xs_map

This maps mesh block ids to cross-section objects.

Every block id used in the transport region should be assigned an XS object.

11.3.5. boundary_conditions

This is the list of transport boundary-condition dictionaries. The names used here must match the boundary ids on the mesh.

11.3.6. point_sources and volumetric_sources

These are optional source objects added at construction time.

11.3.7. uncollided_flux

This is an optional HDF5 file generated by pyopensn.solver.UncollidedSolver. Supplying it enables the first-collision workflow described in Uncollided and First-Collision Transport below.

11.3.8. options

This dictionary holds problem-level behavior that is not naturally part of the mesh, source, or groupset definitions.

11.3.9. sweep_type

This selects the sweep type:

• "AAH": the aggregated sweeper with cycle breaking

• "CBC": the cell-by-cell sweeper

11.3.10. use_gpus

This requests GPU acceleration where supported.

11.3.11. time_dependent

This puts the problem into time-dependent mode at construction time.

11.4. Shared LBS Interface

Because DiscreteOrdinatesProblem derives from pyopensn.solver.LBSProblem, it also exposes the common LBS problem methods such as:

• GetScalarFluxFieldFunction()

• CreateFieldFunction()

• ComputeFissionRate()

• WriteFluxMoments()

• ReadFluxMoments()

• SetPointSources()

• SetVolumetricSources()

• SetBoundaryOptions()

• SetAdjoint()

11.5. Problem Options

The Python API exposes problem options through the options constructor argument. There is currently no general-purpose setter for changing arbitrary problem options after construction; individual aspects of the problem (sources, boundary conditions, adjoint mode, and time-dependent mode) have their own dedicated update methods instead.

The most important options for everyday use are:

• adjoint

• save_angular_flux

• verbose_inner_iterations

• verbose_outer_iterations

• max_ags_iterations

• ags_tolerance

• ags_convergence_check

• field_function_prefix_option

• field_function_prefix

There are also restart, precursor, and message-size options for more specialized workflows.

Restart-related options include:

• restart_writes_enabled: enable restart dump writes.

• write_restart_path: file stem used when writing restart dumps. OpenSn appends the MPI rank and .restart.h5.

• read_restart_path: file stem used when reading a full restart. This is for continuing the same type of solve from restart state.

• read_initial_condition_path: file stem used when reading restart data as an initial condition. pyopensn.solver.TransientSolver can use a steady-state restart this way, then switch the problem to time-dependent mode.

• write_angular_flux_to_restart: include stored angular fluxes in restart dumps when save_angular_flux=True. This is required for full transient continuation restarts, but optional when a steady-state restart is used only as a transient initial condition.

• write_delayed_psi_to_restart: include delayed sweep angular-flux buffers. Full continuation restarts require these buffers whenever the problem has delayed sweep angular state, including partitioned parallel, reflected-boundary, and cyclic-sweep cases. These buffers are optional for the steady-state-restart-as-transient-initial-condition workflow.

Warning

Restart files are rank-layout specific. A restart written with one MPI rank count should be read with the same rank count and a compatible problem definition. Changing from serial to parallel, or from one partition count to another, is not a supported restart workflow.

Note

Problem options are where users should look for global transport behavior. If a setting changes how the entire problem behaves rather than how one groupset behaves, it usually belongs here.

11.5.1. Setting options after construction

Problem-wide options values are fixed at construction time. For parameter studies or staged workflows, rebuild the problem object with a new options dictionary rather than mutating an existing one. Aspects of the problem that do have dedicated post-construction setters are covered elsewhere in this guide: SetPointSources(), SetVolumetricSources(), SetBoundaryOptions(), SetAdjoint(), and SetTimeDependentMode() / SetSteadyStateMode().

11.6. sweep_type

The Python API exposes:

• "AAH"

• "CBC"

Example:

sweep_type="AAH"

If omitted, the default is "AAH".

Operationally, the two sweep types are different:

• "AAH" is the more general aggregated sweeper and should be treated as the default production choice.

• "CBC" is a cell-by-cell sweeper that preserves exact cell-to-cell dependencies.

The practical differences are:

• AAH has explicit delayed-angular-flux machinery for cycle handling.

• AAH can break both local and inter-partition sweep cycles by removing feedback-arc-set edges and lagging the corresponding angular-flux dependencies.

• CBC does not support local sweep cycles.

Note

In the AAH implementation, the lagged data is tied specifically to cycle-breaking dependencies. It is not a general statement that all angular fluxes are always lagged.

Practical recommendation:

• AAH remains the default production choice and is the safer option for most users, particularly for problems with cyclic sweep dependencies.

• Both AAH and CBC support time-dependent (transient) mode.

• Choose CBC only when the sweep graph is known to be acyclic or when you have verified it meets the acyclicity requirement for your specific problem.

11.7. Problem Modes

11.7.1. time_dependent

If time_dependent=True, the problem starts in time-dependent mode.

Important requirement:

• time-dependent operation requires options={"save_angular_flux": True}

Example:

phys = DiscreteOrdinatesProblem(
    ...,
    time_dependent=True,
    options={"save_angular_flux": True},
)

This requirement exists because transient updates need access to angular-flux state from one timestep to the next.

This option changes the problem mode, not just a solver setting. That is why it belongs on the problem object and must be consistent with the solver used later.

11.7.2. use_gpus

use_gpus requests GPU acceleration for supported sweep paths.

Current restrictions:

• both "AAH" and "CBC" sweep types support GPU acceleration,

• curvilinear problems do not support GPU acceleration,

• time-dependent problems do not support GPU acceleration.

Most users should treat this as a deployment choice after the base problem is already running correctly on the CPU.

11.7.3. Adjoint Mode

Adjoint mode is controlled at the problem level because it changes the meaning of the transport problem itself, not just the iterative algorithm.

It can be set at construction time through options or later with SetAdjoint(True).

Adjoint calculations are supported only for Cartesian geometries. Cylindrical and spherical coordinate systems require a dedicated discrete-adjoint treatment of their angular-redistribution terms and are rejected.

Because this is a fundamental change in problem interpretation, users should be deliberate about when they switch it on.

11.8. Uncollided and First-Collision Transport

First-collision transport splits the total angular flux into uncollided and collided components:

\[\Psi = \Psi^u + \Psi^c.\]

pyopensn.solver.UncollidedProblem defines the uncollided transport configuration, and pyopensn.solver.UncollidedSolver computes and stores its flux moments. A subsequent pyopensn.solver.DiscreteOrdinatesProblem reads those moments, constructs the first-collision scattering and fission source, solves for the collided component, and adds the uncollided component back to the converged transport state.

This treatment is useful for localized sources in optically thin or weakly scattering regions where a conventional angular quadrature would otherwise produce strong ray effects.

11.8.1. Two-stage workflow

The uncollided calculation follows the same problem/solver lifecycle as other OpenSn transport calculations. Construct the problem and solver, then call Initialize and Execute:

Uncollided generation is restricted to Cartesian two- and three-dimensional meshes and must run with exactly one MPI rank.

The solver reports source-point progress every 5 percent by default. Each update includes the number of completed source points, elapsed time, and estimated remaining time. Set progress_interval on pyopensn.solver.UncollidedSolver to another percentage, or to zero to disable progress reporting.

Internal threading used by the uncollided solver is capped by the environment variable OPENSN_NUM_THREADS. If the variable is unset or invalid, OpenSn uses 1 thread. The current uncollided implementation applies this cap to reflected-image projection and to the bulk-sweep group solve.

from pyopensn.aquad import GLCProductQuadrature2DXY
from pyopensn.logvol import RPPLogicalVolume
from pyopensn.solver import (
    DiscreteOrdinatesProblem,
    SteadyStateSourceSolver,
    UncollidedProblem,
    UncollidedSolver,
)
from pyopensn.source import PointSource

near_source = RPPLogicalVolume(
    xmin=-0.25,
    xmax=0.25,
    ymin=-0.25,
    ymax=0.25,
    infz=True,
)
point_source = PointSource(
    location=[0.0, 0.0, 0.0],
    strength=[1.0],
)

uncollided_file = "uncollided_flux.h5"
uncollided = UncollidedProblem(
    mesh=mesh,
    num_groups=1,
    groupsets=[{"groups_from_to": [0, 0]}],
    xs_map=xs_map,
    point_sources=[point_source],
    near_source=[near_source],
    scattering_order=1,
)
uncollided_solver = UncollidedSolver(
    problem=uncollided,
    file_name=uncollided_file,
)
uncollided_solver.Initialize()
uncollided_solver.Execute()

quadrature = GLCProductQuadrature2DXY(
    n_polar=2,
    n_azimuthal=24,
    scattering_order=1,
)
problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=1,
    groupsets=[
        {
            "groups_from_to": [0, 0],
            "angular_quadrature": quadrature,
        }
    ],
    xs_map=xs_map,
    uncollided_flux=uncollided_file,
)
solver = SteadyStateSourceSolver(problem=problem, compute_balance=True)
solver.Initialize()
solver.Execute()

Do not also attach point_source to the collided problem. Its contribution is already represented by the uncollided flux and the generated first-collision source. Independent sources that are not represented in the uncollided file may still be attached normally.

11.8.2. Source support

The uncollided generator supports:

• explicit isotropic pyopensn.source.PointSource objects.

The group-wise point-source strengths are total angle-integrated emission rates. Angular distributions are not currently supported by the uncollided problem.

Every point source requires a corresponding entry in near_source. The two lists are matched by position, so near_source is required whenever point_sources is nonempty.

For now, uncollided generation requires each point source to lie strictly inside a single cell. A source located exactly on a face, edge, or vertex is rejected as unsupported.

Important

Point-source placement strongly affects spatial accuracy. Whenever the physical model permits, construct the mesh so that each source lies well inside its containing cell and near the cell centroid. A source near a face, edge, or vertex produces highly uneven local basis weights and makes the near-singular flux and outgoing currents difficult to integrate accurately. The conservation correction restores cell balance, but it cannot recover an unresolved local flux shape.

Do not move a physically fixed source merely to improve the discretization. Instead, align or locally refine the mesh around it. Check mesh convergence of the physical quantities of interest and monitor the reported near-source correction. OpenSn uses fixed fourth-order finite-element volume and face quadratures for this calculation; their order is not currently user-configurable. A large correction indicates that these spatial integrations do not adequately resolve the source on the current mesh.

In two-dimensional problems, a point source represents the two-dimensional transport Green’s function, equivalently a line source per unit out-of-plane depth. Its free-space uncollided scalar flux scales as \(1/(2 \pi r)\), not as the three-dimensional \(1/(4 \pi r^2)\) point-source field.

Finite-volume sources can be approximated by the input author with weighted point sources. For example, a quadrature approximation uses

\[Q_{g,q} = q_g(\mathbf{x}_q) J_q w_q\]

and passes the resulting weighted points through point_sources. The uncollided generator does not consume pyopensn.source.VolumetricSource objects directly. Choose the source quadrature and mesh together so that the weighted points lie well inside their cells. Verify convergence with respect to both the source quadrature and the mesh.

11.8.3. Reflecting boundaries

The uncollided and collided stages must specify the same reflecting boundary conditions. The uncollided generator represents each reflection with image sources. To compute attenuation, OpenSn reflects any portion of each image-source ray outside the physical domain back into the domain, then evaluates attenuation along the resulting physical path. Each image source is ray traced to every finite-element volume quadrature point and projected directly into the spatial discretization. With \(N\) reflecting symmetry planes, each physical source point produces \(2^N-1\) image contributions.

This construction supports up to three planar, mutually orthogonal symmetry planes, such as xmin, ymin, and zmin:

boundary_conditions = [
    {"name": "xmin", "type": "reflecting"},
    {"name": "ymin", "type": "reflecting"},
    {"name": "zmin", "type": "reflecting"},
]

uncollided = UncollidedProblem(
    mesh=mesh,
    num_groups=num_groups,
    groupsets=groupsets,
    xs_map=xs_map,
    point_sources=point_sources,
    near_source=near_source_regions,
    boundary_conditions=boundary_conditions,
    scattering_order=scattering_order,
)
uncollided_solver = UncollidedSolver(
    problem=uncollided,
    file_name="uncollided_flux.h5",
)
uncollided_solver.Initialize()
uncollided_solver.Execute()

The HDF5 file records the reflecting boundary IDs, and the collided problem rejects files generated with a different reflector set.

Physical source points use the cell-wise near-source conservation correction of Woodsford, Ragusa, and Morel (2026) [1]. A shared scale factor is applied to their projected flux and outgoing currents. Its reported magnitude indicates how well the raw volume and face integrations agree; see the derivation in Iterative Solution Algorithms. The correction does not enforce nonnegative individual PWLD coefficients.

Reflected image sources are an OpenSn extension to the published algorithm and are projected without the near-source correction. The HDF5 balance metadata stores absorption computed using the volume quadrature and outflow computed using a separate face quadrature. Because these integrations have no conservation correction, their finite-quadrature errors appear in the reported balance residual.

[1]

C. Woodsford, J. C. Ragusa, and J. E. Morel, “Sweep-based uncollided-flux treatment on unstructured grids,” Progress in Nuclear Energy, vol. 200, 106524, 2026, https://doi.org/10.1016/j.pnucene.2026.106524.

11.8.4. Moment order

scattering_order on pyopensn.solver.UncollidedProblem is the maximum spherical-harmonic order written to the HDF5 file. It must be at least the scattering order used by the collided problem.

The scalar moment is stored as 0,0. Higher moments are stored by ell,m name and are accumulated over all explicit and generated source points. The published sweep-based method considered isotropic scattering and did not require these higher moments. OpenSn extends the method by evaluating \(Y_{\ell,m}(\vec{\Omega}_s)\phi^u_g\) at the finite-element volume quadrature points, where \(\vec{\Omega}_s\) points from the physical or image source to the evaluation point, and projecting that product onto the PWLD basis. This uses the analytic ray direction rather than the collided problem’s \(S_n\) angular quadrature.

11.8.5. Serial generation and parallel reuse

Uncollided generation must run with exactly one MPI rank. The resulting HDF5 file is partition-independent: every rank in the collided calculation reads the same serial file and extracts the cells it owns. Consequently, the collided calculation may use one or multiple MPI ranks.

For a parallel workflow, run the two stages separately:

1. Generate the HDF5 file with one rank.

2. Run the collided input with the desired number of ranks.

When generating the uncollided file, per-rank internal threading remains controlled by OPENSN_NUM_THREADS and defaults to 1. This avoids oversubscribing MPI jobs unless the user explicitly opts in to additional threads.

11.8.6. File compatibility

The collided problem validates the HDF5 file before constructing the first-collision source. The generating and consuming problems must have:

• the same number of energy groups,

• the same global cell IDs and cell count,

• the same cell-node layout and node coordinates,

• matching total cross sections in every cell and group,

• the same reflecting boundary set,

• sufficient stored moment order for the collided scattering order.

The mesh may be repartitioned for the collided calculation because matching is performed using global cell IDs. Cross sections cannot be replaced through SetXSMap after a problem has loaded an uncollided file.

11.8.7. Supported problem mode

The collided use of uncollided_flux is supported for steady-state forward fixed-source calculations. It is not supported for:

• time-dependent calculations,

• adjoint calculations,

• k-eigenvalue solvers.

11.8.8. Flux and balance interpretation

After pyopensn.solver.UncollidedSolver executes, the field functions returned by pyopensn.solver.UncollidedProblem contain the uncollided component \(\Phi^u\).

After the collided solver converges, the scalar flux state and field functions returned by pyopensn.solver.DiscreteOrdinatesProblem contain the recombined total:

\[\Phi = \Phi^u + \Phi^c.\]

The uncollided HDF5 file also stores its production, removal, and outflow rates. The steady-state solver incorporates the uncollided production and outflow when reporting the combined problem balance. Reflected image sources are projected directly from ray traces evaluated at every finite-element volume quadrature point and are not near-source corrected. The generator stores the integrated vacuum outflow for combined balance accounting. The near-source correction rescales the physical source’s uncollided flux moments and outgoing currents together, as specified by the published algorithm.

11.9. Field-Function Interface

The problem object is also where users access transport outputs.

The main field-function accessors are:

• GetScalarFluxFieldFunction()

• CreateFieldFunction()

• GetAngularFieldFunctionList()

Example:

scalar_ffs = phys.GetScalarFluxFieldFunction()
power_ff = phys.CreateFieldFunction("power_generation", "power")

Note

In LBS workflows, users normally work with the field-function objects returned by these accessors directly. They are created from the current transport state when requested. If the transport state changes later, call Update() on an existing updateable field function or create a fresh field function from the current state.

11.10. Angular-Flux Access

The discrete-ordinates problem also exposes angular-flux data directly through GetPsi().

GetPsi() returns a list of NumPy arrays copied from the current angular- flux storage. Mutating the returned arrays does not mutate the problem.

This is a specialized interface and is mainly useful when:

• a workflow needs direct access to angular-flux data,

• a custom analysis step is easier to write in Python,

• the problem was configured with save_angular_flux=True for the needed workflow.

Note

Angular flux is much larger than scalar flux. Many workflows never need it. Only enable and use it when the calculation actually requires it.

11.11. Balance and Leakage

The Cartesian discrete-ordinates problem and its solvers expose two important diagnostics:

• ComputeBalanceTable()

• ComputeLeakage()

11.11.1. ComputeBalanceTable()

This is a solver-level method, exposed on pyopensn.solver.SteadyStateSourceSolver, pyopensn.solver.TransientSolver, pyopensn.solver.PowerIterationKEigenSolver, and pyopensn.solver.NonLinearKEigenSolver. It computes the particle balance for the problem and returns a dictionary of balance terms.

Example:

balance = solver.ComputeBalanceTable()
print(balance["balance"])

It is useful for:

• checking whether the run is physically consistent,

• confirming source, absorption, leakage, and production trends,

• validating test and regression problems.

11.11.2. ComputeLeakage()

This computes boundary leakage and returns a dictionary mapping boundary names to group-wise NumPy arrays.

Example:

leakage = phys.ComputeLeakage(["xmin", "xmax"])

Important requirement:

• leakage computation requires save_angular_flux=True

This requirement exists because leakage is derived from the outgoing angular flux.

11.12. Writing and Reading Transport State

The problem object also exposes file-based state helpers such as:

• WriteFluxMoments()

• ReadFluxMoments()

• CreateAndWriteSourceMoments()

• ReadSourceMoments()

• ReadFluxMomentsAndMakeSourceMoments()

• WriteAngularFluxes()

• ReadAngularFluxes()

• WriteSurfaceAngularFluxes()

• ReadSurfaceAngularFluxes()

These are useful for:

• restart-like workflows,

• response studies,

• source-driven workflows that reuse previously computed fields.

For restartable solver state, prefer the restart options documented above over manual angular-flux files. Manual angular-flux files are useful when the workflow explicitly wants to manage angular fluxes, but restart dumps also carry flux moments, time metadata, precursor state, and solver-specific restart data.

11.13. Updating the Problem In Place

Several parts of the problem can be updated after construction.

The most important methods are:

• SetPointSources()

• SetVolumetricSources()

• SetBoundaryOptions()

• SetAdjoint()

• SetTimeDependentMode()

11.13.1. Updating a Problem After Construction

One of the strengths of the problem API is that major parts of the model can be replaced without rebuilding the entire problem object.

Examples:

phys.SetVolumetricSources(
    clear_volumetric_sources=True,
    volumetric_sources=[new_source],
)

phys.SetBoundaryOptions(
    clear_boundary_conditions=True,
    boundary_conditions=[new_boundary],
)

phys.SetAdjoint(True)

This makes the problem object usable in:

• source studies,

• boundary-condition studies,

• forward/adjoint comparisons,

• transient driver loops that change forcing terms.

11.14. DiscreteOrdinatesCurvilinearProblem

pyopensn.solver.DiscreteOrdinatesCurvilinearProblem is the curvilinear companion to the Cartesian problem class.

It uses the same general construction pattern, but currently requires:

• a suitable curvilinear mesh,

• coord_system=2 for cylindrical coordinates,

• a compatible quadrature and solver setup.

Important current limitations:

• the curvilinear solver only supports cylindrical geometries,

• GPU acceleration is not supported,

• users should treat it as a more specialized path than the standard Cartesian problem.

Example:

phys = DiscreteOrdinatesCurvilinearProblem(
    mesh=mesh,
    coord_system=2,
    num_groups=num_groups,
    groupsets=groupsets,
    xs_map=xs_map,
    sweep_type="AAH",
)

11.15. Typical Construction Patterns

A few patterns show up repeatedly:

• source-driven steady-state problem plus SteadyStateSourceSolver

• time-dependent problem plus TransientSolver

• multiplication problem plus PowerIterationKEigenSolver

• multiplication problem plus NonLinearKEigenSolver

11.16. Practical Guidance

As a rule:

• use DiscreteOrdinatesProblem unless the curvilinear class is specifically needed,

• keep sweep_type="AAH" unless there is a clear reason to choose CBC,

• treat save_angular_flux as a capability switch with a memory cost,

• build the simplest correct problem first, then add extra options,

• remember that the problem object remains central even after the solver is constructed because it owns the transport state and output interface.

Note

Many solver issues are actually problem-definition issues. If something looks unstable or physically wrong, it is often worth reviewing the problem object first: mesh labels, xs_map, groupsets, source definitions, and problem options.