25. Post Processors

For most OpenSn workflows, post processing starts with field functions. A field function is a mesh-based representation of transport data that can be exported or interpolated from Python.

The important practical point is that field functions are created from the current solver state when requested. They are snapshots, not continuously updated solver-owned views. Some field functions can refresh that same snapshot object explicitly through pyopensn.fieldfunc.FieldFunctionGridBased.Update().

25.1. Overview

The OpenSn field function interfaces are:

In practice, most workflows are:

solve the problem,
create the field functions needed from the current state,
export them or evaluate them with interpolation objects.

Note

If the solver state changes, a field function created earlier does not update itself automatically after a later solve or timestep. Either create a new field function from the current state, or call Update() on the existing field function when CanUpdate() returns True.

Note

For LBS workflows, the important field-function creators are GetScalarFluxFieldFunction(), CreateFieldFunction(), and, for discrete ordinates problems, GetAngularFieldFunctionList(). These methods return fresh field functions from the current state; they do not rely on a persistent field-function cache.

25.2. Updating Existing Field Functions

Field functions returned by LBS problem accessors are updateable snapshots. They hold their own field-vector data, but also know how to refresh that data from the owning problem while the owning problem is still alive.

Use pyopensn.fieldfunc.FieldFunctionGridBased.CanUpdate() before refreshing a field function that may have outlived its problem:

scalar_ff = phys.GetScalarFluxFieldFunction()[0]

solver.Advance()
if scalar_ff.CanUpdate():
    scalar_ff.Update()

Calling Update() refreshes the same field-function object. Interpolators and export calls that already reference that object will see the refreshed data the next time they execute.

This is mainly useful in timestep loops or repeated-solve workflows where reusing the same interpolation or export setup is clearer than reconstructing field functions each time.

25.3. Scalar Flux Field Functions

25.3.1. `GetScalarFluxFieldFunction`

Use pyopensn.solver.LBSProblem.GetScalarFluxFieldFunction() to create scalar-flux or flux-moment field functions from the current scalar-flux state.

The common case is scalar flux only:

scalar_ffs = phys.GetScalarFluxFieldFunction()

This returns one field function per energy group, each representing the zeroth moment of the flux.

If higher moments are needed too:

ff_by_group_and_moment = phys.GetScalarFluxFieldFunction(
    only_scalar_flux=False,
)

In that form:

result[g][m] is the field function for group g and moment m

Note

Most visualization and response workflows begin with scalar flux. Higher moments are mainly for diagnostics and specialized analysis.

25.4. Derived Field Functions

25.4.1. `CreateFieldFunction`

Use pyopensn.solver.LBSProblem.CreateFieldFunction() to create a named scalar field function derived from a 1D cross section or from the special case "power".

Examples:

fission_rate_ff = phys.CreateFieldFunction("fission_rate", "sigma_f")
power_ff = phys.CreateFieldFunction("power_generation", "power")

For a built-in or custom 1D XS name, this creates the field:

sum_g xs[g] * phi_g

at each spatial point.

If a power-normalized view is needed:

power_ff = phys.CreateFieldFunction(
    "power_generation_norm",
    "power",
    power_normalization_target=1.0,
)

Important points:

name is the field-function name assigned to the returned object
xs_name can be a built-in 1D XS name, a custom XS name, or "power"
power_normalization_target is optional and affects only the returned field function

Note

power_normalization_target applies a power-based scaling to the returned field function only. It does not rescale the solver’s internal flux vectors, and it does not change field functions created earlier. If you want a normalized power field, use CreateFieldFunction("name", "power", power_normalization_target=...). The same argument can also be used for other derived fields such as sigma_f * phi when a power-normalized post-processed view is desired.

25.5. Angular Flux Field Functions

25.5.1. `GetAngularFieldFunctionList`

Use pyopensn.solver.DiscreteOrdinatesProblem.GetAngularFieldFunctionList() to create field functions for selected angular-flux components.

