Convergence of the compressible solvers

We’ll look at convergence of the 2nd order compressible and 4th order compressible_fv4 solvers using the acoustic_pulse problem and doing simple Richardson convergence testing.

from pyro import Pyro

We want to keep \(\Delta t / \Delta x\) constant as we test convergence so we will use a fixed timestep, following:

\[\Delta t = 3\times 10^{-3} \frac{64}{N}\]

where \(N\) is the number of zones in a dimension.

def timestep(N):
    return 3.e-3 * 64.0 / N

compressible

We’ll run the problem at several different resolutions and store the Pyro simulation objects in a list.

sims = []

for N in [32, 64, 128, 256]:
    dt = timestep(N)
    params = {"driver.fix_dt": dt, "mesh.nx": N, "mesh.ny": N}
    p = Pyro("compressible")
    p.initialize_problem(problem_name="acoustic_pulse", inputs_dict=params)
    p.run_sim()
    sims.append(p)

Now we want to loop over each adjacent pair of simulations, coarsen the finer resolution simulation and compute the norm of the difference. We’ll do this for a single variable.

from itertools import pairwise
var = "density"

for coarse, fine in pairwise(sims):
    cvar = coarse.get_var(var)
    fvar = fine.sim.cc_data.restrict(var)
    e = cvar - fvar
    print(f"{fine.get_grid().nx:3} -> {coarse.get_grid().nx:3} : {e.norm()}")

 64 ->  32 : 0.0002674195946900653
128 ->  64 : 5.7696409241208797e-05
256 -> 128 : 1.3860268814816614e-05

We see that the error is dropping by a factor of ~4 each time, indicating 2nd order convergence.

compressible_fv4

Now we’ll do the same for the 4th order solver. We need to change the Riemann solver to

sims = []

for N in [32, 64, 128, 256]:
    dt = timestep(N)
    params = {"driver.fix_dt": dt, "mesh.nx": N, "mesh.ny": N}
    p = Pyro("compressible_fv4")
    p.initialize_problem(problem_name="acoustic_pulse", inputs_dict=params)
    p.run_sim()
    sims.append(p)

for coarse, fine in pairwise(sims):
    cvar = coarse.get_var(var)
    fvar = fine.sim.cc_data.restrict(var)
    e = cvar - fvar
    print(f"{fine.get_grid().nx:3} -> {coarse.get_grid().nx:3} : {e.norm()}")

 64 ->  32 : 6.519131423273572e-05
128 ->  64 : 4.825569192556014e-06
256 -> 128 : 3.0769222917915304e-07

Now we see that the convergence is close to 4th order, with the error decreasing close to a factor of 16.