3.6. Operator Construction Methods

In discrete-ordinates transport, we need operators to map between the angular flux (defined at discrete directions) and flux moments (expanded in spherical harmonics). The two key operators are:

  • Discrete-to-Moment (D): computes flux moments from the angular flux.

  • Moment-to-Discrete (M): reconstructs the angular flux from flux moments.

OpenSn provides various methods for constructing these operators:

  1. standard

  2. galerkin_one

  3. galerkin_three

This tutorial demonstrates each method and compares their properties.

To run the code, simply type: jupyter nbconvert --to python --execute <basename>.ipynb.

To convert it to a python file (named <basename>.py), simply type: jupyter nbconvert --to python <basename>.ipynb

[11]:

import os
import sys
import numpy as np
from mpi4py import MPI

sys.path.append("../../../..")

from pyopensn.aquad import GLCProductQuadrature3DXYZ

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb Cell 2 line 8
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#W1sdnNjb2RlLXJlbW90ZQ%3D%3D?line=3'>4</a> from mpi4py import MPI
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#W1sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a> sys.path.append("../..")
----> <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#W1sdnNjb2RlLXJlbW90ZQ%3D%3D?line=7'>8</a> from pyopensn.aquad import GLCProductQuadrature3DXYZ
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#W1sdnNjb2RlLXJlbW90ZQ%3D%3D?line=8'>9</a> from pyopensn.context import UseColor, Finalize
     <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#W1sdnNjb2RlLXJlbW90ZQ%3D%3D?line=10'>11</a> UseColor(False)

ModuleNotFoundError: No module named 'pyopensn.aquad'

3.6.1. Helper function

We define a helper function to evaluate the orthogonality of the D2M and M2D operators. Ideally, the product \(\mathbf{D} \cdot \mathbf{M}\) should equal the identity matrix, which means the operators are consistent inverses of each other. For cases where \(N_{mom} \leq N_{dir}\), \(\mathbf{D}\) will be the left inverse of \(\mathbf{M}\).

OpenSn does not store \(\mathbf{D}\) as it is seen in literature. Instead, \(\mathbf{D}^T\) is stored, and you will see various transposes done to convert OpenSn’s \(\mathbf{D}\) into the literature defined \(\mathbf{D}\).

[ ]:

def check_orthogonality(D2M, M2D, label):
    """Check how close D @ M is to the identity matrix."""
    product = D2M.T @ M2D
    n = min(product.shape)
    diag = np.diag(product)
    mask = np.ones(product.shape, dtype=bool)
    np.fill_diagonal(mask, False)
    max_off_diag = np.abs(product[mask]).max() if product[mask].size > 0 else 0.0
    max_diag_dev = np.abs(diag - 1.0).max()
    print(f"{label}:")
    print(f"  D2M shape: {D2M.shape}, M2D shape: {M2D.shape}")
    print(f"  Max off-diagonal of D2M^T @ M2D: {max_off_diag:.6e}")
    print(f"  Max diagonal deviation from 1:   {max_diag_dev:.6e}")
    print()
    return max_off_diag, max_diag_dev

3.6.2. Standard method

The Standard method builds the operators directly from the spherical harmonics evaluated at the quadrature points:

\[\mathbf{D2M}_{n,m} = w_n \, Y_{\ell}^{m}(\hat{\Omega}_n)\]
\[\mathbf{M2D}_{n,m} = \frac{2\ell + 1}{\sum_n w_n} \, Y_{\ell}^{m}(\hat{\Omega}_n)\]

This is the default method and works with any quadrature set and any scattering order. However, for a finite quadrature, \(\mathbf{D} \cdot \mathbf{M}\) is not guaranteed to be the identity.

[ ]:

pquad_std = GLCProductQuadrature3DXYZ(
    n_polar=4, n_azimuthal=8, scattering_order=3, operator_method='standard'
)

D2M_pquad_std = pquad_std.GetDiscreteToMomentOperator()
M2D_pquad_std = pquad_std.GetMomentToDiscreteOperator()