Example:

ang_ffs = phys.GetAngularFieldFunctionList(groups=[0], angles=[0])

Important requirements:

this is available on pyopensn.solver.DiscreteOrdinatesProblem
angular-flux storage must be enabled with save_angular_flux=True
the returned field functions are snapshots created from the currently stored angular flux
if the stored angular flux changes later, call Update() on the existing angular field functions or create new ones from the current state

Important

For transient problems, save_angular_flux=True is not optional. It is required by the transient solver itself, and it is also required if you want to create angular-flux field functions or query GetPsi() during the transient run.

This is mainly useful for:

angular diagnostics,
leakage studies,
checking directional structure in difficult problems.

Note

Angular flux is much larger than scalar flux. Only enable angular-flux storage and create angular field functions when the workflow actually needs them.

Scalar-flux field functions use names based on group and moment, for example:

phi_g000_m00
phi_g001_m00
phi_g000_m01

If the problem uses a field-function prefix, that prefix is applied to these names as well.

25.6. Exporting Field Functions

The main Python export routine is:

pyopensn.fieldfunc.FieldFunctionGridBased.ExportMultipleToPVTU()

Export scalar flux:

from pyopensn.fieldfunc import FieldFunctionGridBased

scalar_ffs = phys.GetScalarFluxFieldFunction()
FieldFunctionGridBased.ExportMultipleToPVTU(
    scalar_ffs,
    "scalar_flux",
)

Export a derived field:

from pyopensn.fieldfunc import FieldFunctionGridBased

power_ff = phys.CreateFieldFunction("power_generation", "power")
FieldFunctionGridBased.ExportMultipleToPVTU(
    [power_ff],
    "power",
)

Export several fields together:

from pyopensn.fieldfunc import FieldFunctionGridBased

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

FieldFunctionGridBased.ExportMultipleToPVTU(
    scalar_ffs + [power_ff],
    "transport_outputs",
)

Note

Export after the solve, not before. For transient problems, export after the completed timestep whose state you want to visualize.

25.7. Field-Function Interpolation

25.7.1. Point Interpolation

Use pyopensn.fieldfunc.FieldFunctionInterpolationPoint to evaluate a field function at a point.

Example:

from pyopensn.fieldfunc import FieldFunctionInterpolationPoint

ffi = FieldFunctionInterpolationPoint()
ffi.AddFieldFunction(phys.GetScalarFluxFieldFunction()[0])
ffi.Initialize()
ffi.Execute()
value = ffi.GetPointValue()

This is useful for detector-like spot checks or comparison against an analytic value.

25.7.2. Line Interpolation

Use pyopensn.fieldfunc.FieldFunctionInterpolationLine to sample a field function along a line.

Example:

from pyopensn.fieldfunc import FieldFunctionInterpolationLine
from pyopensn.math import Vector3

ffi = FieldFunctionInterpolationLine()
ffi.AddFieldFunction(phys.GetScalarFluxFieldFunction()[0])
ffi.SetInitialPoint(Vector3(0.0, 0.0, 0.0))
ffi.SetFinalPoint(Vector3(10.0, 0.0, 0.0))
ffi.SetNumberOfPoints(101)
ffi.Initialize()
ffi.Execute()
ffi.ExportToCSV("centerline")

This is useful for profiles, attenuation curves, and one-dimensional comparisons across runs.

25.7.3. Volume Interpolation

Use pyopensn.fieldfunc.FieldFunctionInterpolationVolume for region averages, sums, maxima, and function-weighted variants over a logical volume.

Example:

from pyopensn.fieldfunc import FieldFunctionInterpolationVolume

ffi = FieldFunctionInterpolationVolume()
ffi.AddFieldFunction(phys.GetScalarFluxFieldFunction()[0])
ffi.SetLogicalVolume(my_lv)
ffi.SetOperationType("avg")
ffi.Initialize()
ffi.Execute()
avg_value = ffi.GetValue()

