diff --git a/development/poly_reconstruction/README.md b/development/poly_reconstruction/README.md new file mode 100644 index 00000000..9bb5a64a --- /dev/null +++ b/development/poly_reconstruction/README.md @@ -0,0 +1,163 @@ +# Polynomial reconstructions in SharpClaw + +SharpClaw offers three linear polynomial reconstructions alongside WENO, selected with `solver.lim_type` (4, 5 or 6). They were written for electromagnetic wave problems, where an option that adds no numerical dissipation is valuable. This note explains what they compute, which cells they use, how accurate they are and how they damp and disperse waves. Everything shown here is regenerated by the scripts in this directory. + +| `lim_type` | name | cells per face value | formal order | numerical dissipation | +|---|---|---|---|---| +| 4 | `'poly'` | N, centred on the face | N | none | +| 5 | `'poly_dissipative'` | N, shifted one cell upwind | N | strongest | +| 6 | `'poly_cell_centered'` | N+1, centred on the cell | N+1 | moderate | + +N is `solver.poly_order`, an even number in {4, 6, 8, 10}. + +## 1. What a reconstruction computes + +PyClaw is a finite-volume code: `q[i]` holds the **average** of the solution over cell $i$, + +$$\bar q_i = \frac{1}{h}\int_{x_{i-1/2}}^{x_{i+1/2}} q(x)\,dx .$$ + +At every step SharpClaw needs the solution **at the faces**. For each cell the reconstruction returns two values: `ql(i)`, the value at its left face $x_{i-1/2}$, and `qr(i)`, the value at its right face $x_{i+1/2}$. Face $x_{i+1/2}$ therefore receives two values, `qr(i)` from the left and `ql(i+1)` from the right, and the Riemann solver resolves the jump between them. That jump is where the upwinding, and so the numerical dissipation, comes from. Nothing computes values at cell centres. + +All three reconstructions are linear combinations of nearby averages with fixed weights: + +$$\texttt{ql}(i) = \sum_j w_j\,\bar q_{i+j}, \qquad \texttt{qr}(i) = \sum_j w_j\,\bar q_{i-j}.$$ + +## 2. Why the averages must be treated as averages + +An average is not the value at the cell centre. For smooth $q$, $\bar q_i = q(x_i) + \tfrac{h^2}{24}q''(x_i) + O(h^4)$. Weights that interpolate point values (Lagrange interpolation through the centres) inherit that $O(h^2)$ error, whatever the stencil width. The first version of `poly.f90` did exactly this, and every order converged at second order: + +![Cell averages vs centre values, and old vs new weights](figures/averages_vs_points.png) + +## 3. How the weights are built + +The primitive function $Q(x)=\int_{x_0}^{x} q\,dx$ is known **exactly** at every face: it is a running sum of $h\,\bar q$. Interpolate $Q$ at the faces of the stencil with a polynomial $P$, and the face value is the slope $P'$ at that face. The weights do not depend on $q$, so this is done once, in exact rational arithmetic, by `poly_weights.py`: + +![Primitive-function construction](figures/primitive_construction.png) + +## 4. Which cells each type uses + +With $k = N/2$, the cells used for `ql(i)` are + +| type | cells for `ql(i)` | cells for `qr(i)` | idea | +|---|---|---|---| +| `'poly'` | $i-k \dots i+k-1$ | $i-k+1 \dots i+k$ | N cells centred on the face | +| `'poly_dissipative'` | $i-k+1 \dots i+k$ | $i-k \dots i+k-1$ | the `'poly'` stencil shifted one cell upwind | +| `'poly_cell_centered'` | $i-k \dots i+k$ | $i-k \dots i+k$ | N+1 cells centred on cell $i$ | + +![Stencils for N = 4 and N = 6](figures/stencils.png) + +The difference that matters is what the two sides of one face see. With `'poly'`, `qr(i)` and `ql(i+1)` use the same cells with the same weights, so they are the same number: the Riemann solver sees no jump and adds no dissipation. The other two types reconstruct each side from a different stencil, so the states differ by $O(h^{N})$ (or $O(h^{N+1})$) and the jump is damped. + +![The two states at one face](figures/face_states.png) + +Every type reaches $k$ cells from the cell being reconstructed, so SharpClaw sets `num_ghost = N/2 + 1` automatically. + +### Weights + +Generated by `python poly_weights.py --table`; `qr(i)` uses the mirror image ($w_j$ applied to $\bar q_{i-j}$). + +**lim_type = 4, `'poly'`** + +| N | order | cells for `ql(i)` | weights | +|---|---|---|---| +| 4 | 4 | i-2 … i+1 | (-1, 7, 7, -1) / 12 | +| 6 | 6 | i-3 … i+2 | (1, -8, 37, 37, -8, 1) / 60 | +| 8 | 8 | i-4 … i+3 | (-3, 29, -139, 533, 533, -139, 29, -3) / 840 | +| 10 | 10 | i-5 … i+4 | (2, -23, 127, -473, 1627, 1627, -473, 127, -23, 2) / 2520 | + +**lim_type = 5, `'poly_dissipative'`** + +| N | order | cells for `ql(i)` | weights | +|---|---|---|---| +| 4 | 4 | i-1 … i+2 | (3, 13, -5, 1) / 12 | +| 6 | 6 | i-2 … i+3 | (-2, 22, 57, -23, 7, -1) / 60 | +| 8 | 8 | i-3 … i+4 | (5, -55, 365, 743, -307, 113, -27, 3) / 840 | +| 10 | 10 | i-4 … i+5 | (-3, 37, -233, 1207, 2131, -893, 367, -113, 22, -2) / 2520 | + +**lim_type = 6, `'poly_cell_centered'`** + +| N | order | cells for `ql(i)` | weights | +|---|---|---|---| +| 4 | 5 | i-2 … i+2 | (-3, 27, 47, -13, 2) / 60 | +| 6 | 7 | i-3 … i+3 | (4, -38, 214, 319, -101, 25, -3) / 420 | +| 8 | 9 | i-4 … i+4 | (-5, 55, -305, 1375, 1879, -641, 199, -41, 4) / 2520 | +| 10 | 11 | i-5 … i+5 | (12, -153, 947, -4003, 15797, 20417, -7303, 2597, -703, 122, -10) / 27720 | + +Two checks on these numbers: `'poly'` with N = 4 is the classical fourth-order face value (-1, 7, 7, -1)/12, and `'poly_cell_centered'` with N = 4 gives exactly the ideal (linear) weights of WENO5. + +## 5. Accuracy + +`sine_reconstruction.py` reconstructs $\sin 2\pi x$ from its exact cell averages and compares the face values with the exact ones. At a coarse 4 cells per wavelength the difference between the types is visible: all three are close to the curve, but only `'poly'` gives each face a single value. + +![Face values of sin 4πx from 8 cells](figures/sine_reconstruction.png) + +Errors fall at the formal order until they reach round-off: + +![Convergence of the interface error](figures/convergence.png) + +Measured with 16 → 32 cells (`python sine_reconstruction.py`): + +| type | N = 4 | N = 6 | N = 8 | N = 10 | +|---|---|---|---|---| +| `'poly'` | 3.98 | 5.97 | 7.95 | 9.94 | +| `'poly_dissipative'` | 3.98 | 5.97 | 7.95 | 9.94 | +| `'poly_cell_centered'` | 4.96 | 6.94 | 8.93 | 10.92 | + +In the full solver the observed order is also limited by the time integrator: the default `'SSP104'` is fourth order, so for N > 4 the time error dominates unless the CFL number is reduced. + +## 6. Dissipation and dispersion + +For the model problem $u_t + u_x = 0$ with the upwind flux, the semi-discrete scheme evolves a Fourier mode $e^{ij\theta}$ ($\theta = kh$) as $e^{\lambda t}$; the exact value is $\lambda = -i\theta/h$. $\mathrm{Re}\,\lambda$ is the damping rate and $-\mathrm{Im}\,\lambda$ the modified wavenumber (below the diagonal, waves travel too slowly). + +![Damping rate and modified wavenumber](figures/fourier.png) + +The same data at a few resolutions, as phase-speed ratio $-\mathrm{Im}\,\lambda/\theta$ (1 is exact) and damping rate $\mathrm{Re}\,\lambda$ (with $h = 1$): + +| points per wavelength | `'poly'` N=4 | `'poly'` N=8 | `'poly_dissipative'` N=4 | `'poly_dissipative'` N=8 | `'poly_cell_centered'` N=4 | `'poly_cell_centered'` N=8 | +|---|---|---|---|---|---|---| +| 16 | 0.9992 / 0 | 1.0000 / 0 | 1.0011 / -0.000 | 1.0000 / -0.000 | 1.0000 / -0.000 | 1.0000 / -0.000 | +| 8 | 0.9882 / 0 | 0.9998 / 0 | 1.0140 / -0.008 | 1.0002 / -0.000 | 0.9985 / -0.003 | 1.0000 / -0.000 | +| 6 | 0.9648 / 0 | 0.9983 / 0 | 1.0337 / -0.042 | 1.0013 / -0.002 | 0.9924 / -0.017 | 0.9996 / -0.001 | +| 4 | 0.8488 / 0 | 0.9701 / 0 | 1.0610 / -0.333 | 1.0065 / -0.057 | 0.9337 / -0.133 | 0.9863 / -0.025 | + +What this means in practice: + +- **`'poly'` does not damp anything.** $\mathrm{Re}\,\lambda = 0$ for every wavenumber, so for linear, constant-coefficient problems the semi-discrete scheme conserves the discrete $L^2$ energy; only the time integrator dissipates (slightly). The flip side is that unresolved modes, including grid-scale ($\theta = \pi$) noise, are never damped. Waves lag (phase speed below 1). +- **`'poly_dissipative'` damps most** and its waves run slightly fast (phase speed above 1). +- **`'poly_cell_centered'` damps less** and its waves lag. It has an exact relation to `'poly'`: the average of its two face states is exactly the `'poly'` face value of order N+2. In other words, `'poly_cell_centered'` of order N **is** `'poly'` of order N+2 plus an upwind correction proportional to the jump, which is why their dispersion curves coincide. + +## 7. Using them + +```python +solver = pyclaw.SharpClawSolver1D(riemann_solver) +solver.lim_type = 'poly' # or 4; 'poly_dissipative' (5), 'poly_cell_centered' (6) +solver.poly_order = 6 # even: 4, 6, 8 or 10 +``` + +The default `poly_order` is 4. Only `lim_type` selects the reconstruction. Each family keeps its own order when switching methods. + +Limitations: + +- **No limiting.** At discontinuities the reconstructions produce Gibbs oscillations; they are meant for smooth solutions. +- **Fortran kernels only**, applied component-wise (`char_decomp` is ignored). They work in 1D, 2D and 3D, since the 2D and 3D solvers apply the 1D reconstruction along each direction. +- **Energy conservation** of `'poly'` is shown here only for linear, constant-coefficient problems; with variable coefficients (e.g. material interfaces stored in `aux`) it is not guaranteed and should be checked for the problem at hand. + +## 8. Files and how to regenerate everything + +| file | purpose | +|---|---| +| `poly_weights.py` | exact weights, Fourier symbols, and the generator of `src/pyclaw/sharpclaw/poly.f90` | +| `sine_reconstruction.py` | the sine example; `--cells N` sets the grid, `--fortran` compares with the compiled code | +| `test_poly_reconstruction.py` | toy tests (run with `python -m pytest` in this directory) | +| `make_figures.py` | regenerates `figures/*.png` (NumPy and matplotlib only) | + +```bash +cd development/poly_reconstruction +python poly_weights.py --check # is poly.f90 up to date? +python poly_weights.py --write # regenerate it after changing the stencils +python sine_reconstruction.py # error and observed-order table +python -m pytest # toy tests +python make_figures.py # rebuild the figures (run by hand) +``` + +CI does not run anything in `development/`, but `src/pyclaw/tests/test_sharpclaw_poly.py` checks that `poly.f90` matches this generator and that the compiled weights equal the exact ones, so the two cannot drift apart unnoticed. diff --git a/development/poly_reconstruction/figures/averages_vs_points.png b/development/poly_reconstruction/figures/averages_vs_points.png new file mode 100644 index 00000000..8cbec865 Binary files /dev/null and b/development/poly_reconstruction/figures/averages_vs_points.png differ diff --git a/development/poly_reconstruction/figures/convergence.png b/development/poly_reconstruction/figures/convergence.png new file mode 100644 index 00000000..ab711fa1 Binary files /dev/null and b/development/poly_reconstruction/figures/convergence.png differ diff --git a/development/poly_reconstruction/figures/face_states.png b/development/poly_reconstruction/figures/face_states.png new file mode 100644 index 00000000..2306fd1f Binary files /dev/null and b/development/poly_reconstruction/figures/face_states.png differ diff --git a/development/poly_reconstruction/figures/fourier.png b/development/poly_reconstruction/figures/fourier.png new file mode 100644 index 00000000..2fd46051 Binary files /dev/null and b/development/poly_reconstruction/figures/fourier.png differ diff --git a/development/poly_reconstruction/figures/primitive_construction.png b/development/poly_reconstruction/figures/primitive_construction.png new file mode 100644 index 00000000..cf9c0a3c Binary files /dev/null and b/development/poly_reconstruction/figures/primitive_construction.png differ diff --git a/development/poly_reconstruction/figures/sine_reconstruction.png b/development/poly_reconstruction/figures/sine_reconstruction.png new file mode 100644 index 00000000..28bb5e2a Binary files /dev/null and b/development/poly_reconstruction/figures/sine_reconstruction.png differ diff --git a/development/poly_reconstruction/figures/stencils.png b/development/poly_reconstruction/figures/stencils.png new file mode 100644 index 00000000..daa26288 Binary files /dev/null and b/development/poly_reconstruction/figures/stencils.png differ diff --git a/development/poly_reconstruction/make_figures.py b/development/poly_reconstruction/make_figures.py new file mode 100644 index 00000000..1ed89af6 --- /dev/null +++ b/development/poly_reconstruction/make_figures.py @@ -0,0 +1,320 @@ +#!/usr/bin/env python +r""" +Regenerate the figures in figures/ used by README.md. + +Run manually after changing the reconstructions; the PNG files are +committed so that readers need nothing installed. Needs only NumPy and +matplotlib:: + + python make_figures.py + +Colour roles, the same in every figure: blue is the left-face value ql(i), +orange the right-face value qr(i); the orders N = 4, 6, 8, 10 use a +light-to-dark blue ramp; the three reconstruction types get one panel +each rather than a colour. +""" +import os +from fractions import Fraction + +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt # noqa: E402 +import numpy as np # noqa: E402 +from matplotlib.patches import Rectangle # noqa: E402 + +import poly_weights as pw # noqa: E402 +from sine_reconstruction import interface_error, sine_averages # noqa: E402 + +HERE = os.path.dirname(os.path.abspath(__file__)) +OUT = os.path.join(HERE, 'figures') + +# Palette: categorical slots 1-2 and an ordinal blue ramp, validated +# (CVD and contrast) on the light surface below. +DPI = 150 + +SURFACE = '#fcfcfb' +INK = '#0b0b0b' +INK_2 = '#52514e' +MUTED = '#a3a29c' +GRID = '#e6e5e0' +QL = '#2a78d6' # left-face value ql(i) +QR = '#eb6834' # right-face value qr(i) +QL_FILL = '#cde2fb' +QR_FILL = '#fbd9c9' +ORDER_RAMP = {4: '#86b6ef', 6: '#3987e5', 8: '#1c5cab', 10: '#0d366b'} + +TITLES = {4: "lim_type 4, 'poly'", + 5: "lim_type 5, 'poly_dissipative'", + 6: "lim_type 6, 'poly_cell_centered'"} + +plt.rcParams.update({ + 'figure.facecolor': SURFACE, 'axes.facecolor': SURFACE, 'savefig.facecolor': SURFACE, + 'axes.edgecolor': MUTED, 'axes.labelcolor': INK_2, 'axes.titlecolor': INK, + 'xtick.color': INK_2, 'ytick.color': INK_2, 'text.color': INK, + 'axes.grid': True, 'grid.color': GRID, 'grid.linewidth': 0.8, + 'axes.spines.top': False, 'axes.spines.right': False, + 'font.size': 10, 'axes.titlesize': 10.5, 'axes.titleweight': 'bold', + 'lines.linewidth': 2, 'legend.frameon': False, +}) + + +def save(fig, name): + os.makedirs(OUT, exist_ok=True) + path = os.path.join(OUT, name) + # no Software tag, so the files do not change with the matplotlib version + fig.savefig(path, dpi=DPI, metadata={'Software': None}, bbox_inches='tight') + plt.close(fig) + print('wrote', os.path.relpath(path, HERE)) + + +# ---------------------------------------------------------------------------- +def fig_averages_vs_points(): + """Cell averages are not centre values, and treating them as such + caps the old weights at second order.""" + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 3.8), gridspec_kw={'width_ratios': [1, 1]}) + + n = 4 + h = 1. / n + x = np.linspace(0, 1, 400) + avg = sine_averages(n) + xc = (np.arange(n) + 0.5) * h + ax1.plot(x, np.sin(2*np.pi*x), color=MUTED, lw=1.5, label='q(x) = sin 2πx') + for j in range(n): + ax1.add_patch(Rectangle((j*h, 0), h, avg[j], facecolor=QL_FILL, edgecolor=SURFACE, lw=2)) + ax1.plot(xc, avg, 's', color=QL, ms=6, label='cell average (what PyClaw stores)') + ax1.plot(xc, np.sin(2*np.pi*xc), 'o', color=INK, ms=4, label='value at the cell centre') + ax1.set_xlabel('x') + ax1.set_xlim(0, 1) + ax1.set_ylim(-1.15, 1.15) + for j in range(n): + ax1.annotate('', (xc[j], avg[j]), (xc[j], np.sin(2*np.pi*xc[j])), + arrowprops=dict(arrowstyle='<->', color=INK_2, lw=1, shrinkA=4, shrinkB=4)) + ax1.legend(loc='lower left', fontsize=8.5) + ax1.set_title('Average ≠ centre value: gap ≈ h²q″/24', loc='left') + + cells_n = [8, 16, 32, 64, 128] + for order in pw.ORDERS: + cells = list(pw.stencil(6, order)) + w_pt = [float(v) for v in pw.point_weights(cells, Fraction(-1, 2))] + errors = [] + for m in cells_n: + hh = 1. / m + xx = (np.arange(m) + 0.5) * hh + ql = sum(wj * np.roll(sine_averages(m), -j) for j, wj in zip(cells, w_pt)) + errors.append(np.abs(ql - np.sin(2*np.pi*(xx - hh/2))).max()) + ax2.loglog(cells_n, errors, color=QR, lw=1.5, alpha=0.9) + fv = np.array([interface_error(m, 6, order) for m in cells_n]) + keep = fv > 1e-13 + ax2.loglog(np.array(cells_n)[keep], fv[keep], color=ORDER_RAMP[order], marker='o', ms=4) + ax2.text(140, 2e-5, 'point-value weights\n(old), all N:\nslope 2', color=QR, fontsize=8.5, va='center') + for order in pw.ORDERS: + ax2.plot([], [], color=ORDER_RAMP[order], marker='o', ms=4, + label='cell-average, N=%d: slope %d' % (order, order + 1)) + ax2.legend(loc='center left', bbox_to_anchor=(1.0, 0.45), fontsize=8) + ax2.set_xlim(7, 600) + ax2.set_ylim(1e-14, 1) + ax2.set_xlabel('number of cells') + ax2.set_ylabel('max interface error') + ax2.set_title('Old vs new weights (2k+1-cell stencils)', loc='left') + save(fig, 'averages_vs_points.png') + + +# ---------------------------------------------------------------------------- +def fig_primitive_construction(): + """Interpolate the primitive at the faces, differentiate at the face.""" + cells = list(pw.stencil(6, 4)) # i-2 .. i+2 + q = lambda s: 1.0 + 0.5*np.sin(0.9*s + 0.4) + 0.08*s # noqa: E731 + big_q = lambda s: 1.0*s - 0.5/0.9*np.cos(0.9*s + 0.4) + 0.04*s**2 # noqa: E731 + faces = np.arange(cells[0] - 0.5, cells[-1] + 1.0) + avg = [big_q(j + 0.5) - big_q(j - 0.5) for j in cells] + prim = np.concatenate([[0], np.cumsum(avg)]) + target = -0.5 + xs = np.linspace(faces[0] - 0.3, faces[-1] + 0.3, 400) + + fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(7.5, 6.2), sharex=True, + gridspec_kw={'height_ratios': [1, 1.15]}) + for j, a in zip(cells, avg): + ax1.add_patch(Rectangle((j - 0.5, 0), 1, a, facecolor=QL_FILL, edgecolor=SURFACE, lw=2)) + ax1.text(j, a + 0.05, r'$\bar{q}_{i%s}$' % ('' if j == 0 else '%+d' % j), ha='center', fontsize=10, color=INK_2) + ax1.plot(xs, q(xs), color=MUTED, lw=1.5) + ax1.plot([target], [q(target)], 'o', color=QL, ms=8, zorder=5) + ax1.annotate('wanted: ql(i) = q(x$_{i-1/2}$)', (target, q(target)), (target - 2.3, 1.85), + color=QL, fontsize=9, arrowprops=dict(arrowstyle='->', color=QL, lw=1)) + ax1.set_ylim(0, 2.1) + ax1.set_ylabel('q') + ax1.set_title('1. Known: the cell averages over the stencil', loc='left') + + coeffs = np.polyfit(faces, prim, len(faces) - 1) + ax2.plot(xs, np.polyval(coeffs, xs), color=INK_2, lw=1.5, label='interpolating polynomial P(x)') + ax2.plot(faces, prim, 'o', color=INK, ms=6, label=r'Q at the faces = running sum of $h\,\bar{q}$') + slope = np.polyval(np.polyder(coeffs), target) + t = np.linspace(target - 1.1, target + 1.1, 2) + ax2.plot(t, np.polyval(coeffs, target) + slope*(t - target), color=QL, lw=2.5, + label="tangent at x$_{i-1/2}$: slope P′ = ql(i)") + ax2.set_ylabel('Q(x) = ∫ q dx') + ax2.set_title('2. Interpolate the primitive at the faces; its slope at the face is ql(i)', loc='left') + ax2.legend(loc='upper left', fontsize=8.5) + ax2.set_xticks(faces) + ax2.set_xticklabels(['x$_{i%s}$' % ('-1/2' if f == -0.5 else '%+g' % f).replace('+0.5', '+1/2') + .replace('-1.5', '-3/2').replace('+1.5', '+3/2').replace('-2.5', '-5/2') + .replace('+2.5', '+5/2') for f in faces]) + for ax in (ax1, ax2): + ax.axvline(target, color=QL, lw=0.8, ls=':') + fig.text(0.01, -0.02, 'Shown for the 5-cell stencil of lim_type 6, N = 4. The weights do not depend ' + 'on q: they follow from this construction once, exactly, in poly_weights.py.', + fontsize=8.5, color=INK_2) + save(fig, 'primitive_construction.png') + + +# ---------------------------------------------------------------------------- +def draw_strip(ax, y, used, colour, fill, label, face=None, span=5, centre_cell=True): + """Row of cells i-span..i+span; cells in `used` are filled.""" + for j in range(-span, span + 1): + on = j in used + ax.add_patch(Rectangle((j - 0.5, y - 0.35), 1, 0.7, facecolor=fill if on else SURFACE, + edgecolor=colour if on else GRID, lw=1.6 if on else 1)) + if on: + ax.text(j, y, 'i' if j == 0 else 'i%+d' % j, ha='center', va='center', fontsize=7.5, color=INK) + if centre_cell: + ax.add_patch(Rectangle((-0.5, y - 0.35), 1, 0.7, fill=False, edgecolor=INK, lw=1.6)) + if face is not None: + ax.plot([face, face], [y - 0.5, y + 0.5], color=colour, lw=3, solid_capstyle='butt') + ax.text(-span - 0.8, y, label, ha='right', va='center', fontsize=9, color=colour) + + +def fig_stencils(): + """Which cells each type uses for ql(i) and qr(i).""" + fig, axes = plt.subplots(3, 2, figsize=(11, 5.6)) + for row, lim_type in enumerate(pw.LIM_TYPES): + for col, order in enumerate((4, 6)): + ax = axes[row, col] + cells = list(pw.stencil(lim_type, order)) + draw_strip(ax, 1, cells, QL, QL_FILL, 'ql(i)', face=-0.5) + draw_strip(ax, 0, [-j for j in cells], QR, QR_FILL, 'qr(i)', face=0.5) + ax.set_xlim(-8.6, 5.7) + ax.set_ylim(-0.6, 1.6) + ax.axis('off') + ax.set_title('%s, N = %d: %d cells, order %d' + % (TITLES[lim_type], order, len(cells), pw.formal_order(lim_type, order)), + loc='left', fontsize=9.5) + fig.text(0.01, -0.02, 'Bold outline: cell i. Thick bar: the face whose value is computed ' + '(x$_{i-1/2}$ for ql, x$_{i+1/2}$ for qr). qr(i) uses the mirror image of the ql(i) cells.', + fontsize=8.5, color=INK_2) + fig.tight_layout() + save(fig, 'stencils.png') + + +def fig_face_states(): + """The two states the Riemann solver sees at one face.""" + fig, axes = plt.subplots(3, 1, figsize=(7.5, 5.6)) + notes = {4: 'same cells, same weights: one value, no jump, no dissipation', + 5: 'each side leans upwind: the states differ, the jump dissipates', + 6: 'each side centred on its own cell: the states differ, the jump dissipates'} + for ax, lim_type in zip(axes, pw.LIM_TYPES): + cells = list(pw.stencil(lim_type, 4)) + # face x_{i+1/2}: left state qr(i), right state ql(i+1) + draw_strip(ax, 1, [-j for j in cells], QR, QR_FILL, 'left state qr(i)', face=0.5) + draw_strip(ax, 0, [1 + j for j in cells], QL, QL_FILL, 'right state ql(i+1)', face=0.5, + centre_cell=False) + ax.set_xlim(-11.5, 5.7) + ax.set_ylim(-0.6, 1.6) + ax.axis('off') + ax.set_title('%s (N = 4): %s' % (TITLES[lim_type], notes[lim_type]), loc='left', fontsize=9) + fig.text(0.01, -0.02, 'Both rows reconstruct the value at the same face x$_{i+1/2}$ (thick bar), ' + 'from the left and from the right.', fontsize=8.5, color=INK_2) + fig.tight_layout() + save(fig, 'face_states.png') + + +# ---------------------------------------------------------------------------- +def fig_sine_reconstruction(): + """ql and qr of sin(4 pi x) from 8 cell averages (4 per wavelength), N = 4.""" + n = 8 + h = 1. / n + x = np.linspace(0, 1, 400) + avg = sine_averages(n, kwave=2) + fig, axes = plt.subplots(1, 3, figsize=(11, 3.6), sharey=True) + for ax, lim_type in zip(axes, pw.LIM_TYPES): + ql, qr = pw.reconstruct(avg, lim_type, 4) + for j in range(n): + ax.add_patch(Rectangle((j*h, min(0, avg[j])), h, abs(avg[j]), facecolor='#f0efec', + edgecolor=SURFACE, lw=2)) + ax.plot(x, np.sin(4*np.pi*x), color=MUTED, lw=1.5) + faces = np.arange(n) * h + ax.plot(faces, ql, '>', color=QL, ms=7, label='ql(i) at x$_{i-1/2}$') + ax.plot(faces + h, qr, '<', color=QR, ms=7, label='qr(i) at x$_{i+1/2}$') + jump = np.abs(np.roll(ql, -1) - qr).max() + ax.set_title('%s, N = 4' % TITLES[lim_type], loc='left', fontsize=9.5) + ax.text(0.03, -1.32, 'max |ql(i+1) − qr(i)| = %.1e' % jump, fontsize=8.5, color=INK_2) + ax.set_xlabel('x') + ax.set_xlim(0, 1) + ax.set_ylim(-1.45, 1.3) + axes[0].legend(loc='upper right', fontsize=8.5) + fig.text(0.01, -0.04, 'Grey: exact sin 4πx and its 8 cell averages, a coarse 4 cells per wavelength. ' + 'Each face gets a value from each neighbouring cell; with \'poly\' they coincide.', + fontsize=8.5, color=INK_2) + fig.tight_layout() + save(fig, 'sine_reconstruction.png') + + +def fig_convergence(): + """Interface error against the number of cells.""" + cells_n = np.array([4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128]) + fig, axes = plt.subplots(1, 3, figsize=(11, 3.8), sharey=True) + for ax, lim_type in zip(axes, pw.LIM_TYPES): + for order in pw.ORDERS: + errors = np.array([interface_error(m, lim_type, order) for m in cells_n]) + keep = errors > 1e-13 + ax.loglog(cells_n[keep], errors[keep], color=ORDER_RAMP[order], marker='o', ms=3.5) + last = np.flatnonzero(keep)[-1] + ax.text(cells_n[last] * 1.12, errors[last], 'N=%d' % order, fontsize=8, + color=INK_2, va='center') + ax.set_title(TITLES[lim_type], loc='left', fontsize=9.5) + ax.set_xlabel('number of cells') + ax.set_xlim(3.5, 250) + ax.set_ylim(1e-14, 2) + axes[0].set_ylabel('max interface error, sin 2πx') + fig.text(0.01, -0.04, 'Slopes: N for \'poly\' and \'poly_dissipative\', N+1 for ' + '\'poly_cell_centered\'. Curves stop where round-off takes over.', + fontsize=8.5, color=INK_2) + fig.tight_layout() + save(fig, 'convergence.png') + + +def fig_fourier(): + """Dissipation and dispersion for u_t + u_x = 0 (semi-discrete).""" + theta = np.linspace(1e-4, np.pi, 600) + fig, axes = plt.subplots(2, 3, figsize=(11, 6.2), sharex=True, sharey='row') + for col, lim_type in enumerate(pw.LIM_TYPES): + top, bottom = axes[0, col], axes[1, col] + for order in pw.ORDERS: + lam = pw.fourier_symbol(lim_type, order, theta) + top.plot(theta, lam.real, color=ORDER_RAMP[order], label='N=%d' % order) + bottom.plot(theta, -lam.imag, color=ORDER_RAMP[order], label='N=%d' % order) + bottom.plot(theta, theta, color=MUTED, lw=1.2, ls='--', label='exact') + top.set_title(TITLES[lim_type], loc='left', fontsize=9.5) + if lim_type == 4: + top.text(0.1, -0.35, 'Re λ = 0 for every N:\nno numerical dissipation', fontsize=9, color=INK_2) + bottom.set_xlabel('θ = kh (points per wavelength = 2π/θ)') + bottom.set_xticks([0, np.pi/4, np.pi/2, 3*np.pi/4, np.pi]) + bottom.set_xticklabels(['0', 'π/4\n(8)', 'π/2\n(4)', '3π/4\n(2.7)', 'π\n(2)']) + axes[0, 0].set_ylabel('Re λ (damping rate)') + axes[1, 0].set_ylabel('−Im λ (modified wavenumber)') + axes[0, 1].legend(loc='lower left', fontsize=8.5) + axes[1, 0].legend(loc='upper left', fontsize=8.5) + axes[0, 0].set_ylim(-2.8, 0.15) + fig.text(0.01, -0.03, 'Mode exp(ijθ) with h = 1 and upwind flux evolves as exp(λt); exact: ' + 'λ = −iθ. Re λ = 0 means no damping; −Im λ below the diagonal means waves lag.', + fontsize=8.5, color=INK_2) + fig.tight_layout() + save(fig, 'fourier.png') + + +if __name__ == '__main__': + fig_averages_vs_points() + fig_primitive_construction() + fig_stencils() + fig_face_states() + fig_sine_reconstruction() + fig_convergence() + fig_fourier() diff --git a/development/poly_reconstruction/poly_weights.py b/development/poly_reconstruction/poly_weights.py new file mode 100644 index 00000000..c07e7319 --- /dev/null +++ b/development/poly_reconstruction/poly_weights.py @@ -0,0 +1,331 @@ +#!/usr/bin/env python +r""" +Exact weights for the SharpClaw polynomial reconstructions, and the +generator of ``src/pyclaw/sharpclaw/poly.f90``. + +PyClaw stores cell averages. A reconstruction computes, for each cell i, +the value ql(i) at its left face x_{i-1/2} and qr(i) at its right face +x_{i+1/2}, as a linear combination of the averages of nearby cells: + + ql(i) = sum_j w_j q(i+j), qr(i) = sum_j w_j q(i-j). + +The weights are found by interpolating the primitive function +Q(x) = int_{x_lo}^{x} q dx, whose values at the faces are exact sums of +cell averages, and differentiating the interpolant at the face. They are +exact rationals; ``poly.f90`` stores each set as integers over a common +denominator. + +lim_type selects the cells (offsets j) used for ql(i); with k = N/2 and +N = reconstruction_order in {4, 6, 8, 10}: + + 4 'poly' j = -k .. k-1 (N cells centred on the face) + 5 'poly_dissipative' j = -k+1 .. k (central stencil shifted + one cell upwind) + 6 'poly_cell_centered' j = -k .. k (N+1 cells centred on the cell) + +Usage:: + + python poly_weights.py --table # print the weights as Markdown tables + python poly_weights.py --check # exit 1 if poly.f90 is out of date + python poly_weights.py --write # regenerate poly.f90 +""" +import argparse +import os +import sys +from fractions import Fraction +from math import lcm + +import numpy as np + +LIM_TYPES = {4: 'poly', 5: 'poly_dissipative', 6: 'poly_cell_centered'} +ORDERS = (4, 6, 8, 10) +DESCRIPTIONS = {4: 'central, non-dissipative', + 5: 'upwind-biased, dissipative', + 6: 'cell-centred, dissipative'} + +HERE = os.path.dirname(os.path.abspath(__file__)) +POLY_F90 = os.path.join(HERE, '..', '..', 'src', 'pyclaw', 'sharpclaw', 'poly.f90') + + +# ---------------------------------------------------------------------------- +# Weights +# ---------------------------------------------------------------------------- +def formal_order(lim_type, order): + """Formal order of accuracy of the interface values.""" + return order + 1 if lim_type == 6 else order + + +def stencil(lim_type, order): + """Offsets j of the cells used for ql(i) (value at x_{i-1/2}).""" + if lim_type not in LIM_TYPES or order not in ORDERS: + raise ValueError('lim_type must be in %s and order in %s' + % (tuple(LIM_TYPES), ORDERS)) + k = order // 2 + return {4: range(-k, k), 5: range(-k + 1, k + 1), 6: range(-k, k + 1)}[lim_type] + + +def face_weights(cells, x): + """Weights giving the value at x from the averages of the unit cells + centred at the integers in ``cells`` (cell j spans [j-1/2, j+1/2]). + + Interpolates the primitive function at the cell faces with a Lagrange + polynomial and returns the weights of its derivative at x. + """ + cells = list(cells) + x = Fraction(x) + faces = [Fraction(2*cells[0] - 1, 2) + m for m in range(len(cells) + 1)] + # d/dx of the Lagrange basis polynomial of each face, at x + dbasis = [] + for k, fk in enumerate(faces): + total = Fraction(0) + for a, fa in enumerate(faces): + if a == k: + continue + term = Fraction(1) / (fk - fa) + for b, fb in enumerate(faces): + if b not in (a, k): + term *= (x - fb) / (fk - fb) + total += term + dbasis.append(total) + # Q(face_k) = sum of the averages of the cells left of face k + return [sum(dbasis[k] for k in range(len(faces)) if k > j) + for j in range(len(cells))] + + +def point_weights(cells, x): + """Lagrange weights that treat q as point values at the cell centres + (the incorrect, second-order form used before; for comparison only).""" + cells = list(cells) + x = Fraction(x) + weights = [] + for j in cells: + w = Fraction(1) + for m in cells: + if m != j: + w *= (x - m) / Fraction(j - m) + weights.append(w) + return weights + + +def ql_weights(lim_type, order): + """(offsets, exact weights) for ql(i).""" + cells = list(stencil(lim_type, order)) + return cells, face_weights(cells, Fraction(-1, 2)) + + +def integer_weights(lim_type, order): + """(lo, hi, numerators, denominator) for ql(i).""" + cells, w = ql_weights(lim_type, order) + denominator = lcm(*[v.denominator for v in w]) + numerators = [int(v * denominator) for v in w] + return cells[0], cells[-1], numerators, denominator + + +# ---------------------------------------------------------------------------- +# Numerical helpers (used by the example, the tests and the figures) +# ---------------------------------------------------------------------------- +def reconstruct(q, lim_type, order): + """ql, qr for a periodic array of cell averages q (any length).""" + cells, w = ql_weights(lim_type, order) + w = [float(v) for v in w] + ql = sum(wj * np.roll(q, -j) for j, wj in zip(cells, w)) + qr = sum(wj * np.roll(q, j) for j, wj in zip(cells, w)) + return ql, qr + + +def fourier_symbol(lim_type, order, theta): + r"""Semi-discrete symbol lambda(theta) for u_t + u_x = 0 with h = 1 and + the upwind flux, so that dq/dt = lambda q for q = exp(i j theta). + The exact value is -i theta; Re(lambda) < 0 is dissipation.""" + cells, w = ql_weights(lim_type, order) + theta = np.asarray(theta, dtype=float) + # left state at face i+1/2 is qr(i) = sum_j w_j q(i-j) + qr = sum(float(wj) * np.exp(-1j * j * theta) for j, wj in zip(cells, w)) + return -(1 - np.exp(-1j * theta)) * qr + + +# ---------------------------------------------------------------------------- +# Fortran generator +# ---------------------------------------------------------------------------- +HEADER = """\ +! =================================================================== +module poly +! =================================================================== +! Linear (unlimited) polynomial reconstruction of interface values +! from cell averages, for the SharpClaw finite-volume solvers. +! +! q(:,i) holds the average of q over cell i. For each cell, ql(:,i) is +! the reconstructed value at the left face x_(i-1/2) and qr(:,i) the +! value at the right face x_(i+1/2). The weights are exact: they are +! obtained by differentiating the polynomial that interpolates the +! primitive function of q at the cell faces, so they are correct for +! cell averages (not point values). They are applied component-wise. +! +! lim_type selects the stencil; order = solver.poly_order, +! an even number N in {4, 6, 8, 10}. With k = N/2: +! +! 4 ('poly') ql(i) uses cells i-k..i+k-1 (N cells, +! centred on the face). Order N. ql(i+1) +! equals qr(i), so the scheme adds no +! numerical dissipation. +! 5 ('poly_dissipative') ql(i) uses cells i-k+1..i+k (the central +! stencil shifted one cell upwind for each +! face state). Order N, dissipative. +! 6 ('poly_cell_centered') ql(i) uses cells i-k..i+k (N+1 cells, +! centred on the cell, polynomial degree N). +! Order N+1, dissipative; equal to WENO of +! order N+1 with its ideal linear weights. +! +! qr uses the mirror image of the ql stencil. Every variant reaches k +! cells from the cell being reconstructed, so num_ghost = k + 1. +! +! Coefficients generated with exact rational arithmetic. +! =================================================================== + + implicit none + +contains + + ! =================================================================== + subroutine poly_weights(lim_type, order, lo, hi, w, ok) + ! =================================================================== + ! Weights for the left-face value ql(i) = sum_j w(j) * q(i+j), + ! j = lo..hi. + + integer, intent(in) :: lim_type, order + integer, intent(out) :: lo, hi + double precision, intent(out) :: w(-5:5) + logical, intent(out) :: ok + + ok = .true. + lo = 0; hi = 0 + w = 0.d0 + + select case (lim_type) +""" + +FOOTER = """\ + case default + ok = .false. + end select + + end subroutine poly_weights + + ! =================================================================== + subroutine poly_recon(q,ql,qr,num_eqn,maxnx,num_ghost,lim_type,order) + ! =================================================================== + + integer, intent(in) :: num_eqn, maxnx, num_ghost, lim_type, order + double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost) + double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost) + + double precision :: w(-5:5) + integer :: lo, hi, reach, i, j, mx2 + logical :: ok + + call poly_weights(lim_type, order, lo, hi, w, ok) + if (.not. ok) then + print *, 'ERROR: polynomial reconstruction needs lim_type 4, 5 or 6' + print *, ' and an even reconstruction_order between 4 and 10.' + error stop + endif + + mx2 = size(q,2) + reach = max(-lo, hi) + if (reach >= num_ghost) then + print *, 'ERROR: polynomial reconstruction needs num_ghost >= ', reach+1 + error stop + endif + + ql = 0.d0 + qr = 0.d0 + do i = 1+reach, mx2-reach + do j = lo, hi + ql(:,i) = ql(:,i) + w(j)*q(:,i+j) + ! qr mirrors the ql stencil about cell i + qr(:,i) = qr(:,i) + w(j)*q(:,i-j) + end do + end do + + end subroutine poly_recon + +end module poly +""" + + +def _fortran_array(numerators, denominator): + values = ['%d.d0' % v for v in numerators] + if len(values) <= 8: + return '[' + ', '.join(values) + '] / %d.d0' % denominator + half = len(values) // 2 + return ('[' + ', '.join(values[:half]) + ', &\n' + ' '*22 + + ', '.join(values[half:]) + '] / %d.d0' % denominator) + + +def fortran_source(): + """The full text of poly.f90.""" + cases = [] + for lim_type in LIM_TYPES: + orders = [] + for order in ORDERS: + lo, hi, numerators, denominator = integer_weights(lim_type, order) + orders.append(' case (%d)\n' + ' lo = %d; hi = %d\n' + ' w(%d:%d) = %s' + % (order, lo, hi, lo, hi, _fortran_array(numerators, denominator))) + cases.append(' case (%d) ! %s\n' + ' select case (order)\n' + '%s\n' + ' case default\n' + ' ok = .false.\n' + ' end select' + % (lim_type, DESCRIPTIONS[lim_type], '\n'.join(orders))) + return HEADER + '\n'.join(cases) + '\n' + FOOTER + + +def markdown_tables(): + """The weights of every scheme as Markdown tables.""" + out = [] + for lim_type, name in LIM_TYPES.items(): + out.append("**lim_type = %d, `'%s'`** (weights for ql(i), the value at " + "xi-1/2; qr(i) uses the mirror image)\n" % (lim_type, name)) + out.append('| N | order | cells | weights |') + out.append('|---|---|---|---|') + for order in ORDERS: + lo, hi, numerators, denominator = integer_weights(lim_type, order) + cells = 'i%+d … i%+d' % (lo, hi) if lo else 'i … i%+d' % hi + out.append('| %d | %d | %s | (%s) / %d |' + % (order, formal_order(lim_type, order), cells, + ', '.join(str(v) for v in numerators), denominator)) + out.append('') + return '\n'.join(out) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__.split('\n\n')[0]) + group = parser.add_mutually_exclusive_group(required=True) + group.add_argument('--write', action='store_true', help='regenerate poly.f90') + group.add_argument('--check', action='store_true', + help='exit with status 1 if poly.f90 is out of date') + group.add_argument('--table', action='store_true', help='print Markdown tables') + args = parser.parse_args(argv) + + if args.table: + print(markdown_tables()) + return 0 + if args.write: + with open(POLY_F90, 'w') as f: + f.write(fortran_source()) + print('Wrote %s' % os.path.normpath(POLY_F90)) + return 0 + with open(POLY_F90) as f: + current = f.read() + if current == fortran_source(): + print('poly.f90 is up to date.') + return 0 + print('poly.f90 differs from the generator; run with --write.') + return 1 + + +if __name__ == '__main__': + sys.exit(main()) diff --git a/development/poly_reconstruction/sine_reconstruction.py b/development/poly_reconstruction/sine_reconstruction.py new file mode 100644 index 00000000..45d8a8b2 --- /dev/null +++ b/development/poly_reconstruction/sine_reconstruction.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python +r""" +Toy example: reconstruct sin(2 pi x) on a periodic grid of N cells with +every polynomial reconstruction and order, and report the interface error +and the observed order of accuracy. + +The input is the exact cell averages of sin(2 pi x); the output is the +value at every cell face, compared with the exact sin(2 pi x_{i -/+ 1/2}). + +Usage:: + + python sine_reconstruction.py # N = 16, NumPy weights + python sine_reconstruction.py --cells 32 + python sine_reconstruction.py --fortran # also check the compiled + # Fortran (clawpack built) +""" +import argparse + +import numpy as np + +from poly_weights import LIM_TYPES, ORDERS, formal_order, reconstruct + + +def sine_averages(num_cells, kwave=1): + """Exact averages of sin(2 pi kwave x) over the cells of [0, 1].""" + h = 1. / num_cells + xc = (np.arange(num_cells) + 0.5) * h + w = 2 * np.pi * kwave + return (np.cos(w * (xc - h/2)) - np.cos(w * (xc + h/2))) / (w * h) + + +def interface_error(num_cells, lim_type, order): + """Max error of ql and qr against the exact face values.""" + h = 1. / num_cells + xc = (np.arange(num_cells) + 0.5) * h + ql, qr = reconstruct(sine_averages(num_cells), lim_type, order) + return max(np.abs(ql - np.sin(2*np.pi*(xc - h/2))).max(), + np.abs(qr - np.sin(2*np.pi*(xc + h/2))).max()) + + +def fortran_difference(num_cells, lim_type, order): + """Max difference between the NumPy and the compiled Fortran result.""" + from clawpack.pyclaw.sharpclaw import sharpclaw1 + ng = order // 2 + 1 + q = sine_averages(num_cells) + padded = np.concatenate([q[-ng:], q, q[:ng]])[None, :] + ql_f, qr_f = sharpclaw1.poly.poly_recon(padded, num_cells, ng, lim_type, order) + ql, qr = reconstruct(q, lim_type, order) + real = slice(ng, num_cells + ng) + return max(np.abs(ql_f[0, real] - ql).max(), np.abs(qr_f[0, real] - qr).max()) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__.split('\n\n')[0]) + parser.add_argument('--cells', type=int, default=16, help='number of cells N (default 16)') + parser.add_argument('--fortran', action='store_true', + help='compare with the compiled Fortran reconstruction') + args = parser.parse_args(argv) + n = args.cells + + print('Interface error for sin(2 pi x) from %d cell averages; the observed order' + ' compares %d and %d cells.\n' % (n, n, 2*n)) + header = '%-20s %3s %7s %12s %9s' % ('lim_type', 'N', 'formal', 'max error', 'observed') + if args.fortran: + header += ' %14s' % '|NumPy-Fortran|' + print(header) + print('-' * len(header)) + for lim_type, name in LIM_TYPES.items(): + for order in ORDERS: + e1 = interface_error(n, lim_type, order) + e2 = interface_error(2*n, lim_type, order) + line = '%-20s %3d %7d %12.3e %9.2f' % ('%d %s' % (lim_type, name), order, + formal_order(lim_type, order), e1, + np.log2(e1 / e2)) + if args.fortran: + line += ' %14.1e' % fortran_difference(n, lim_type, order) + print(line) + print('\nObserved orders fall below the formal order once the error nears ' + 'round-off (about 1e-14); use fewer cells to see the high orders.') + + +if __name__ == '__main__': + main() diff --git a/development/poly_reconstruction/test_poly_reconstruction.py b/development/poly_reconstruction/test_poly_reconstruction.py new file mode 100644 index 00000000..0657cbaa --- /dev/null +++ b/development/poly_reconstruction/test_poly_reconstruction.py @@ -0,0 +1,71 @@ +""" +Toy tests for the polynomial reconstructions, using the exact weights from +poly_weights.py. Run from this directory with ``python -m pytest``. + +CI does not collect development/; the Fortran and the generator are +covered there by src/pyclaw/tests/test_sharpclaw_poly.py. +""" +from fractions import Fraction + +import numpy as np +import pytest + +import poly_weights as pw +from sine_reconstruction import interface_error + +SCHEMES = [(lim_type, order) for lim_type in pw.LIM_TYPES for order in pw.ORDERS] + + +@pytest.mark.parametrize("lim_type, order", SCHEMES) +def test_sine_order_of_accuracy(lim_type, order): + grids = (16, 32) if order < 8 else (8, 16) + coarse, fine = (interface_error(n, lim_type, order) for n in grids) + assert np.log2(coarse / fine) > pw.formal_order(lim_type, order) - 0.4 + + +@pytest.mark.parametrize("lim_type, order", SCHEMES) +def test_weights_are_consistent(lim_type, order): + """Weights sum to one and reproduce averages of polynomials of the + degree the stencil supports exactly.""" + cells, w = pw.ql_weights(lim_type, order) + assert sum(w) == 1 + degree = len(cells) - 1 + for p in range(degree + 1): + # average of x^p over cell j is ((j+1/2)^(p+1) - (j-1/2)^(p+1)) / (p+1) + averages = [(Fraction(2*j + 1, 2)**(p + 1) - Fraction(2*j - 1, 2)**(p + 1)) / (p + 1) + for j in cells] + assert sum(wj * a for wj, a in zip(w, averages)) == Fraction(-1, 2)**p + + +@pytest.mark.parametrize("order", pw.ORDERS) +def test_poly_is_non_dissipative(order): + theta = np.linspace(0, np.pi, 201) + assert np.abs(pw.fourier_symbol(4, order, theta).real).max() < 1e-13 + + +@pytest.mark.parametrize("lim_type", (5, 6)) +@pytest.mark.parametrize("order", pw.ORDERS) +def test_dissipative_variants_are_stable_and_dissipative(lim_type, order): + theta = np.linspace(1e-3, np.pi, 201) + real = pw.fourier_symbol(lim_type, order, theta).real + assert real.max() < 1e-13 + assert real[-1] < -0.5 + + +def test_point_value_weights_are_second_order(): + """The old weights treated averages as point values: second order.""" + cells = list(pw.stencil(6, 4)) + w = [float(v) for v in pw.point_weights(cells, Fraction(-1, 2))] + errors = [] + for n in (16, 32): + h = 1. / n + xc = (np.arange(n) + 0.5) * h + avg = (np.cos(2*np.pi*(xc - h/2)) - np.cos(2*np.pi*(xc + h/2))) / (2*np.pi*h) + ql = sum(wj * np.roll(avg, -j) for j, wj in zip(cells, w)) + errors.append(np.abs(ql - np.sin(2*np.pi*(xc - h/2))).max()) + assert abs(np.log2(errors[0] / errors[1]) - 2) < 0.1 + + +def test_generated_fortran_is_current(): + with open(pw.POLY_F90) as f: + assert f.read() == pw.fortran_source() diff --git a/examples/advection_1d/advection_1d.py b/examples/advection_1d/advection_1d.py index 84765d21..71c5ca2d 100755 --- a/examples/advection_1d/advection_1d.py +++ b/examples/advection_1d/advection_1d.py @@ -19,7 +19,7 @@ import numpy as np from clawpack import riemann -def setup(nx=400, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, +def setup(nx=400, kernel_language='Python', use_petsc=False, solver_type='classic', weno_order=5, time_integrator='SSP104', outdir='./_output'): if use_petsc: diff --git a/examples/euler_1d/Makefile b/examples/euler_1d/Makefile index a4ad15d5..40e42fae 100644 --- a/examples/euler_1d/Makefile +++ b/examples/euler_1d/Makefile @@ -2,7 +2,7 @@ F2PY = f2py SHARPCLAW = ../../src/pyclaw/sharpclaw -ONE_D_SHARPCLAW_SOURCES = evec.f90 $(SHARPCLAW)/flux1.f90 $(SHARPCLAW)/ClawParams.f90 $(SHARPCLAW)/workspace.f90 $(SHARPCLAW)/weno.f90 $(SHARPCLAW)/reconstruct.f90 +ONE_D_SHARPCLAW_SOURCES = evec.f90 $(SHARPCLAW)/flux1.f90 $(SHARPCLAW)/ClawParams.f90 $(SHARPCLAW)/workspace.f90 $(SHARPCLAW)/weno.f90 $(SHARPCLAW)/poly.f90 $(SHARPCLAW)/reconstruct.f90 all: make sharpclaw1.so diff --git a/examples/euler_1d/setup.py b/examples/euler_1d/setup.py index 1dbad212..944aef9f 100644 --- a/examples/euler_1d/setup.py +++ b/examples/euler_1d/setup.py @@ -21,6 +21,7 @@ def configuration(parent_package='',top_path=None): sharpclaw_srcs = [pjoin(sharpclaw_dir, src) for src in ['ClawParams.f90', 'workspace.f90', 'weno.f90', + 'poly.f90', 'reconstruct.f90','flux1.f90']] config.add_extension('sharpclaw1', diff --git a/src/pyclaw/meson.build b/src/pyclaw/meson.build index 6c443052..86c8e6a7 100644 --- a/src/pyclaw/meson.build +++ b/src/pyclaw/meson.build @@ -236,11 +236,11 @@ fortran_extensions = { 'classic2_sw_sphere': ['limiter.f90','philim.f90','flux2.f90','step2ds.f90','../../../examples/shallow_sphere/step2qcor.f90','../../../examples/shallow_sphere/qcor.f90'], }, 'sharpclaw': { - 'sharpclaw1': ['ClawParams.f90','weno.f90','reconstruct.f90','evec.f90','workspace.f90','flux1.f90'], - 'sharpclaw2': ['ClawParams.f90','weno.f90','reconstruct.f90','evec.f90','workspace.f90','flux2.f90','flux1.f90'], - 'sharpclaw3': ['ClawParams.f90','weno.f90','reconstruct.f90','evec.f90','workspace.f90','flux3.f90','flux1.f90'], + 'sharpclaw1': ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90','evec.f90','workspace.f90','flux1.f90'], + 'sharpclaw2': ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90','evec.f90','workspace.f90','flux2.f90','flux1.f90'], + 'sharpclaw3': ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90','evec.f90','workspace.f90','flux3.f90','flux1.f90'], # 1D Euler with characteristic decomposition for SharpClaw limiting - 'euler_sharpclaw1': ['ClawParams.f90','weno.f90','reconstruct.f90','../../../examples/euler_1d/evec.f90','workspace.f90','flux1.f90'], + 'euler_sharpclaw1': ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90','../../../examples/euler_1d/evec.f90','workspace.f90','flux1.f90'], 'euler_tfluct1': ['../../../examples/euler_1d/euler_tfluct.f90'], }, } diff --git a/src/pyclaw/sharpclaw/ClawParams.f90 b/src/pyclaw/sharpclaw/ClawParams.f90 index b1e28357..75f67171 100644 --- a/src/pyclaw/sharpclaw/ClawParams.f90 +++ b/src/pyclaw/sharpclaw/ClawParams.f90 @@ -9,7 +9,7 @@ module ClawParams integer :: num_dim, num_waves, index_capa ! Method-related parameters: - integer :: char_decomp, lim_type, weno_order + integer :: char_decomp, lim_type, reconstruction_order integer, allocatable :: mthlim(:) logical :: fwave, tfluct_solver diff --git a/src/pyclaw/sharpclaw/flux1.f90 b/src/pyclaw/sharpclaw/flux1.f90 index e631adaa..689ac85a 100644 --- a/src/pyclaw/sharpclaw/flux1.f90 +++ b/src/pyclaw/sharpclaw/flux1.f90 @@ -115,9 +115,11 @@ subroutine flux1(q1d,dq1d,aux,dt,cfl,t,ixyz,num_aux,num_eqn,mx,num_ghost,maxnx,r end select case(3) call weno5(q1d,ql,qr,num_eqn,maxnx,num_ghost) + case(4:6) ! Linear polynomial reconstruction (module poly) + call poly_comp(q1d,ql,qr,num_eqn,maxnx,num_ghost) case default write(*,*) 'ERROR: Unrecognized limiter type option' - write(*,*) 'You should set 1<=lim_type<=3' + write(*,*) 'You should set 1<=lim_type<=6' stop end select diff --git a/src/pyclaw/sharpclaw/poly.f90 b/src/pyclaw/sharpclaw/poly.f90 new file mode 100644 index 00000000..72231682 --- /dev/null +++ b/src/pyclaw/sharpclaw/poly.f90 @@ -0,0 +1,154 @@ +! =================================================================== +module poly +! =================================================================== +! Linear (unlimited) polynomial reconstruction of interface values +! from cell averages, for the SharpClaw finite-volume solvers. +! +! q(:,i) holds the average of q over cell i. For each cell, ql(:,i) is +! the reconstructed value at the left face x_(i-1/2) and qr(:,i) the +! value at the right face x_(i+1/2). The weights are exact: they are +! obtained by differentiating the polynomial that interpolates the +! primitive function of q at the cell faces, so they are correct for +! cell averages (not point values). They are applied component-wise. +! +! lim_type selects the stencil; order = solver.poly_order, +! an even number N in {4, 6, 8, 10}. With k = N/2: +! +! 4 ('poly') ql(i) uses cells i-k..i+k-1 (N cells, +! centred on the face). Order N. ql(i+1) +! equals qr(i), so the scheme adds no +! numerical dissipation. +! 5 ('poly_dissipative') ql(i) uses cells i-k+1..i+k (the central +! stencil shifted one cell upwind for each +! face state). Order N, dissipative. +! 6 ('poly_cell_centered') ql(i) uses cells i-k..i+k (N+1 cells, +! centred on the cell, polynomial degree N). +! Order N+1, dissipative; equal to WENO of +! order N+1 with its ideal linear weights. +! +! qr uses the mirror image of the ql stencil. Every variant reaches k +! cells from the cell being reconstructed, so num_ghost = k + 1. +! +! Coefficients generated with exact rational arithmetic. +! =================================================================== + + implicit none + +contains + + ! =================================================================== + subroutine poly_weights(lim_type, order, lo, hi, w, ok) + ! =================================================================== + ! Weights for the left-face value ql(i) = sum_j w(j) * q(i+j), + ! j = lo..hi. + + integer, intent(in) :: lim_type, order + integer, intent(out) :: lo, hi + double precision, intent(out) :: w(-5:5) + logical, intent(out) :: ok + + ok = .true. + lo = 0; hi = 0 + w = 0.d0 + + select case (lim_type) + case (4) ! central, non-dissipative + select case (order) + case (4) + lo = -2; hi = 1 + w(-2:1) = [-1.d0, 7.d0, 7.d0, -1.d0] / 12.d0 + case (6) + lo = -3; hi = 2 + w(-3:2) = [1.d0, -8.d0, 37.d0, 37.d0, -8.d0, 1.d0] / 60.d0 + case (8) + lo = -4; hi = 3 + w(-4:3) = [-3.d0, 29.d0, -139.d0, 533.d0, 533.d0, -139.d0, 29.d0, -3.d0] / 840.d0 + case (10) + lo = -5; hi = 4 + w(-5:4) = [2.d0, -23.d0, 127.d0, -473.d0, 1627.d0, & + 1627.d0, -473.d0, 127.d0, -23.d0, 2.d0] / 2520.d0 + case default + ok = .false. + end select + case (5) ! upwind-biased, dissipative + select case (order) + case (4) + lo = -1; hi = 2 + w(-1:2) = [3.d0, 13.d0, -5.d0, 1.d0] / 12.d0 + case (6) + lo = -2; hi = 3 + w(-2:3) = [-2.d0, 22.d0, 57.d0, -23.d0, 7.d0, -1.d0] / 60.d0 + case (8) + lo = -3; hi = 4 + w(-3:4) = [5.d0, -55.d0, 365.d0, 743.d0, -307.d0, 113.d0, -27.d0, 3.d0] / 840.d0 + case (10) + lo = -4; hi = 5 + w(-4:5) = [-3.d0, 37.d0, -233.d0, 1207.d0, 2131.d0, & + -893.d0, 367.d0, -113.d0, 22.d0, -2.d0] / 2520.d0 + case default + ok = .false. + end select + case (6) ! cell-centred, dissipative + select case (order) + case (4) + lo = -2; hi = 2 + w(-2:2) = [-3.d0, 27.d0, 47.d0, -13.d0, 2.d0] / 60.d0 + case (6) + lo = -3; hi = 3 + w(-3:3) = [4.d0, -38.d0, 214.d0, 319.d0, -101.d0, 25.d0, -3.d0] / 420.d0 + case (8) + lo = -4; hi = 4 + w(-4:4) = [-5.d0, 55.d0, -305.d0, 1375.d0, & + 1879.d0, -641.d0, 199.d0, -41.d0, 4.d0] / 2520.d0 + case (10) + lo = -5; hi = 5 + w(-5:5) = [12.d0, -153.d0, 947.d0, -4003.d0, 15797.d0, & + 20417.d0, -7303.d0, 2597.d0, -703.d0, 122.d0, -10.d0] / 27720.d0 + case default + ok = .false. + end select + case default + ok = .false. + end select + + end subroutine poly_weights + + ! =================================================================== + subroutine poly_recon(q,ql,qr,num_eqn,maxnx,num_ghost,lim_type,order) + ! =================================================================== + + integer, intent(in) :: num_eqn, maxnx, num_ghost, lim_type, order + double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost) + double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost) + + double precision :: w(-5:5) + integer :: lo, hi, reach, i, j, mx2 + logical :: ok + + call poly_weights(lim_type, order, lo, hi, w, ok) + if (.not. ok) then + print *, 'ERROR: polynomial reconstruction needs lim_type 4, 5 or 6' + print *, ' and an even reconstruction_order between 4 and 10.' + error stop + endif + + mx2 = size(q,2) + reach = max(-lo, hi) + if (reach >= num_ghost) then + print *, 'ERROR: polynomial reconstruction needs num_ghost >= ', reach+1 + error stop + endif + + ql = 0.d0 + qr = 0.d0 + do i = 1+reach, mx2-reach + do j = lo, hi + ql(:,i) = ql(:,i) + w(j)*q(:,i+j) + ! qr mirrors the ql stencil about cell i + qr(:,i) = qr(:,i) + w(j)*q(:,i-j) + end do + end do + + end subroutine poly_recon + +end module poly diff --git a/src/pyclaw/sharpclaw/reconstruct.f90 b/src/pyclaw/sharpclaw/reconstruct.f90 index 4db69532..fb75d099 100644 --- a/src/pyclaw/sharpclaw/reconstruct.f90 +++ b/src/pyclaw/sharpclaw/reconstruct.f90 @@ -83,8 +83,12 @@ subroutine dealloc_recon_workspace(lim_type,char_decomp) deallocate(hh) deallocate(uh) end select - recon_alloc = .False. + case(3) + deallocate(uu) + deallocate(dq1m) end select + recon_alloc = .False. + end subroutine dealloc_recon_workspace ! =================================================================== @@ -97,14 +101,14 @@ subroutine weno_comp(q,ql,qr,num_eqn,maxnx,num_ghost) ! It does no characteristic decomposition use weno - use clawparams, only: weno_order + use clawparams, only: reconstruction_order implicit none integer, intent(in) :: num_eqn, maxnx, num_ghost double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost) double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost) - select case(weno_order) + select case(reconstruction_order) case (5) call weno5(q,ql,qr,num_eqn,maxnx,num_ghost) case (7) @@ -120,7 +124,7 @@ subroutine weno_comp(q,ql,qr,num_eqn,maxnx,num_ghost) case (17) call weno17(q,ql,qr,num_eqn,maxnx,num_ghost) case default - print *, 'ERROR: weno_order must be an odd number between 5 and 17 (inclusive).' + print *, 'ERROR: reconstruction_order in WENO must be an odd number between 5 and 17 (inclusive).' stop end select @@ -821,4 +825,22 @@ subroutine tvd2_wave(q,ql,qr,wave,s,mthlim,num_eqn,num_ghost) return end subroutine tvd2_wave + ! =================================================================== + subroutine poly_comp(q,ql,qr,num_eqn,maxnx,num_ghost) + ! =================================================================== + ! Linear polynomial reconstruction; see module poly for the + ! meaning of lim_type 4, 5 and 6. + + use poly, only: poly_recon + use clawparams, only: lim_type, reconstruction_order + implicit none + + integer, intent(in) :: num_eqn, maxnx, num_ghost + double precision, intent(in) :: q(num_eqn,maxnx+2*num_ghost) + double precision, intent(out) :: ql(num_eqn,maxnx+2*num_ghost),qr(num_eqn,maxnx+2*num_ghost) + + call poly_recon(q,ql,qr,num_eqn,maxnx,num_ghost,lim_type,reconstruction_order) + + end subroutine poly_comp + end module reconstruct diff --git a/src/pyclaw/sharpclaw/setup.py b/src/pyclaw/sharpclaw/setup.py index ece05310..404afec2 100644 --- a/src/pyclaw/sharpclaw/setup.py +++ b/src/pyclaw/sharpclaw/setup.py @@ -5,16 +5,16 @@ def configuration(parent_package='',top_path=None): config = Configuration('sharpclaw', parent_package, top_path) config.add_extension('sharpclaw1', - ['ClawParams.f90','weno.f90','reconstruct.f90', + ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90', 'evec.f90','workspace.f90','flux1.f90']) config.add_extension('sharpclaw2', - ['ClawParams.f90','weno.f90','reconstruct.f90', + ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90', 'evec.f90','workspace.f90','flux2.f90', 'flux1.f90']) config.add_extension('sharpclaw3', - ['ClawParams.f90','weno.f90','reconstruct.f90', + ['ClawParams.f90','weno.f90','poly.f90','reconstruct.f90', 'evec.f90','workspace.f90','flux3.f90', 'flux1.f90']) diff --git a/src/pyclaw/sharpclaw/solver.py b/src/pyclaw/sharpclaw/solver.py index 7f991e6d..762bbefb 100644 --- a/src/pyclaw/sharpclaw/solver.py +++ b/src/pyclaw/sharpclaw/solver.py @@ -5,6 +5,7 @@ # Created: 2010-03-20 # Author: David Ketcheson """ + # Solver superclass from clawpack.pyclaw.solver import Solver from clawpack.pyclaw.util import add_parent_doc @@ -36,18 +37,58 @@ class SharpClawSolver(Solver): .. attribute:: lim_type - Limiter(s) to be used. - - 0: No limiting. + Reconstruction used to compute interface values from cell + averages. Can be set with an integer or, equivalently, a name: + + - -1: first-order Godunov (Python kernels only). + - 0: no limiting (not implemented). - 1: TVD reconstruction. - - 2: WENO reconstruction. + - 2: WENO reconstruction (odd + ``weno_order``, 5 to 17). + - 3: original fifth-order WENO routine. + - 4 or ``'poly'``: non-dissipative polynomial reconstruction. + - 5 or ``'poly_dissipative'``: dissipative (upwind-biased) + polynomial reconstruction. + - 6 or ``'poly_cell_centered'``: dissipative, cell-centred + polynomial reconstruction. + + The value is always stored, and read back, as an integer. + + The polynomial reconstructions (4, 5 and 6) are linear and + unlimited, so they suit smooth solutions, e.g. electromagnetic + waves in smooth media; at discontinuities they produce Gibbs + oscillations. They take an even ``poly_order`` N in + {4, 6, 8, 10}, are available with Fortran kernels only, and act + component-wise (``char_decomp`` is ignored). With k = N/2: + + - ``'poly'`` uses the N cells centred on each face, so both + sides of a face get the same value: no numerical dissipation + is added (only the time integrator dissipates); for linear, + constant-coefficient problems the semi-discrete scheme + conserves the discrete L2 energy. Formal order N. Unresolved (grid-scale) + modes are not damped. + - ``'poly_dissipative'`` uses the same width, shifted one cell + upwind for each face state. Formal order N, with damping + that grows with wavenumber. + - ``'poly_cell_centered'`` uses the N+1 cells centred on each + cell (a polynomial of degree N). Formal order N+1, and less + dissipative than ``'poly_dissipative'``; it equals WENO of + order N+1 with its ideal linear weights. + + All three need k+1 ghost cells, which is set automatically. The + time integrator ``'SSP104'`` is fourth order, so with N > 4 the + time error dominates unless the CFL number is reduced. ``Default = 2`` .. attribute:: weno_order - Order of the WENO reconstruction. From 1st to 17th order (PyWENO) + Order of WENO reconstruction for ``lim_type`` 2. Default is 5. + + .. attribute:: poly_order - ``Default = 5`` + Polynomial stencil parameter N (4, 6, 8 or 10). Default is 4. + The formal order is N for types 4 and 5, and N+1 for type 6. .. attribute:: time_integrator @@ -150,6 +191,33 @@ class SharpClawSolver(Solver): 'SSPLMM53': [0.24, 0.25] } + _lim_type_names = { + 'poly': 4, + 'poly_dissipative': 5, + 'poly_cell_centered': 6, + } + + _poly_lim_types = (4, 5, 6) + _poly_orders = (4, 6, 8, 10) + + @property + def lim_type(self): + return self._lim_type + + @lim_type.setter + def lim_type(self, value): + if isinstance(value, str): + try: + value = self._lim_type_names[value.lower()] + except KeyError: + raise ValueError("Unknown lim_type '%s'; valid names are %s" + % (value, ', '.join(self._lim_type_names))) + self._lim_type = value + + _weno_orders = tuple(range(5, 18, 2)) + # lim_types the Fortran kernels (flux1.f90) implement + _fortran_lim_types = (1, 2, 3, 4, 5, 6) + # ======================================================================== # Initialization routines # ======================================================================== @@ -160,6 +228,7 @@ def __init__(self,riemann_solver=None,claw_package=None): self.limiters = [1] self.lim_type = 2 self.weno_order = 5 + self.poly_order = 4 self.time_integrator = 'SSP104' self.char_decomp = 0 self.tfluct_solver = False @@ -203,8 +272,32 @@ def setup(self,solution): """ Allocate RK stage arrays or previous step solutions and fortran routine work arrays. """ + # Validate here: the Fortran kernels answer a bad lim_type or order + # with a bare 'stop', which ends the process without an exception. + if self.lim_type not in (-1, 0) + self._fortran_lim_types: + raise ValueError('lim_type must be between -1 and 6, or one of the names %s; got %s.' + % (', '.join(self._lim_type_names), self.lim_type)) + if self.lim_type == 0: + raise NotImplementedError('Unlimited reconstruction (lim_type 0) is not implemented.') + if self.kernel_language == 'Fortran' and self.lim_type not in self._fortran_lim_types: + raise NotImplementedError('lim_type %s is implemented in Python kernels only.' + % self.lim_type) + if self.lim_type == 2: - self.num_ghost = (self.weno_order+1)//2 + if self.weno_order not in self._weno_orders: + raise ValueError('WENO (lim_type 2) needs an odd weno_order from 5 to 17, ' + 'not %s.' % self.weno_order) + self.num_ghost = int((self.weno_order+1)/2) + + if self.lim_type in self._poly_lim_types: + if self.poly_order not in self._poly_orders: + raise ValueError('Polynomial reconstruction (lim_type %s) needs an even ' + 'poly_order in %s, not %s.' + % (self.lim_type, self._poly_orders, self.poly_order)) + if self.kernel_language != 'Fortran': + raise NotImplementedError('Polynomial reconstruction (lim_type %s) is only ' + 'implemented in Fortran kernels.' % self.lim_type) + self.num_ghost = int(self.poly_order)//2 + 1 if self.lim_type == 2 and self.weno_order != 5 and self.kernel_language == 'Python': raise Exception('Only 5th-order WENO reconstruction is implemented in Python kernels. \ @@ -270,7 +363,7 @@ def __del__(self): Deallocate F90 module arrays. Also delete Fortran objects, which otherwise tend to persist in Python sessions. """ - if self.kernel_language=='Fortran': + if self.kernel_language == 'Fortran' and self.fmod is not None: self.fmod.clawparams.dealloc_clawparams() self.fmod.workspace.dealloc_workspace(self.char_decomp) self.fmod.reconstruct.dealloc_recon_workspace(self.fmod.clawparams.lim_type,self.fmod.clawparams.char_decomp) @@ -568,7 +661,8 @@ def _set_fortran_parameters(self,state,clawparams,workspace,reconstruct): grid = state.grid clawparams.num_dim = grid.num_dim clawparams.lim_type = self.lim_type - clawparams.weno_order = self.weno_order + clawparams.reconstruction_order = (self.poly_order if self.lim_type in self._poly_lim_types + else self.weno_order) clawparams.char_decomp = self.char_decomp clawparams.tfluct_solver = self.tfluct_solver clawparams.fwave = self.fwave diff --git a/src/pyclaw/tests/test_sharpclaw_poly.py b/src/pyclaw/tests/test_sharpclaw_poly.py new file mode 100644 index 00000000..2da9ff9a --- /dev/null +++ b/src/pyclaw/tests/test_sharpclaw_poly.py @@ -0,0 +1,278 @@ +""" +Tests for the SharpClaw polynomial reconstructions (lim_type 4, 5 and 6), +including independent WENO and polynomial order settings. +""" +import warnings + +import numpy as np +import pytest + +from clawpack import pyclaw, riemann + +sharpclaw1 = pytest.importorskip("clawpack.pyclaw.sharpclaw.sharpclaw1") + +POLY = {4: 'poly', 5: 'poly_dissipative', 6: 'poly_cell_centered'} +ORDERS = (4, 6, 8, 10) + + +def formal_order(lim_type, order): + """Formal order of accuracy of the interface values.""" + return order + 1 if lim_type == 6 else order + + +def sine_averages(xc, h, kwave=1): + """Exact cell averages of sin(2*pi*kwave*x).""" + w = 2*np.pi*kwave + return (np.cos(w*(xc - h/2)) - np.cos(w*(xc + h/2)))/(w*h) + + +def reconstruct(num_cells, lim_type, order): + """Reconstruct interface values of sin(2*pi*x) on a periodic grid; + return the max interface error and the max jump ql(i+1) - qr(i).""" + num_ghost = order//2 + 1 + h = 1./num_cells + xc = (np.arange(-num_ghost, num_cells + num_ghost) + 0.5)*h + q = sine_averages(xc, h)[None, :] + ql, qr = sharpclaw1.poly.poly_recon(q, num_cells, num_ghost, lim_type, order) + real = slice(num_ghost, num_cells + num_ghost) + err_l = np.abs(ql[0, real] - np.sin(2*np.pi*(xc[real] - h/2))).max() + err_r = np.abs(qr[0, real] - np.sin(2*np.pi*(xc[real] + h/2))).max() + jump = np.abs(ql[0, num_ghost+1:num_cells+num_ghost+1] + - qr[0, num_ghost:num_cells+num_ghost]).max() + return max(err_l, err_r), jump + + +def advect(num_cells, lim_type, order, kwave=1, cfl=0.4, kernel_language='Fortran'): + """Advect sin(2*pi*kwave*x) once around a periodic unit domain; + return the L1 error and the ratio of final to initial L2 norm.""" + solver = pyclaw.SharpClawSolver1D(riemann.advection_1D) + solver.kernel_language = kernel_language + solver.lim_type = lim_type + solver.poly_order = order + solver.bc_lower[0] = pyclaw.BC.periodic + solver.bc_upper[0] = pyclaw.BC.periodic + solver.cfl_desired = cfl + solver.cfl_max = 1.05*cfl + + x = pyclaw.Dimension(0., 1., num_cells, name='x') + domain = pyclaw.Domain(x) + state = pyclaw.State(domain, 1) + state.problem_data['u'] = 1. + state.q[0, :] = sine_averages(x.centers, x.delta, kwave) + + claw = pyclaw.Controller() + claw.solution = pyclaw.Solution(state, domain) + claw.solver = solver + claw.tfinal = 1. + claw.num_output_times = 1 + claw.output_format = None + claw.keep_copy = True + claw.verbosity = 0 + claw.run() + + q0 = claw.frames[0].q[0] + qf = claw.frames[-1].q[0] + error = x.delta*np.abs(qf - q0).sum() + energy_ratio = np.linalg.norm(qf)/np.linalg.norm(q0) + return error, energy_ratio + + +class TestReconstruction: + + @pytest.mark.parametrize("lim_type", sorted(POLY)) + @pytest.mark.parametrize("order", ORDERS) + def test_order_of_accuracy(self, lim_type, order): + # Coarse grids for the high orders keep the errors above round-off + grids = (16, 32) if order < 8 else (8, 16) + coarse, fine = (reconstruct(n, lim_type, order)[0] for n in grids) + observed = np.log2(coarse/fine) + assert observed > formal_order(lim_type, order) - 0.4 + + @pytest.mark.parametrize("order", ORDERS) + def test_poly_is_single_valued_at_faces(self, order): + """lim_type 4 gives the same value on both sides of each face, + so the Riemann solver adds no dissipation.""" + assert reconstruct(32, 4, order)[1] < 1e-14 + + @pytest.mark.parametrize("lim_type", (5, 6)) + def test_dissipative_variants_jump_at_faces(self, lim_type): + assert reconstruct(32, lim_type, 4)[1] > 1e-6 + + +class TestSolver: + + @pytest.mark.parametrize("lim_type", sorted(POLY)) + def test_convergence(self, lim_type): + coarse, fine = (advect(n, lim_type, 4)[0] for n in (40, 80)) + assert np.log2(coarse/fine) > formal_order(lim_type, 4) - 0.3 + + @pytest.mark.parametrize("lim_type", sorted(POLY)) + @pytest.mark.parametrize("order", (8, 10)) + def test_high_orders_are_stable(self, lim_type, order): + error, energy_ratio = advect(40, lim_type, order) + assert error < 1e-5 + assert energy_ratio <= 1. + 1e-12 + + def test_poly_conserves_energy(self): + """An unresolved mode (4 cells per wavelength) keeps its energy + with 'poly' but is damped by the dissipative variants.""" + ratios = {name: advect(16, name, 6, kwave=4)[1] for name in POLY.values()} + assert ratios['poly'] > 0.999 + assert ratios['poly_dissipative'] < ratios['poly_cell_centered'] < 0.5 + + @pytest.mark.parametrize("order", (5, 7, 12)) + def test_rejects_invalid_order(self, order): + with pytest.raises(ValueError): + advect(20, 'poly', order) + + def test_rejects_python_kernels(self): + with pytest.raises(NotImplementedError): + advect(20, 'poly', 4, kernel_language='Python') + + +class TestSetupValidation: + """setup() rejects what the Fortran kernels would answer with a bare + 'stop', which ends the process with exit status 0 and no exception.""" + + def setup_solver(self, lim_type, order, kernel_language='Fortran'): + solver = pyclaw.SharpClawSolver1D(riemann.advection_1D) + solver.kernel_language = kernel_language + solver.lim_type = lim_type + solver.weno_order = order + solver.bc_lower[0] = pyclaw.BC.periodic + solver.bc_upper[0] = pyclaw.BC.periodic + domain = pyclaw.Domain(pyclaw.Dimension(0., 1., 20, name='x')) + state = pyclaw.State(domain, 1) + state.problem_data['u'] = 1. + solver.setup(pyclaw.Solution(state, domain)) + return solver + + @pytest.mark.parametrize("order", (4, 6, 19)) + def test_weno_rejects_even_or_out_of_range_order(self, order): + with pytest.raises(ValueError): + self.setup_solver(2, order) + + @pytest.mark.parametrize("order", (5, 17)) + def test_weno_accepts_odd_orders(self, order): + assert self.setup_solver(2, order).num_ghost == (order + 1)//2 + + def test_godunov_is_python_only(self): + with pytest.raises(NotImplementedError): + self.setup_solver(-1, 5) + assert self.setup_solver(-1, 5, kernel_language='Python').lim_type == -1 + + def test_unlimited_is_not_implemented(self): + with pytest.raises(NotImplementedError): + self.setup_solver(0, 5, kernel_language='Python') + + @pytest.mark.parametrize("lim_type", (-2, 7, 1000)) + def test_unknown_lim_type(self, lim_type): + with pytest.raises(ValueError): + self.setup_solver(lim_type, 5) + +class TestLimTypeNames: + + @pytest.mark.parametrize("name, value", [ + ('poly', 4), ('poly_dissipative', 5), ('poly_cell_centered', 6)]) + def test_name_maps_to_integer(self, name, value): + solver = pyclaw.SharpClawSolver1D(riemann.advection_1D) + solver.lim_type = name + assert solver.lim_type == value + solver.lim_type = name.upper() + assert solver.lim_type == value + solver.lim_type = value + assert solver.lim_type == value + + def test_unknown_name_raises(self): + solver = pyclaw.SharpClawSolver1D(riemann.advection_1D) + with pytest.raises(ValueError): + solver.lim_type = 'spline' + + def test_poly_name_and_integer_give_same_result(self): + assert advect(20, 'poly', 6) == advect(20, 4, 6) + + +@pytest.mark.parametrize("select_poly_first", [True, False]) +def test_independent_orders(select_poly_first): + solver = pyclaw.SharpClawSolver1D(riemann.advection_1D) + assert (solver.lim_type, solver.weno_order, solver.poly_order) == (2, 5, 4) + with warnings.catch_warnings(): + warnings.simplefilter("error", FutureWarning) + if select_poly_first: + solver.lim_type = 'poly' + solver.weno_order = 11 + solver.poly_order = 6 + assert solver.lim_type == (4 if select_poly_first else 2) + solver.lim_type = 'poly' + assert (solver.weno_order, solver.poly_order) == (11, 6) + solver.lim_type = 2 + assert (solver.weno_order, solver.poly_order) == (11, 6) + + +def test_poly_2d(): + """The 2D solver reuses the 1D reconstruction along each direction.""" + solver = pyclaw.SharpClawSolver2D(riemann.advection_2D) + solver.lim_type = 'poly' + solver.poly_order = 6 + solver.bc_lower = [pyclaw.BC.periodic]*2 + solver.bc_upper = [pyclaw.BC.periodic]*2 + solver.cfl_desired = 0.4 + solver.cfl_max = 0.42 + + x = pyclaw.Dimension(0., 1., 24, name='x') + y = pyclaw.Dimension(0., 1., 24, name='y') + domain = pyclaw.Domain([x, y]) + state = pyclaw.State(domain, 1) + state.problem_data['u'] = 1. + state.problem_data['v'] = 1. + X, Y = state.grid.p_centers + h = x.delta + state.q[0, ...] = sine_averages(X, h)*sine_averages(Y, h) + + claw = pyclaw.Controller() + claw.solution = pyclaw.Solution(state, domain) + claw.solver = solver + claw.tfinal = 1. + claw.num_output_times = 1 + claw.output_format = None + claw.keep_copy = True + claw.verbosity = 0 + claw.run() + + assert solver.num_ghost == 4 + q0 = claw.frames[0].q[0] + qf = claw.frames[-1].q[0] + assert np.abs(qf - q0).max() < 1e-4 + + +@pytest.fixture(scope='module') +def generator(): + """development/poly_reconstruction/poly_weights.py, loaded by path.""" + import importlib.util + import os + path = os.path.join(os.path.dirname(__file__), '..', '..', '..', 'development', + 'poly_reconstruction', 'poly_weights.py') + if not os.path.exists(path): + pytest.skip('development/ is not available') + spec = importlib.util.spec_from_file_location('poly_weights', path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +class TestGenerator: + """poly.f90 is generated by development/poly_reconstruction/poly_weights.py; + these tests stop the two from drifting apart.""" + + def test_poly_f90_matches_generator(self, generator): + with open(generator.POLY_F90) as f: + assert f.read() == generator.fortran_source(), \ + 'poly.f90 is out of date: run development/poly_reconstruction/poly_weights.py --write' + + @pytest.mark.parametrize("lim_type", sorted(POLY)) + @pytest.mark.parametrize("order", ORDERS) + def test_compiled_weights_match_exact_weights(self, generator, lim_type, order): + lo, hi, w, ok = sharpclaw1.poly.poly_weights(lim_type, order) + cells, exact = generator.ql_weights(lim_type, order) + assert ok and (lo, hi) == (cells[0], cells[-1]) + np.testing.assert_allclose(w[lo+5:hi+6], [float(v) for v in exact], rtol=0, atol=1e-15)