def execute_lambdify(ui, spacing=0.01, a=0.5, timesteps=500):
"""Execute diffusion stencil using vectorised numpy array accesses."""
nx, ny = ui.shape
dx2, dy2 = spacing**2, spacing**2
dt = dx2 * dy2 / (2 * a * (dx2 + dy2))
u = np.concatenate((ui, np.zeros_like(ui))).reshape((2, nx, ny))
def diffusion_stencil():
"""Create stencil and substitutions for the diffusion equation"""
p = sympy.Function('p')
x, y, t, h, s = sympy.symbols('x y t h s')
dx2 = p(x, y, t).diff(x, x).as_finite_difference([x - h, x, x + h])
dy2 = p(x, y, t).diff(y, y).as_finite_difference([y - h, y, y + h])
dt = p(x, y, t).diff(t).as_finite_difference([t, t + s])
eqn = sympy.Eq(dt, a * (dx2 + dy2))
stencil = sympy.solve(eqn, p(x, y, t + s))[0]
return stencil, (p(x, y, t), p(x + h, y, t), p(x - h, y, t),
p(x, y + h, t), p(x, y - h, t), s, h)
stencil, subs = diffusion_stencil()
kernel = sympy.lambdify(subs, stencil, 'numpy')
# Execute timestepping loop with alternating buffers
tstart = time.time()
for ti in range(timesteps):
t0 = ti % 2
t1 = (ti + 1) % 2
u[t1, 1:-1, 1:-1] = kernel(u[t0, 1:-1, 1:-1], u[t0, 2:, 1:-1],
u[t0, :-2, 1:-1], u[t0, 1:-1, 2:],
u[t0, 1:-1, :-2], dt, spacing)
runtime = time.time() - tstart
log("Lambdify: Diffusion with dx=%0.4f, dy=%0.4f, executed %d timesteps in %f seconds"
% (spacing, spacing, timesteps, runtime))
return u[ti % 2, :, :], runtime
评论列表
文章目录