Available operation types are:

"sum"
"avg"
"max"
"sum_func"
"avg_func"
"max_func"

The *_func variants use a scalar material function supplied with SetOperationFunction.

25.8. Volume Postprocessor

The pyopensn.post.VolumePostprocessor computes scalar-flux integrals, maxima, minima, or volume-weighted averages over spatial regions and energy groups. It produces a single result value per region and group, making it useful for reaction rates and monitoring quantities across geometry subsets or energy ranges.

25.8.1. Basic Usage

Create and execute a volume postprocessor:

from pyopensn.post import VolumePostprocessor

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
)
pps.Execute()
result = pps.GetValue()

Available operation types are:

"integral" — volume-weighted integral of scalar flux
"avg" — volume-weighted average of scalar flux
"max" — maximum scalar flux in region
"min" — minimum scalar flux in region

25.8.2. Spatial Restriction

By default, a postprocessor operates over the entire domain. Restrict computation to mesh blocks:

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    block_ids=[1, 2],
)

Or restrict to one or more logical volumes:

from pyopensn.logvol import RPPLogicalVolume

lv = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=0.5)
pps = VolumePostprocessor(
    problem=phys,
    value_type="avg",
    logical_volumes=[lv],
)

Combine block and logical-volume restrictions — the postprocessor uses only cells inside both:

lv1 = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=0.5)
lv2 = RPPLogicalVolume(xmin=1.0, xmax=2.0, ymin=0.0, ymax=0.5)
pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    logical_volumes=[lv1, lv2],
    block_ids=[1],
)

Each logical volume produces one row of results.

25.8.3. Energy Restriction

By default, a postprocessor returns results for all energy groups. Restrict to a single group:

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    group=6,
)

Or restrict to a single groupset:

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    groupset=1,
)

The group and groupset parameters are mutually exclusive.

25.8.4. Multipliers and Cross-Section Weighting

By default, the postprocessor multiplies each group’s scalar flux by 1.0. Apply a uniform multiplier:

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    multiplier=2.5,
)

Apply group-specific multipliers:

group_mults = [1.0, 1.5, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0]  # one per group
pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    group_multipliers=group_mults,
)

Weight by a cross section (for example, fission rate):

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    xs_multiplier="sigma_f",
)

The cross-section name must exist in the problem’s XS definitions. Only one of multiplier, group_multipliers, and xs_multiplier may be specified.

25.8.5. Results

After calling Execute(), retrieve results with GetValue(). The return value is a 2D array indexed as result[region][group]:

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
)
pps.Execute()
result = pps.GetValue()

# Single region, all groups
# result[0] is a vector of values, one per group
for group_index, value in enumerate(result[0]):
    print(f"Group {group_index}: {value}")

With multiple logical volumes:

lv1 = RPPLogicalVolume(...)
lv2 = RPPLogicalVolume(...)
pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    logical_volumes=[lv1, lv2],
)
pps.Execute()
result = pps.GetValue()

# result[0] is values for lv1
# result[1] is values for lv2

With energy restriction:

pps = VolumePostprocessor(
    problem=phys,
    value_type="integral",
    group=3,
)
pps.Execute()
result = pps.GetValue()

# result[0] is a vector with one element (the single group)
print(result[0][0])

25.9. Cross-Section Sensitivity Postprocessor

The pyopensn.post.CrossSectionSensitivityPostprocessor computes adjoint-weighted sensitivities of a response to a total, scattering, or production (nu_sigma_f) cross-section coefficient, from forward and adjoint angular fluxes or flux moments. It is intended for discrete ordinates problems that have already been solved in a forward and in an adjoint configuration.

25.9.1. Basic Usage

