2.6.1. Material Source Selected by Block ID

This tutorial applies a volumetric source to every cell carrying a selected mesh block ID. This form is useful when a source occupies an entire material region.

2.6.1.1. Build the mesh and source block

The slab is two centimeters long. Cells in its left half are assigned block ID 1, while the remaining cells retain block ID 0. VolumetricSource(block_ids=[1], ...) then activates the source only in the block-1 cells. Both blocks use the same material so the example isolates source selection from material effects.

from mpi4py import MPI
from pyopensn.aquad import GLProductQuadrature1DSlab
from pyopensn.context import Finalize
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 = [i / 20.0 for i in range(41)]
mesh = OrthogonalMeshGenerator(node_sets=[nodes]).Execute()
mesh.SetUniformBlockID(0)

source_block = RPPLogicalVolume(infx=True, infy=True, zmin=0.0, zmax=1.0)
mesh.SetBlockIDFromLogicalVolume(source_block, 1, True)

xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.0)
quadrature = GLProductQuadrature1DSlab(n_polar=16, scattering_order=0)

2.6.1.2. Define the source

The complete source definition is the block selection and one strength value for each energy group.

source = VolumetricSource(block_ids=[1], group_strength=[1.0])

2.6.1.3. Configure and solve the problem

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,
            "l_max_its": 100,
        }
    ],
    xs_map=[{"block_ids": [0, 1], "xs": xs}],
    volumetric_sources=[source],
    boundary_conditions=[
        {"name": "zmin", "type": "reflecting"},
        {"name": "zmax", "type": "reflecting"},
    ],
)
solver = SteadyStateSourceSolver(problem=problem, compute_balance=True)
solver.Initialize()
solver.Execute()

2.6.1.4. Verify the selected source region

The active block has unit length and a source strength of one, giving a total production rate of one. With reflecting boundaries and pure absorption, the absorption rate must also be one.

balance = solver.ComputeBalanceTable()
production = balance["production_rate"]
absorption = balance["absorption_rate"]
balance_error = abs(production - absorption)
if rank == 0:
    print(f"Block-source production rate={production:.8e}")
    print(f"Block-source absorption rate={absorption:.8e}")
    print(f"Block-source balance error={balance_error:.8e}")
assert abs(production - 1.0) < 1.0e-8
assert balance_error < 1.0e-8
if "opensn_console" not in globals():
    from IPython import get_ipython
    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()