4.5.1. Running a Source-Strength Study

This tutorial varies a source strength directly in Python and records a detector-region response without launching separate terminal commands.

4.5.1.1. Initialize the reusable model

The mesh, materials, groupset, and solver remain fixed. Only the volumetric source is replaced between solves, which keeps the study self-contained and avoids environment variables or generated input files.

from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
from pyopensn.fieldfunc import FieldFunctionInterpolationVolume
from pyopensn.logvol import RPPLogicalVolume
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver
from pyopensn.source import VolumetricSource
from pyopensn.xs import MultiGroupXS

rank = MPI.COMM_WORLD.rank
nodes = [6.0 * i / 120.0 for i in range(121)]
mesh = OrthogonalMeshGenerator(node_sets=[nodes]).Execute()
mesh.SetUniformBlockID(0)
source_region = RPPLogicalVolume(infx=True, infy=True, zmin=1.0, zmax=2.0)
detector_region = RPPLogicalVolume(infx=True, infy=True, zmin=4.0, zmax=5.0)

xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.5)
quadrature = GLProductQuadrature1DSlab(n_polar=32, scattering_order=0)
initial_source = VolumetricSource(logical_volume=source_region, group_strength=[0.5])
problem = DiscreteOrdinatesProblem(
    mesh=mesh,
    num_groups=1,
    groupsets=[
        {
            "groups_from_to": (0, 0),
            "angular_quadrature": quadrature,
            "inner_linear_method": "petsc_gmres",
            "l_abs_tol": 1.0e-10,
        }
    ],
    xs_map=[{"block_ids": [0], "xs": xs}],
    volumetric_sources=[initial_source],
)
solver = SteadyStateSourceSolver(problem=problem)
solver.Initialize()

4.5.1.2. Sweep the source strength

Before each solve, the previous volumetric source is cleared and replaced. Because the transport equation is linear, the detector response divided by source strength should remain constant.

strengths = [0.5, 1.0, 2.0]
responses = []
for strength in strengths:
    problem.SetVolumetricSources(clear_volumetric_sources=True)
    source = VolumetricSource(logical_volume=source_region, group_strength=[strength])
    problem.SetVolumetricSources(volumetric_sources=[source])
    solver.Execute()

    scalar_flux = problem.GetScalarFluxFieldFunction(only_scalar_flux=True)[0]
    detector = FieldFunctionInterpolationVolume()
    detector.SetOperationType("sum")
    detector.SetLogicalVolume(detector_region)
    detector.AddFieldFunction(scalar_flux)
    detector.Execute()
    responses.append(float(detector.GetValue()))

normalized = [response / strength for response, strength in zip(responses, strengths)]
linearity_error = max(normalized) - min(normalized)
if rank == 0:
    for strength, response in zip(strengths, responses):
        print(f"Source strength {strength:.1f}: response={response:.6e}")
    print(f"Source-study linearity error={linearity_error:.6e}")
assert linearity_error < 1.0e-8

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