2.6.3. Point Source

This tutorial places an isotropic point source inside a two-dimensional transport domain.

2.6.3.1. Build the transport domain

PointSource accepts a physical location and one strength per energy group. The source is placed away from mesh lines so that it belongs unambiguously to one cell. Unlike a volumetric source, its total strength does not scale with cell volume.

from mpi4py import MPI
from pyopensn.aquad import GLCProductQuadrature2DXY
from pyopensn.context import Finalize
from pyopensn.mesh import OrthogonalMeshGenerator
from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSourceSolver
from pyopensn.source import PointSource
from pyopensn.xs import MultiGroupXS

rank = MPI.COMM_WORLD.rank
nodes = [-1.0 + i / 10.0 for i in range(21)]
mesh = OrthogonalMeshGenerator(node_sets=[nodes, nodes]).Execute()
mesh.SetOrthogonalBoundaries()
mesh.SetUniformBlockID(0)

xs = MultiGroupXS()
xs.CreateSimpleOneGroup(sigma_t=1.0, c=0.0)
quadrature = GLCProductQuadrature2DXY(
    n_polar=2, n_azimuthal=8, scattering_order=0
)

2.6.3.2. Define the source

The point source is fully specified by its physical location and one strength for each energy group.

source = PointSource(
    location=[0.037, -0.041, 0.0], strength=[1.0]
)

2.6.3.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}],
    point_sources=[source],
    boundary_conditions=[
        {"name": "xmin", "type": "reflecting"},
        {"name": "xmax", "type": "reflecting"},
        {"name": "ymin", "type": "reflecting"},
        {"name": "ymax", "type": "reflecting"},
    ],
)
solver = SteadyStateSourceSolver(problem=problem, compute_balance=True)
solver.Initialize()
solver.Execute()

2.6.3.4. Verify the point-source strength

The point source has total strength one. With reflecting boundaries and a purely absorbing material, all source particles must eventually contribute to absorption.

balance = solver.ComputeBalanceTable()
production = balance["production_rate"]
absorption = balance["absorption_rate"]
balance_error = abs(production - absorption)
if rank == 0:
    print(f"Point-source production rate={production:.8e}")
    print(f"Point-source absorption rate={absorption:.8e}")
    print(f"Point-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()