2.6.2. Volumetric Source over a Logical Volume

This tutorial uses a logical volume to place a source independently of the mesh block assignments.

2.6.2.1. Define the source region

Every cell belongs to block 0, but the source is active only from \(z=0.5\) to \(z=1.5\). Passing the RPPLogicalVolume directly to VolumetricSource avoids creating a separate material block for that region.

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_region = RPPLogicalVolume(
    infx=True, infy=True, zmin=0.5, zmax=1.5
)
xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.0)
quadrature = GLProductQuadrature1DSlab(n_polar=16, scattering_order=0)

2.6.2.2. Define the source

The source definition directly associates the strength with the previously constructed logical volume.

source = VolumetricSource(
    logical_volume=source_region, group_strength=[2.0]
)

2.6.2.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], "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.2.4. Verify the logical-volume source

The selected region has unit length and source strength two, so its total production rate is two. Reflecting boundaries require the same total absorption rate.

balance = solver.ComputeBalanceTable()
production = balance["production_rate"]
absorption = balance["absorption_rate"]
balance_error = abs(production - absorption)
if rank == 0:
    print(f"Logical-volume production rate={production:.8e}")
    print(f"Logical-volume absorption rate={absorption:.8e}")
    print(f"Logical-volume balance error={balance_error:.8e}")
assert abs(production - 2.0) < 1.0e-8
assert balance_error < 1.0e-8

2.6.2.5. Replace the source

Source-setting methods add to the existing source collection unless it is cleared explicitly. The following calls remove the original volumetric source, install a replacement over the same logical volume, and solve the updated problem without rebuilding its mesh, materials, or groupset.

replacement_source = VolumetricSource(
    logical_volume=source_region, group_strength=[1.0]
)
problem.SetVolumetricSources(clear_volumetric_sources=True)
problem.SetVolumetricSources(volumetric_sources=[replacement_source])
solver.Execute()

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