12. 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

12.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.

12.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.

12.3. Constructor Inputs

12.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.

12.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.

12.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.

12.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.

12.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.

12.3.6. point_sources and volumetric_sources

These are optional source objects added at construction time.

12.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.

12.3.8. options

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

12.3.9. sweep_type

This selects the sweep type:

• "AAH": the aggregated sweeper with cycle breaking

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

12.3.10. use_gpus

This requests GPU acceleration where supported.

12.3.11. time_dependent

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

12.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()

12.5. Problem Options

The Python API exposes problem options through the constructor and SetOptions().

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.

12.5.1. Setting options after construction

Problem options can also be changed after construction:

phys.SetOptions(
    verbose_inner_iterations=False,
    max_ags_iterations=200,
    ags_tolerance=1.0e-8,
)

This is useful for parameter studies or staged workflows where the problem definition is reused.

12.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.

12.7. Problem Modes

12.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.

12.7.2. use_gpus

use_gpus requests GPU acceleration for supported sweep paths.

Current restrictions:

• only "AAH" is supported for GPU use,

• curvilinear problems do not support GPU acceleration,

• time-dependent problems do not support GPU acceleration,

• adjoint 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.

12.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).

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

12.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.

12.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.

12.8.2. Source support

The uncollided generator supports:

• explicit pyopensn.source.PointSource objects.

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.

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.

12.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 and folds attenuation paths back through the physical mesh. 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.

The near-source calculation independently projects the ray-traced volume flux and integrates ray-traced face currents. These two quadratures generally do not satisfy exact cell balance at finite resolution. OpenSn reports their relative mismatch but preserves both projections; rescaling outgoing currents cell by cell can recursively amplify quadrature error along long streaming paths. The HDF5 balance metadata uses the projected removal and the conservative global outflow remainder. The directly integrated vacuum outflow is also printed as a consistency diagnostic.

12.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.

12.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.

12.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.

12.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.

12.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. The generator reports the integrated and conservative effective outflows and stores the conservative value used for combined balance accounting; this correction does not rescale the uncollided flux moments.

12.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.

12.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.

12.11. Balance and Leakage

The Cartesian discrete-ordinates problem exposes two important diagnostics:

• ComputeBalance()

• ComputeLeakage()

12.11.1. ComputeBalance()

This computes the particle balance for the problem and returns a dictionary of balance terms.

It is useful for:

• checking whether the run is physically consistent,

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

• validating test and regression problems.

12.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.

12.12. Writing and Reading Transport State

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

• WriteFluxMoments()

• ReadFluxMoments()

• CreateAndWriteSourceMoments()

• ReadSourceMoments()

• ReadFluxMomentsAndMakeSourceMoments()

• WriteAngularFluxes()

• ReadAngularFluxes()

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.

12.13. Updating the Problem In Place

Several parts of the problem can be updated after construction.

The most important methods are:

• SetOptions()

• SetPointSources()

• SetVolumetricSources()

• SetBoundaryOptions()

• SetAdjoint()

• SetTimeDependentMode()

12.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.

12.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",
)

12.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

12.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.