1.1.2. Your First 1D Fixed-Source Calculation

This notebook builds a small one-group slab problem from start to finish. It is adapted from the one-dimensional transport regression inputs, with the mesh and material data reduced for an introductory serial calculation.

1.1.2.1. Learning objectives

By the end, you will be able to create an orthogonal mesh, define one-group cross sections and a volumetric source, select a slab angular quadrature, solve a steady-state problem, and verify a scalar-flux value.

Prerequisite: Anatomy of an OpenSn Calculation.

1.1.2.2. Import the OpenSn objects

Each import corresponds to one stage of the input: mesh construction, material data, angular discretization, sources, problem and solver construction, and result interrogation.

if "opensn_console" not in globals():
    from mpi4py import MPI
    from pyopensn.aquad import GLProductQuadrature1DSlab
    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.2.3. Create the slab mesh

A one-dimensional OpenSn mesh is aligned with the \(z\) axis. The 41 node locations below define 40 equal cells on \(0 \le z \le 4\). Every cell receives block ID 0 so that one material can be assigned to the whole slab.

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

1.1.2.4. Define material, source, and angular data

The material is purely absorbing: its total cross section is \(\Sigma_t=1\) and its scattering ratio is \(c=0\). A uniform isotropic source of strength 1 acts throughout block 0. The Gauss-Legendre product quadrature represents direction with eight polar ordinates.

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

1.1.2.5. Assemble and solve the problem

A groupset connects energy group 0 to the angular quadrature and supplies the inner iterative-solver controls. Vacuum boundaries are the default, so no incoming angular flux is prescribed at either end of the slab.

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.2.6. Verify the scalar flux

In an infinite purely absorbing medium, the scalar flux would approach source strength divided by total cross section, or 1 here. Vacuum leakage lowers the solution near the slab ends, while the maximum near the center remains close to 1.

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_1D_MAX_FLUX={max_flux:.8e}")

1.1.2.7. Interpretation and exercises

The center is least affected by vacuum leakage and therefore has the largest scalar flux. Try halving the slab length, increasing the total cross section, or changing the scattering ratio and predict how the maximum will respond before rerunning the notebook.

Next: Your First 2D Fixed-Source Calculation. See also the Python API pages for OrthogonalMeshGenerator, MultiGroupXS, and DiscreteOrdinatesProblem.

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

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