1.1.3. Your First 2D Fixed-Source Calculation

This notebook extends the one-dimensional slab calculation to a square Cartesian domain.

1.1.3.1. Learning objectives

You will create a two-dimensional mesh, name its orthogonal boundaries, choose a quadrature appropriate to 2D Cartesian transport, solve a one-group fixed-source problem, and verify the scalar flux.

Prerequisite: Your First 1D Fixed-Source Calculation.

1.1.3.2. 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.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
    UseColor(False)

1.1.3.3. Create the square mesh

Supplying two node sets creates a Cartesian mesh in the \(x\)-\(y\) plane. SetOrthogonalBoundaries assigns the names xmin, xmax, ymin, and ymax; the problem uses vacuum conditions on all four by default.

num_cells = 16
length = 4.0
nodes = [length * i / num_cells for i in range(num_cells + 1)]
mesh = OrthogonalMeshGenerator(node_sets=[nodes, nodes]).Execute()
mesh.SetOrthogonalBoundaries()
mesh.SetUniformBlockID(0)

1.1.3.4. Define the physical and angular model

As in the 1D example, the material is purely absorbing and a uniform isotropic source fills the domain. The two-dimensional product quadrature combines polar and azimuthal directions so particles can travel throughout the \(x\)-\(y\) plane.

cross_sections = MultiGroupXS()
cross_sections.CreateSimpleOneGroup(sigma_t=1.0, c=0.0)
source = VolumetricSource(block_ids=[0], group_strength=[1.0])
quadrature = GLCProductQuadrature2DXY(
    n_polar=2, n_azimuthal=8, scattering_order=0
)

1.1.3.5. Assemble and solve the problem

The same problem-and-solver pattern applies in any dimension. Only the mesh and angular quadrature have changed from the slab example.

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,
            "gmres_restart_interval": 30,
        }
    ],
    xs_map=[{"block_ids": [0], "xs": cross_sections}],
    volumetric_sources=[source],
)

solver = SteadyStateSourceSolver(problem=problem)
solver.Initialize()
solver.Execute()

1.1.3.6. Verify and interpret the scalar flux

The scalar flux is symmetric about both centerlines because the square, material, source, and boundary conditions are symmetric. Its maximum is near the center and below the infinite-medium value of 1 because particles can leak through four vacuum boundaries.

whole_domain = RPPLogicalVolume(infx=True, infy=True, infz=True)
interpolator = FieldFunctionInterpolationVolume()
interpolator.SetOperationType("max")
interpolator.SetLogicalVolume(whole_domain)
interpolator.AddFieldFunction(problem.GetScalarFluxFieldFunction()[0])
interpolator.Execute()
max_flux = interpolator.GetValue()

if rank == 0:
    print(f"FOUNDATION_2D_MAX_FLUX={max_flux:.8e}")

1.1.3.7. Exercises and next steps

Try using a rectangular domain, increasing n_azimuthal, or making one boundary reflecting. Before rerunning, predict whether the maximum scalar flux should rise or fall.

Next: Running in Serial and Parallel. The modeling tutorials cover mesh, logical-volume, material, source, boundary, and groupset choices in more detail.

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

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