Convergence of the Euler Solver

Test the convergence of the Euler solver by running the acoustic pulse problem and doing Richardson extrapolation to estimate the convergence.

We’ll use the acoustic pulse problem described in McCorquodale & Colella 2011.

from ppmpy import Euler
from ppmpy.initial_conditions import acoustic_pulse

import numpy as np

Here we run simulations at different resolutions and store the simulation object in a list.

simulations = []
for nx in [32, 64, 128, 256]:
    dt = 1.0 / nx * 0.192
    e = Euler(nx, 0.5, fixed_dt=dt, init_cond=acoustic_pulse,
              use_flattening=False, use_limiting=False,
              bc_left_type="periodic", bc_right_type="periodic")
    e.evolve(0.24, verbose=False)
    print(f"nx = {nx}, number of steps = {e.nstep}")
    simulations.append(e)

nx = 32, number of steps = 40

nx = 64, number of steps = 80

nx = 128, number of steps = 160

nx = 256, number of steps = 320

Now we’ll loop over pairs of simulations and coarsen the finer simulation by a factor of 2, so it is at the same resolution as the coarser simulation, and then compute the norm of the difference.

from itertools import pairwise
ivar = 0

for coarse, fine in pairwise(simulations):
    _, cd = fine.grid.coarsen(fine.U[:, ivar])
    err = coarse.grid.norm(coarse.U[:, ivar] - cd)
    print(f"{fine.grid.nx:3d} -> {coarse.grid.nx:3d} : {err}")

 64 ->  32 : 0.0001187035470491678
128 ->  64 : 2.622374113240233e-05
256 -> 128 : 6.732654910238016e-06

We see that we are converging second order as we increase the resolution.