"sigma_t" sensitivities use angular fluxes (psi); "scatter" and "production" sensitivities use flux moments (phi). A single problem object only holds one state at a time, so a typical pattern solves forward, writes that state to disk, solves adjoint, then builds the postprocessor against the current (adjoint) state and the saved forward files:

from pyopensn.post import CrossSectionSensitivityPostprocessor

# ... solve forward, then phys.WriteAngularFluxes("forward_psi")
# ... solve adjoint (phys now holds the adjoint state)

sens = CrossSectionSensitivityPostprocessor(
    problem=phys,
    sensitivity_type="sigma_t",
    forward_angular_fluxes="forward_psi",
)
sens.Execute()
result = sens.GetValue()

25.9.2. Sensitivity Types

"sigma_t" (default) — total cross-section sensitivity. One column per selected group.
"scatter" — scattering-transfer coefficient sensitivity for a single from_group -> to_group pair. One column per selected scattering moment (Legendre order).
"production" — fission production (nu_sigma_f) sensitivity for a single group. Always a single column.

# Scattering sensitivity for group 2 -> group 1, moment 0
scatter_sens = CrossSectionSensitivityPostprocessor(
    problem=phys,
    sensitivity_type="scatter",
    from_group=2,
    to_group=1,
    moment=0,
)

# Production sensitivity for group 0
prod_sens = CrossSectionSensitivityPostprocessor(
    problem=phys,
    sensitivity_type="production",
    group=0,
)

from_group and to_group are required for "scatter". group is required for "production". For "sigma_t", omitting group computes sensitivities for every group; for "scatter", omitting moment (equivalently ell) computes sensitivities for every scattering moment up to the problem’s scattering order. moment and ell are aliases for the same input and cannot both be supplied.

25.9.3. Spatial Restriction

As with pyopensn.post.VolumePostprocessor, block_ids and logical_volumes restrict the cells included in the sensitivity integral, and can be combined:

from pyopensn.logvol import RPPLogicalVolume

lv = RPPLogicalVolume(xmin=0.0, xmax=1.0, ymin=0.0, ymax=0.5)
sens = CrossSectionSensitivityPostprocessor(
    problem=phys,
    sensitivity_type="sigma_t",
    logical_volumes=[lv],
    block_ids=[1],
)

Each logical volume produces one row of results; with no logical volumes, there is a single row over the (optionally block-restricted) domain.

25.9.4. Forward and Adjoint State

By default, both the forward and the adjoint terms are read from the problem’s current state. To compare an explicit forward and adjoint solve, write one of them to disk and point the postprocessor at it:

forward_angular_fluxes / adjoint_angular_fluxes — file prefixes written by WriteAngularFluxes; used for "sigma_t".
forward_flux_moments / adjoint_flux_moments — file prefixes written by WriteFluxMoments; used for "scatter" and "production", and by ApplyKEigenvalueScaling below.
flux_moments_single_file — set True if the flux moments were written to a single file instead of one file per rank.

Leaving a prefix empty (the default) uses the problem’s current state for that term.

25.9.5. Relative Sensitivities

By default, results are absolute sensitivities dR/dx. Set relative=True to scale by the cross-section value itself, x * dR/dx, which is dimensionless and easier to compare across coefficients:

sens = CrossSectionSensitivityPostprocessor(
    problem=phys,
    sensitivity_type="sigma_t",
    relative=True,
)

25.9.6. k-Eigenvalue Scaling

For k-eigenvalue problems, call pyopensn.post.CrossSectionSensitivityPostprocessor.ApplyKEigenvalueScaling() after Execute() to convert the computed sensitivities into first-order k-eigenvalue sensitivities, normalized by the fission-weighted inner product <psi_adj, F psi> built from the same forward/adjoint flux-moment sources configured at construction. As above, this requires an actual forward and adjoint pair of states — leaving both prefixes empty would take the “forward” and “adjoint” terms from the same current state and would not produce a meaningful sensitivity:

