import numpy as np import numpy.linalg as la # noqa: F401 import pyopencl as cl # noqa: F401 import pyopencl.array # noqa import pyopencl.tools # noqa import pyopencl.clrandom # noqa import loopy as lp # noqa ndim = 3 nvars = 3 gamma = 1.4 def rho(state): return state[0] def velocity(state): return state[1:ndim+1]/rho(state) def energy(state): return state[ndim+1] def kinetic_energy(state): vel = velocity(state) return 0.5*rho(state)*np.dot(vel, vel) def internal_energy(state): return energy(state) - kinetic_energy(state) def pressure(state): return (gamma - 1)*internal_energy(state) def sound_speed(state): return np.sqrt(gamma*pressure(state)/rho(state)) def flux(state): vel = velocity(state) p = pressure(state) result = np.outer(state, vel) for i in range(ndim): result[i+1,i] += p result[ndim+1,i] += vel[i]*p return result