3.2. Product Quadrature

This set uses a Gauss-Legendre quadrature along the polar angle, and a Gauss-Chebyshev quadrature along the azimuthal angle.

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

import os
import sys
import numpy as np

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

from pyopensn.aquad import GLCProductQuadrature2DXY
from pyopensn.context import UseColor, Finalize

UseColor(False)

3.2.1. Quadrature parameters

Here, we use a Product Quadrature for two-dimensional geometries.

We pick a quadrature set with 4 polar angles (whose cosines are between -1 and 1), and 12 azimuthal angles (between 0 and \(2\pi\)). Hence, in 2D, we will have a total of \(2 \times 12 = 24\) directions.

# Create a 2D angular quadrature with 4 polar and 12 azimuthal angles.
pquad = GLCProductQuadrature2DXY(4, 12)

3.2.2. Retrieve directions

vec3_omegas = pquad.omegas

n_directions = len(vec3_omegas)
print('number of directions =',n_directions)
omegas = np.zeros((n_directions,3))

for d in range(n_directions):
    omegas[d,:] = [vec3_omegas[d].x, vec3_omegas[d].y, vec3_omegas[d].z]

dim = 2

3.2.3. Create a function to plot the quadrature

import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
import itertools

colors = itertools.cycle(["r", "g", "b", "k", "m", "c", "y", "crimson"])

# Create figure with options
fig = plt.figure(dpi=150)
ax = fig.add_subplot(111, projection='3d')

# Transparent sphere data
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
# Plot the surface
ax.plot_surface(x, y, z, color='b', alpha=0.1)

# Create a custom 3D arrow class
class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        # Initialize with dummy 2D points
        super().__init__((0, 0), (0, 0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def do_3d_projection(self, renderer=None):
        # Project 3D coordinates to 2D using the axes projection matrix
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.get_proj())
        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
        # Return a numerical depth value (for example, the mean of the projected z-values)
        return np.mean(zs)

    def draw(self, renderer):
        # Perform projection using the axes' projection matrix
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.get_proj())
        self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
        return super().draw(renderer)


start = -1.0

# Create axes arrows
a = Arrow3D(
    [start, 1.1],
    [0, 0],
    [0, 0],
    mutation_scale=10,
    lw=0.5,
    arrowstyle="-|>",
    color="darkorange",
)
ax.add_artist(a)

a = Arrow3D(
    [0, 0],
    [start, 1.1],
    [0, 0],
    mutation_scale=10,
    lw=0.5,
    arrowstyle="-|>",
    color="darkorange",
)
ax.add_artist(a)

a = Arrow3D(
    [0, 0],
    [0, 0],
    [start, 1.1],
    mutation_scale=10,
    lw=0.5,
    arrowstyle="-|>",
    color="darkorange",
)
ax.add_artist(a)

# The following code plots quadrature directions.
# Variables 'dim', 'n_directions', and 'omegas' must be defined appropriately elsewhere.
d_skip = 0
n_octants = int(2 ** dim)
n_dir_per_oct = int(n_directions / n_octants)
print("n octants=", n_octants, "ndir per octant=", n_dir_per_oct)
for oc in range(n_octants):
    clr = next(colors)
    d_skip = oc * n_dir_per_oct
    for d in range(n_dir_per_oct):
        om = omegas[d + d_skip, :]
        ax.plot3D([0, om[0]], [0, om[1]], [0, om[2]], c=clr, linewidth=0.75)

mu = omegas[:, -1]
polar_level = np.unique(mu)
for r in polar_level:
    x = np.sqrt(1 - r ** 2) * np.cos(u)
    y = np.sqrt(1 - r ** 2) * np.sin(u)
    z = r * np.ones_like(u)
    ax.plot3D(x, y, z, 'grey', linestyle="dashed", linewidth=0.5)

ax.view_init(30, 70)

# uncomment this line for interactive plot
# plt.show()

3.2.4. 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)