off_pquad_std, diag_pquad_std = check_orthogonality(D2M_pquad_std, M2D_pquad_std, "Standard")

3.6.3. Galerkin One method

The Galerkin One method builds the M2D operator the same way as the Standard method, then computes D2M by directly inverting the M2D matrix. This requires the operator to be square, i.e., the number of directions must equal the number of moments.

When using galerkin_one, the scattering_order parameter is optional. If omitted, OpenSn automatically selects the scattering order so that the number of spherical harmonic moments equals the number of quadrature directions, yielding a square (and thus invertible) system.

This method produces operators that satisfy \(\mathbf{D} \cdot \mathbf{M} = \mathbf{I}\) to machine precision, as long as \(\mathbf{M}\) is well-conditioned.

[6]:

pquad_g1 = GLCProductQuadrature3DXYZ(
    n_polar=4, n_azimuthal=8, operator_method='galerkin_one'
)

D2M_pquad_g1 = pquad_g1.GetDiscreteToMomentOperator()
M2D_pquad_g1 = pquad_g1.GetMomentToDiscreteOperator()

off_pquad_g1, diag_pquad_g1 = check_orthogonality(D2M_pquad_g1, M2D_pquad_g1, "Galerkin One")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb Cell 8 line 1
----> <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a> pquad_g1 = GLCProductQuadrature3DXYZ(
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=1'>2</a>     n_polar=4, n_azimuthal=8, operator_method='galerkin_one'
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=2'>3</a> )
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a> D2M_g1 = pquad_g1.GetDiscreteToMomentOperator()
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X10sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a> M2D_g1 = pquad_g1.GetMomentToDiscreteOperator()

NameError: name 'GLCProductQuadrature3DXYZ' is not defined

3.6.4. Galerkin Three method

The Galerkin Three method provides a middle ground. It orthogonalizes the spherical harmonics to form a set of approximate spherical harmonics that are orthogonal with respect to the given quadrature rule. Both D2M and M2D are built from these approximate, orthogonalized harmonics.

Unlike galerkin_one, this method does not require a square operator, so it works with any scattering order. Due to the construction of \(\mathbf{D}\) and \(\mathbf{M}\), it guarantees \(\mathbf{D} \cdot \mathbf{M} = \mathbf{I}\).

[ ]:

pquad_g3 = GLCProductQuadrature3DXYZ(
    n_polar=4, n_azimuthal=8, scattering_order=3, operator_method='galerkin_three'
)

D2M_pquad_g3 = pquad_g3.GetDiscreteToMomentOperator()
M2D_pquad_g3 = pquad_g3.GetMomentToDiscreteOperator()

off_pquad_g3, diag_pquad_g3 = check_orthogonality(D2M_pquad_g3, M2D_pquad_g3, "Galerkin Three")

3.6.5. Summary comparison

The table below compares the orthogonality of \(\mathbf{D} \cdot \mathbf{M}\) for each method using a 3D-Product quadrature set with 4 Polar angles and 8 Azimuthal angles. Smaller values indicate operators that are closer to being exact inverses of each other.

[ ]:

print(f"{'Method':<20} {'Max off-diagonal':>20} {'Max diag deviation':>20}")
print("-" * 62)
print(f"{'Standard':<20} {off_pquad_std:>20.6e} {diag_pquad_std:>20.6e}")
print(f"{'Galerkin One':<20} {off_pquad_g1:>20.6e} {diag_pquad_g1:>20.6e}")
print(f"{'Galerkin Three':<20} {off_pquad_g3:>20.6e} {diag_pquad_g3:>20.6e}")

3.6.6. Comparison 2

We will now show how these operators can differ from eachother using a different quadrature set.

First, we apply the Standard method to an LDFE quadrature set with 1 level of global refinement.

[12]:

from pyopensn.aquad import SLDFEsqQuadrature3DXYZ

ldfe_std = SLDFEsqQuadrature3DXYZ(
    level=1, scattering_order=10, operator_method="standard"
)

D2M_ldfe_std = ldfe_std.GetDiscreteToMomentOperator()
M2D_ldfe_std = ldfe_std.GetMomentToDiscreteOperator()

off_ldfe_std, diag_ldfe_std = check_orthogonality(D2M_ldfe_std, M2D_ldfe_std, "Standard")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb Cell 15 line 1
----> <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X20sdnNjb2RlLXJlbW90ZQ%3D%3D?line=0'>1</a> ldfe_std = SLDFEsqQuadrature3DXY(
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X20sdnNjb2RlLXJlbW90ZQ%3D%3D?line=1'>2</a>     level=1, scattering_order=20, operator_method="standard"
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X20sdnNjb2RlLXJlbW90ZQ%3D%3D?line=2'>3</a> )
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X20sdnNjb2RlLXJlbW90ZQ%3D%3D?line=4'>5</a> D2M_ldfe_std = ldfe_std.GetDiscreteToMomentOperator()
      <a href='vscode-notebook-cell://wsl%2Bubuntu/root/opensn/doc/source/tutorials/aquad/operator_methods.ipynb#X20sdnNjb2RlLXJlbW90ZQ%3D%3D?line=5'>6</a> M2D_ldfe_std = ldfe_std.GetMomentToDiscreteOperator()

NameError: name 'SLDFEsqQuadrature3DXY' is not defined

Next, we do the same with the Galerkin One method.

[ ]:

ldfe_gq1 = SLDFEsqQuadrature3DXYZ(
    level=1, operator_method="galerkin_one"
)

D2M_ldfe_gq1 = ldfe_gq1.GetDiscreteToMomentOperator()
M2D_ldfe_gq1 = ldfe_gq1.GetMomentToDiscreteOperator()

off_ldfe_gq1, diag_ldfe_gq1 = check_orthogonality(D2M_ldfe_gq1, M2D_ldfe_gq1, "Galerkin One")

Finally, we apply the Galerkin Three method.

[ ]:

ldfe_gq3 = SLDFEsqQuadrature3DXYZ(
    level=1, scattering_order=10, operator_method="galerkin_three"
)

D2M_ldfe_gq3 = ldfe_gq3.GetDiscreteToMomentOperator()
M2D_ldfe_gq3 = ldfe_gq3.GetMomentToDiscreteOperator()

off_ldfe_gq3, diag_ldfe_gq3 = check_orthogonality(D2M_ldfe_gq3, M2D_ldfe_gq3, "Galerkin Three")

3.6.7. Summary Comparison

The table below compares the orthogonality of \(\mathbf{D} \cdot \mathbf{M}\). It can be seen that the deviation of each matrix varies, depending on the construction chosen.

[ ]:

print(f"{'Method':<20} {'Max off-diagonal':>20} {'Max diag deviation':>20}")
print("-" * 62)
print(f"{'Standard':<20} {off_ldfe_std:>20.6e} {diag_ldfe_std:>20.6e}")
print(f"{'Galerkin One':<20} {off_ldfe_gq1:>20.6e} {diag_ldfe_gq1:>20.6e}")
print(f"{'Galerkin Three':<20} {off_ldfe_gq3:>20.6e} {diag_ldfe_gq3:>20.6e}")

3.6.8. Finalize (for Jupyter Notebook only)

In Python script mode, PyOpenSn automatically handles environment termination. However, this automatic finalization does not occur when running in a Jupyter notebook, so explicit finalization of the environment at the end of the notebook is required. Do not call the finalization in Python script mode, or in console mode.

Note that PyOpenSn’s finalization must be called before MPI’s finalization.

[ ]:

from IPython import get_ipython

def finalize_env():
    Finalize()
    MPI.Finalize()

ipython_instance = get_ipython()
if ipython_instance is not None:
    ipython_instance.events.register("post_execute", finalize_env)