richtmyer_meshkov_instability

import numpy
from numba import jit
import matplotlib.pyplot as plt
from tqdm import tqdm

@jit(nopython=True)
def calc_vel(p):

    x = p[0]
    y = p[1]
    vx = numpy.sign(y)*numpy.sin(x)*numpy.exp(-numpy.abs(y))
    vy = numpy.cos(x)*numpy.exp(-numpy.abs(y))
    return numpy.array([vx,vy])

@jit(nopython=True)
def calc_all_vel(p_list):

    res = numpy.zeros_like(p_list)
    for i in range(len(res)):
        res[i] = calc_vel(p_list[i])
    return res

x0_list = numpy.linspace(0,20,2000)
y0_list = 1e-3*numpy.sin(x0_list)
p_list = numpy.vstack((x0_list, y0_list)).T

fig, ax = plt.subplots()
ax.set_ylim((-1,1))
image, = ax.plot(x0_list, y0_list)
plt.draw()

dt = 0.2e-2
v_list = numpy.zeros_like(p_list)
counter = 0
for n in tqdm(range(800)):
    p_list = p_list + calc_all_vel(p_list)*dt
    if n%2 == 0:
        image.set_data(p_list.T[0],
                       p_list.T[1])
        plt.draw()
        fig.savefig('output/snapshot_'+str(counter).zfill(3)+'.png')
        counter += 1