4.1.3. Cross-Section Sensitivities

This tutorial uses CrossSectionSensitivityPostprocessor to compute first-order sensitivities of a two-group, volume-integrated scalar-flux response to total and isotropic-scattering cross sections. It demonstrates the complete forward/adjoint workflow, including saving the forward solution data needed after the problem is switched to adjoint mode.

4.1.3.1. Define a problem with analytic sensitivities

Consider a homogeneous, two-group slab of length \(L\) with reflecting boundaries. A uniform source \(q_0\) is applied in group 0, which can scatter into group 1. There is no upscatter. The adjacent cross_section_sensitivity_2g.xs file supplies all material coefficients used by the transport solve:

Coefficient	Value
\(\Sigma_{t,0}\)	1.0
\(\Sigma_{t,1}\)	0.8
\(\Sigma_{s,0\rightarrow0}\)	0.2
\(\Sigma_{s,0\rightarrow1}\)	0.3
\(\Sigma_{s,1\rightarrow1}\)	0.1

Using those cross sections, define the within-group removal coefficients

\[a_0=\Sigma_{t,0}-\Sigma_{s,0\rightarrow0}, \qquad a_1=\Sigma_{t,1}-\Sigma_{s,1\rightarrow1}.\]

The uniform group fluxes and the group-1 response are then

\[\phi_0=\frac{q_0}{a_0}, \qquad \phi_1=\frac{\Sigma_{s,0\rightarrow1}\phi_0}{a_1}, \qquad R=\int_0^L\phi_1(z)\,dz =\frac{Lq_0\Sigma_{s,0\rightarrow1}}{a_0a_1}.\]

Consequently, \(\partial R/\partial\Sigma_{t,g}=-R/a_g\), while \(\partial R/\partial\Sigma_{s,0\rightarrow1}=R/\Sigma_{s,0\rightarrow1}\). This retains a transparent analytic answer while making the group-selection arguments meaningful.

import os
from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.post import CrossSectionSensitivityPostprocessor
from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver
from pyopensn.source import VolumetricSource
from pyopensn.xs import MultiGroupXS

rank = MPI.COMM_WORLD.rank
xs = MultiGroupXS()
xs.LoadFromOpenSn("cross_section_sensitivity_2g.xs")

# Mirror the loaded transfer coefficients only for the analytic check below.
# The transport problem itself obtains them from the cross-section file.
sigma_t = xs.sigma_t
self_scatter = [0.2, 0.1]
downscatter = 0.3
removal = [sigma_t[g] - self_scatter[g] for g in range(xs.num_groups)]
source_strength = 1.0
length = 1.0
num_cells = 24

nodes = [length * i / num_cells for i in range(num_cells + 1)]
mesh = OrthogonalMeshGenerator(node_sets=[nodes]).Execute()
mesh.SetUniformBlockID(0)

forward_source = VolumetricSource(
    block_ids=[0], group_strength=[source_strength, 0.0]
)
boundaries = [
    {"name": "zmin", "type": "reflecting"},
    {"name": "zmax", "type": "reflecting"},
]
quadrature = GLProductQuadrature1DSlab(n_polar=32, scattering_order=0)
problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=xs.num_groups,
    groupsets=[
        {
            "groups_from_to": (0, 1),
            "angular_quadrature": quadrature,
            "inner_linear_method": "petsc_gmres",
            "l_abs_tol": 1.0e-11,
            "l_max_its": 200,
        }
    ],
    xs_map=[{"block_ids": [0], "xs": xs}],
    volumetric_sources=[forward_source],
    boundary_conditions=boundaries,
    options={
        "save_angular_flux": True,
        "verbose_inner_iterations": False,
        "verbose_outer_iterations": False,
    },
)
solver = SteadyStateSourceSolver(problem=problem)
solver.Initialize()
solver.Execute()

forward_flux_prefix = "xs_sensitivity_tutorial_forward_phi_"
forward_angular_prefix = "xs_sensitivity_tutorial_forward_psi_"
problem.WriteFluxMoments(forward_flux_prefix)
problem.WriteAngularFluxes(forward_angular_prefix)

4.1.3.2. Solve the adjoint response problem

The response is the group-1 scalar flux integrated over the entire slab, so the adjoint problem uses a uniform unit source in group 1 and zero source in group 0. SetAdjoint(True) clears the previous source, boundary, and flux state; the response source and reflecting boundaries must therefore be applied again before executing the solver.

Only the forward state is written to files. When an adjoint prefix is omitted from a sensitivity postprocessor, it reads the current adjoint state from problem.

response_source = VolumetricSource(block_ids=[0], group_strength=[0.0, 1.0])
problem.SetAdjoint(True)
problem.SetBoundaryOptions(boundary_conditions=boundaries)
problem.SetVolumetricSources(volumetric_sources=[response_source])
solver.Execute()

4.1.3.3. Evaluate and interpret the sensitivities

A sigma_t sensitivity uses the forward and adjoint angular fluxes because total-cross-section perturbations act on the angular transport operator. Omitting group returns one column per energy group. A scatter sensitivity uses flux moments; moment=0, from_group=0, and to_group=1 select the isotropic downscattering coefficient. Both postprocessors return a two-dimensional array indexed first by spatial region and then by the selected group or scattering moment.

These calls compute absolute derivatives \(\partial R/\partial x\). Set relative=True to obtain \(x\,\partial R/\partial x\) instead. The optional block_ids and logical_volumes arguments can restrict a result to the contribution from selected parts of a heterogeneous model. For a k-eigenvalue problem, sensitivity_type="production" selects a fission-production coefficient; call ApplyKEigenvalueScaling(k_eff) after Execute() to convert the raw bilinear form to a first-order eigenvalue sensitivity.

total_sensitivity = CrossSectionSensitivityPostprocessor(
    problem=problem,
    sensitivity_type="sigma_t",
    forward_angular_fluxes=forward_angular_prefix,
)
total_sensitivity.Execute()
d_response_d_sigma_t = [float(value) for value in total_sensitivity.GetValue()[0]]

scatter_sensitivity = CrossSectionSensitivityPostprocessor(
    problem=problem,
    sensitivity_type="scatter",
    moment=0,
    from_group=0,
    to_group=1,
    forward_flux_moments=forward_flux_prefix,
)
scatter_sensitivity.Execute()
d_response_d_downscatter = float(scatter_sensitivity.GetValue()[0][0])

response = length * source_strength * downscatter / (removal[0] * removal[1])
expected_total_sensitivity = [-response / removal[g] for g in range(2)]
expected_scatter_sensitivity = response / downscatter
sensitivity_errors = [
    abs(computed - expected)
    for computed, expected in zip(
        d_response_d_sigma_t, expected_total_sensitivity
    )
]
sensitivity_errors.append(
    abs(d_response_d_downscatter - expected_scatter_sensitivity)
)
max_sensitivity_error = max(sensitivity_errors)
if rank == 0:
    print(f"Group-0 sigma-t sensitivity={d_response_d_sigma_t[0]:.8e}")
    print(f"Group-1 sigma-t sensitivity={d_response_d_sigma_t[1]:.8e}")
    print(f"Group-0-to-1 scatter sensitivity={d_response_d_downscatter:.8e}")
    print(f"Cross-section sensitivity max error={max_sensitivity_error:.8e}")

MPI.COMM_WORLD.Barrier()
for prefix in (forward_flux_prefix, forward_angular_prefix):
    try:
        os.remove(f"{prefix}{rank}.h5")
    except FileNotFoundError:
        pass

assert max_sensitivity_error < 1.0e-8

if "opensn_console" not in globals():
    from IPython import get_ipython
    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()