# ... solve forward, then phys.WriteFluxMoments("forward_phi")
# ... solve adjoint (phys now holds the adjoint state)

sens = CrossSectionSensitivityPostprocessor(
    problem=phys,
    sensitivity_type="production",
    group=0,
    forward_flux_moments="forward_phi",
)
sens.Execute()
sens.ApplyKEigenvalueScaling(k_eff=1.00123)
result = sens.GetValue()

adjoint_flux_moments is left unset here because phys already holds the adjoint state at this point, so the default (current state) is correct for that term; only the forward term needs to be pointed at the saved file.

25.9.7. Results

After calling Execute(), retrieve results with GetValue(). The return value is a 2D array indexed as result[region][column], where column runs over groups ("sigma_t"), scattering moments ("scatter"), or is a single entry ("production"):

sens.Execute()
result = sens.GetValue()

# Single region, sigma_t sensitivities for every group
for group_index, value in enumerate(result[0]):
    print(f"Group {group_index}: {value}")

25.10. Other Useful Post-Processing Paths

For discrete ordinates problems, other useful output paths include:

GetPsi() for direct access to stored angular-flux arrays
ComputeLeakage() for boundary leakage
WriteAngularFluxes() and ReadAngularFluxes() for angular flux storage and restart-style workflows

These are not field-function exports, but they are often part of the same post-processing workflow.

25.11. Transient Workflows

For transient problems, the usual pattern is:

advance the timestep,
update or create the desired field functions,
export or interpolate them,
repeat as needed.

If the transient workflow needs angular-flux output as well, the problem must have been created with options={"save_angular_flux": True} from the start.

For example:

scalar_ff = phys.GetScalarFluxFieldFunction()[0]

while not solver.Finished():
    solver.Advance()
    scalar_ff.Update()
    ...

This keeps one field-function object and refreshes it from each completed timestep. Creating a new field function after each Advance() is also valid, but is usually unnecessary when the existing object supports Update().

25.12. Practical Guidance

For most workflows:

use GetScalarFluxFieldFunction() for scalar flux,
use CreateFieldFunction() for power or XS-weighted derived outputs,
use GetAngularFieldFunctionList() only when angular information is actually needed,
use ExportMultipleToPVTU() for visualization output,
use point interpolation for spot checks,
use line interpolation for profiles,
use volume interpolation for averages and integrated responses.
call Update() before reusing a field function after the solver state changes, or create a fresh field function from the current state.

If a script is becoming complicated, it is often worth separating the solve and the post-processing logic into different helper functions. That keeps the transport setup readable and makes the output workflow easier to reuse.

25. Post Processors

25.1. Overview

25.2. Updating Existing Field Functions

25.3. Scalar Flux Field Functions

25.3.1. GetScalarFluxFieldFunction

25.4. Derived Field Functions

25.4.1. CreateFieldFunction

25.5. Angular Flux Field Functions

25.5.1. GetAngularFieldFunctionList

25.6. Exporting Field Functions

25.7. Field-Function Interpolation

25.7.1. Point Interpolation

25.7.2. Line Interpolation

25.7.3. Volume Interpolation

25.8. Volume Postprocessor

25.8.1. Basic Usage

25.8.2. Spatial Restriction

25.8.3. Energy Restriction

25.8.4. Multipliers and Cross-Section Weighting

25.8.5. Results

25.9. Cross-Section Sensitivity Postprocessor

25.9.1. Basic Usage

25.9.2. Sensitivity Types

25.9.3. Spatial Restriction

25.9.4. Forward and Adjoint State

25.9.5. Relative Sensitivities

25.9.6. k-Eigenvalue Scaling

25.9.7. Results

25.10. Other Useful Post-Processing Paths

25.11. Transient Workflows

25.12. Practical Guidance

25.3.1. `GetScalarFluxFieldFunction`

25.4.1. `CreateFieldFunction`

25.5.1. `GetAngularFieldFunctionList`