3.1.5. A Nested-Square Fixed-Source Problem

This two-dimensional problem demonstrates transport through three nested material regions. A unit isotropic source fills the lower-left square, a highly scattering shell surrounds it, and a more absorbing background fills the remainder of the domain. Particles are removed by absorption or escape through the vacuum exterior boundaries.

3.1.5.1. Import the OpenSn objects

[ ]:

if "opensn_console" not in globals():
    from mpi4py import MPI
    from pyopensn.aquad import GLCProductQuadrature2DXY
    from pyopensn.context import Finalize, UseColor
    from pyopensn.logvol import RPPLogicalVolume
    from pyopensn.mesh import OrthogonalMeshGenerator
    from pyopensn.post import VolumePostprocessor
    from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver
    from pyopensn.source import VolumetricSource
    from pyopensn.xs import MultiGroupXS

    rank = MPI.COMM_WORLD.rank
    UseColor(False)

3.1.5.2. Create the mesh

The domain is the square \([0,10]\times[0,10]\) cm. Ten cells per centimeter place mesh nodes on the material interfaces at 2 and 4 cm. As with the other two-dimensional tutorials, the calculation represents a unit thickness in the \(z\) direction.

[ ]:

domain_length = 10.0
cells_per_cm = 10
num_cells = int(domain_length * cells_per_cm)
nodes = [i / cells_per_cm for i in range(num_cells + 1)]

mesh = OrthogonalMeshGenerator(node_sets=[nodes, nodes]).Execute()
mesh.SetOrthogonalBoundaries()

3.1.5.3. Assign the nested material regions

Region I occupies \([0,2]\times[0,2]\) cm. Region II is the part of \([0,4]\times[0,4]\) cm outside Region I, and Region III fills the rest of the domain.

Nested-square material layout

The mesh is initialized as Region III. Region II is then assigned over the 4 cm square, followed by Region I over the 2 cm square. The final assignment overwrites the cells in the source region.

[ ]:

REGION_I = 0
REGION_II = 1
REGION_III = 2

mesh.SetUniformBlockID(REGION_III)

region_ii_extent = RPPLogicalVolume(
    xmin=0.0, xmax=4.0,
    ymin=0.0, ymax=4.0,
    infz=True,
)
mesh.SetBlockIDFromLogicalVolume(region_ii_extent, REGION_II, True)

region_i_extent = RPPLogicalVolume(
    xmin=0.0, xmax=2.0,
    ymin=0.0, ymax=2.0,
    infz=True,
)
mesh.SetBlockIDFromLogicalVolume(region_i_extent, REGION_I, True)

3.1.5.4. Define the materials and source

For each material, CreateSimpleOneGroup receives the total cross section \(\sigma_t\) and scattering ratio \(c=\sigma_s/\sigma_t\).

Region	\(\sigma_t\) (cm\(^{-1}\))	\(c\)	\(\sigma_s\) (cm\(^{-1}\))	Source (cm\(^{-3}\) s\(^{-1}\))
I	1.0	0.90	0.90	1.0
II	1.5	0.96	1.44	0.0
III	1.0	0.30	0.30	0.0

Only Region I is selected by the volumetric source.

[ ]:

material_data = [
    (REGION_I, 1.0, 0.90),
    (REGION_II, 1.5, 0.96),
    (REGION_III, 1.0, 0.30),
]

xs_map = []
for block_id, sigma_t, scattering_ratio in material_data:
    xs = MultiGroupXS()
    xs.CreateSimpleOneGroup(sigma_t=sigma_t, c=scattering_ratio)
    xs_map.append({"block_ids": [block_id], "xs": xs})

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

3.1.5.5. Assemble and solve the transport problem

A Gauss–Legendre–Chebyshev quadrature represents particle directions in the \(xy\) plane. All four exterior faces are vacuum boundaries.

[ ]:

quadrature = GLCProductQuadrature2DXY(
    n_polar=2, n_azimuthal=32, scattering_order=0
)

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-9,
            "l_max_its": 300,
            "gmres_restart_interval": 30,
        }
    ],
    xs_map=xs_map,
    volumetric_sources=[source],
    boundary_conditions=[
        {"name": "xmin", "type": "vacuum"},
        {"name": "xmax", "type": "vacuum"},
        {"name": "ymin", "type": "vacuum"},
        {"name": "ymax", "type": "vacuum"},
    ],
)

solver = SteadyStateSourceSolver(problem=problem, compute_balance=True)
solver.Initialize()
solver.Execute()

3.1.5.6. Check the source and reaction-rate balance

The source occupies a \(2\times2\) cm square and has unit strength, so its integrated rate is 4 particles per second per unit depth. The region-wise absorption rates are

\[R_{a,i}=\int_{V_i}\sigma_{a,i}\phi(\mathbf{r})\,dV.\]

Their sum and the outward leakage should reproduce the source rate. The source-region average flux provides a compact regression quantity for the spatial solution.

[ ]:

absorption_rates = {}
for block_id, _, _ in material_data:
    postprocessor = VolumePostprocessor(
        problem=problem,
        value_type="integral",
        block_ids=[block_id],
        xs_multiplier="sigma_a",
    )
    postprocessor.Execute()
    absorption_rates[block_id] = postprocessor.GetValue()[0][0]

source_flux_postprocessor = VolumePostprocessor(
    problem=problem, value_type="avg", block_ids=[REGION_I]
)
source_flux_postprocessor.Execute()
source_region_average_flux = source_flux_postprocessor.GetValue()[0][0]

balance = solver.ComputeBalanceTable()
source_rate = balance["production_rate"] + balance["inflow_rate"]
loss_rate = balance["absorption_rate"] + balance["outflow_rate"]
relative_balance_error = abs(source_rate - loss_rate) / source_rate
total_absorption_rate = sum(absorption_rates.values())

if rank == 0:
    print(f"NESTED_SQUARE_SOURCE_RATE={source_rate:.8e}")
    print(f"NESTED_SQUARE_SOURCE_AVG_FLUX={source_region_average_flux:.8e}")
    print(f"NESTED_SQUARE_ABSORPTION_RATE={total_absorption_rate:.8e}")
    print(f"NESTED_SQUARE_BALANCE_ERROR={relative_balance_error:.8e}")

assert abs(source_rate - 4.0) < 1.0e-8
assert abs(total_absorption_rate - balance["absorption_rate"]) < 1.0e-10
assert relative_balance_error < 1.0e-8

3.1.5.7. Export and visualize the scalar flux

After running the generated Python input, use the following code to export the scalar flux for visualization. It is kept in Markdown so the regression test does not create output files.

from pyopensn.fieldfunc import FieldFunctionGridBased

scalar_flux = problem.GetScalarFluxFieldFunction()[0]
FieldFunctionGridBased.ExportMultipleToPVTU(
    [scalar_flux], "Flux/NestedSquare_Phi"
)

Nested square source flux

3.1.5.8. Next steps

Compare the absorption rate in each material, refine the mesh near the reentrant corners, or increase the azimuthal resolution to examine angular discretization effects.

[ ]:

if "opensn_console" not in globals():
    from IPython import get_ipython

    if get_ipython() is not None:
        Finalize()
        MPI.Finalize()