SHARE
TWEET

Gray Scott reaction-diffusion

a guest Sep 2nd, 2019 187 in 267 days
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #!/usr/bin/env python3
  2.  
  3. import numba
  4. import numpy as np
  5.  
  6. @numba.stencil
  7. def laplace(arr):
  8.     return -arr[0,0] + 0.25 * (arr[-1,0] + arr[0,-1] + arr[1,0] + arr[0,1])
  9.  
  10. @numba.njit(parallel=True)
  11. def grayscott(niters, u, v, du, dv, f, k, dt):
  12.     for i in range(niters):
  13.         conv = u * v * v
  14.         u += dt * (du * laplace(u) - conv + f * (1 - u))
  15.         v += dt * (dv * laplace(v) + conv - (f + k) * v)
  16.     return u, v
  17.  
  18.  
  19. if __name__ == '__main__':
  20.     import time
  21.  
  22.     m, n = 320, 320  # Grid shape
  23.     u = np.ones((m, n))
  24.     v = np.zeros((m, n))
  25.  
  26.     du = 2.0  # diffusion coefficient of u
  27.     dv = 1.0  # diffusion coefficient of v
  28.     f = 0.037  # feed rate
  29.     k = 0.061  # kill rate
  30.     dt = 0.45  # discrete time step (dt * du should be less than 1.0, try lowering if things explode or diminish)
  31.  
  32.     # Initialize seed_num random points used to seed v
  33.     seed_num = 10
  34.     y = np.random.randint(0, m, seed_num)
  35.     x = np.random.randint(0, n, seed_num)
  36.  
  37.     # Seed v for a couple hundred iterations, updating v after every 5 iterations
  38.     for _ in range(40):
  39.         v[y, x] = 0.5 * (v[y, x] + 1)
  40.         u, v = grayscott(5, u, v, du, dv, f, k, dt)
  41.  
  42.     # On an old 4th gen i7, 5000 iterations takes ~10-15 seconds, and scales (roughly) linearly with grid size
  43.     tic = time.time()
  44.     u, v = grayscott(5000, u, v, du, dv, f, k, dt)
  45.     toc = time.time()
  46.     print(f'Done, took {toc - tic :.3f} s')
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top