Symbolics

PyWENO contains a symbolic module to help authors develop their own WENO methods and/or explore the basics of WENO methods. Below are a few quick examples demonstrating how the symbolic routines of PyWENO are used.

Interpolating polynomials

First, let’s build some grid points and y-values:

>>> import sympy
>>> import pyweno
>>> (x0, x1, x2) = sympy.var('x0 x1 x2')
>>> (y0, y1, y2) = sympy.var('y0 y1 y2')

Then, the polynomial that interpolates the points (x0, y0), (x1, y1), and (x2, y2) is given by:

>>> p = pyweno.symbolic.polynomial_interpolator([x0, x1, x2], [y0, y1, y2])
>>> p
y0*(x - x1)*(x - x2)/((x0 - x1)*(x0 - x2)) + y1*(x - x0)*(x - x2)/((x1 - x0)*(x1 - x2)) + y2*(x - x0)*(x - x1)/((x2 - x0)*(x2 - x1))

and is a function of the SymPy variable x. For example:

>>> x = sympy.var('x')
>>> p.subs(x, x2)
y2

The polynomial that interpolates the primitive function of \(f\) such that

\[f(x_i) = \sum_j y_j (x_{j+1} - x_{j})\]

is given by:

>>> P = pyweno.symbolic.primitive_polynomial_interpolator([x0, x1, x2], [y1, y2])
>>> P
y1*(x - x0)*(x - x2)/(x1 - x2) + (x - x0)*(x - x1)*(y1*(x1 - x0) + y2*(x2 - x1))/((x2 - x0)*(x2 - x1))

and is also a function of the SymPy variable x. For example:

>>> P.subs(x, x1)
y1*(x1 - x0)

For uniform grids, one could define the grid points by:

>>> (x, dx) = sympy.var('x dx')
>>> xs = [ dx, 2*dx, 3*dx ]
>>> p = pyweno.symbolic.polynomial_interpolator(xs, [y0, y1, y2])
>>> p
y0*(x - 3*dx)*(x - 2*dx)/(2*dx**2) + y2*(x - dx)*(x - 2*dx)/(2*dx**2) - y1*(x - dx)*(x - 3*dx)/dx**2

Reconstruction coefficients

Hereafter we assume that the grid is uniform. Furthermore, to specify a point within a cell, the interval [-1, 1] is used as a reference.

The reconstruction coefficients for a 5th (=2k-1 where k=3) order WENO scheme corresponding to the reconstruction point at the left side (\(\xi\) = -1) of each grid cell are given by:

>>> c = pyweno.symbolic.reconstruction_coefficients(k=3, xi=[ -1 ])
>>> c
{'k': 3,
 'n': 1,
 (0, 0, 0): 11/6,
 (0, 0, 1): -7/6,
 (0, 0, 2): 1/3,
 (0, 1, 0): 1/3,
 (0, 1, 1): 5/6,
 (0, 1, 2): -1/6,
 (0, 2, 0): -1/6,
 (0, 2, 1): 5/6,
 (0, 2, 2): 1/3

Note that the return value c is a dictionary of SymPy objects, indexed according to c[l,r,j] where l is the index of the reconstruction point and r is the left-shift of the stencil.

Recall that the reconstruction coefficients \(c\) are used to reconstruct the original (unknown) function \(f\) at each point \(\xi_l\) in xi according to

\[f^r(\xi^l) \approx \sum_{j=0}^{k-1} c^l_{r,j} \, \bar{f}_{i-r+j}\]

for each \(l\) from 0 to len(xi), where \(\bar{f}_{i-r+j}\) is the cell average of \(f\) in the cell \(i-r+j\).

Optimal weights

The optimal weights for a 5th (=2k-1 where k=3) order WENO scheme corresponding to the reconstruction point at the left side of each grid cell are given by:

>>> w = pyweno.symbolic.optimal_weights(3, [ -1 ])
>>> w
({'k': 3, 'n': 1, (0, 0): 1/10, (0, 1): 3/5, (0, 2): 3/10}, {0: False, 'n': 1})

Note that the return value w is a tuple of dictionaries of SymPy objects. The first dictionary contains the weights, and is indexed according to w[l,r]. the second dictionary contains boolean values determining if the weights are split (negative).

Recall that the optimal weights are used to obtain an optimally high-order reconstruction of the original function \(f\) given the low-order reconstructions \(f^r(\xi^l)\) according to

\[f(\xi^l) \approx \sum_{r=0}^{k-1} \varpi^{l,r} f^r(\xi_l).\]

Smoothness coefficients

The Jiang-Shu smoothness coefficients for a 5th (=2k-1 where k=3) order WENO scheme are given by:

>>> beta = pyweno.symbolic.jiang_shu_smoothness_coefficients(3)

The return value beta is a dictionary of SymPy objects, and is indexed according to beta[r,m,n] (see the reference documentation for details).

Recall that the smoothness coefficients beta[r, m, n] are used to compute the non-linear weights \(\omega^{l,r}\) (used in place of \(\varpi^{l,r}\) in non-smooth regions) according to

\[\omega^{l,r} = \frac{\alpha^{l,r}}{\alpha^{l,0} + \cdots + \alpha_{l,k-1}}\]

where

\[\alpha^{l,r} = \frac{\varpi^{l,r}}{(\epsilon + \sigma^r)^p}\]

and

\[\sigma^r = \sum_{m=1}^{2k-1} \sum_{n=1}^{2k-1} \beta_{r,m,n}\, \overline{f}_{i-k+m}\, \overline{f}_{i-k+n}.\]