[OpenGL] Basic OpenGL backend #492

cmake/TaichiCore.cmake

@@ -24,6 +24,7 @@ endif()
 
 option(USE_STDCPP "Use -stdlib=libc++" OFF)
 option(TI_WITH_CUDA "Build with CUDA support" OFF)
+option(TI_WITH_OPENGL "Build with OpenGL backend" ON)
 
 include_directories(${CMAKE_SOURCE_DIR})
 include_directories(external/xxhash)
@@ -53,6 +54,25 @@ if (TI_WITH_CUDA)
     endif()
 endif()
 
+if (TI_WITH_OPENGL)
+  if(NOT GLEW_VERSION)
+    set(GLEW_VERSION 2.0.0)
+  endif()
+  find_package(GLEW ${GLEW_VERSION})
+  if (GLEW_FOUND)
+    message("Building with GLEW ${GLEW_VERSION}")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DTI_WITH_OPENGL=1")
+    #if (MSVC)
+    # TODO(archibate): figure out how GL is installed on Windows...
+    #else()
+    target_link_libraries(${LIBRARY_NAME} /usr/lib/libGLEW.so GLEW)
+    target_link_libraries(${LIBRARY_NAME} /usr/lib/libglfw.so glfw)
+    #endif()
+  else()
+    message(FATAL_ERROR "GLEW not found")
+  endif()
+endif()
+
 # http://llvm.org/docs/CMake.html#embedding-llvm-in-your-project
 find_package(LLVM REQUIRED CONFIG 8.0)
 message(STATUS "Found LLVM ${LLVM_PACKAGE_VERSION}")

examples/opengl_backend.py

@@ -0,0 +1,17 @@
+import taichi as ti
+
+ti.init(arch=ti.opengl)
+
+x = ti.var(ti.f32, shape=(2, 4))
+
+@ti.kernel
+def func():
+  x[1, 2] = 666
+  x[2, 1] = 527
+
+x[1, 2] = 233
+x[1, 3] = 888
+func()
+print(x[1, 2])
+print(x[1, 3])
+print(x[2, 1])

examples/opengl_mpm99.py

@@ -0,0 +1,104 @@
+import taichi as ti
+import numpy as np
+from time import perf_counter
+ti.init(arch=ti.opengl) # Try to run on GPU
+quality = 1 # Use a larger value for higher-res simulations
+n_particles, n_grid = 9000 * quality ** 2, 128 * quality
+dx, inv_dx = 1 / n_grid, float(n_grid)
+dt = 1e-4 / quality
+p_vol, p_rho = (dx * 0.5)**2, 1
+p_mass = p_vol * p_rho
+E, nu = 0.1e4, 0.2 # Young's modulus and Poisson's ratio
+mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1+nu) * (1 - 2 * nu)) # Lame parameters
+
+x = ti.Vector(2, dt=ti.f32, shape=n_particles) # position
+v = ti.Vector(2, dt=ti.f32, shape=n_particles) # velocity
+C = ti.Matrix(2, 2, dt=ti.f32, shape=n_particles) # affine velocity field
+F = ti.Matrix(2, 2, dt=ti.f32, shape=n_particles) # deformation gradient
+material = ti.var(dt=ti.i32, shape=n_particles) # material id
+Jp = ti.var(dt=ti.f32, shape=n_particles) # plastic deformation
+grid_v = ti.Vector(2, dt=ti.f32, shape=(n_grid, n_grid)) # grid node momentum/velocity
+grid_m = ti.var(dt=ti.f32, shape=(n_grid, n_grid)) # grid node mass
+
+@ti.kernel
+def substep():
+  for i, j in grid_m:
+    grid_v[i, j] = [0, 0]
+    grid_m[i, j] = 0
+  for p in x: # Particle state update and scatter to grid (P2G)
+    base = (x[p] * inv_dx - 0.5).cast(int)
+    fx = x[p] * inv_dx - base.cast(float)
+    # Quadratic kernels  [http://mpm.graphics   Eqn. 123, with x=fx, fx-1,fx-2]
+    w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]
+    F[p] = (ti.Matrix.identity(ti.f32, 2) + dt * C[p]) @ F[p] # deformation gradient update
+    h = ti.exp(10 * (1.0 - Jp[p])) # Hardening coefficient: snow gets harder when compressed
+    if material[p] == 1: # jelly, make it softer
+      h = 0.3
+    mu, la = mu_0 * h, lambda_0 * h
+    if material[p] == 0: # liquid
+      mu = 0.0
+    U, sig, V = ti.svd(F[p])
+    J = 1.0
+    for d in ti.static(range(2)):
+      new_sig = sig[d, d]
+      if material[p] == 2:  # Snow
+        new_sig = min(max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3)  # Plasticity
+      Jp[p] *= sig[d, d] / new_sig
+      sig[d, d] = new_sig
+      J *= new_sig
+    if material[p] == 0:  # Reset deformation gradient to avoid numerical instability
+      F[p] = ti.Matrix.identity(ti.f32, 2) * ti.sqrt(J)
+    elif material[p] == 2:
+      F[p] = U @ sig @ V.T() # Reconstruct elastic deformation gradient after plasticity
+    stress = 2 * mu * (F[p] - U @ V.T()) @ F[p].T() + ti.Matrix.identity(ti.f32, 2) * la * J * (J - 1)
+    stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress
+    affine = stress + p_mass * C[p]
+    for i, j in ti.static(ti.ndrange(3, 3)): # Loop over 3x3 grid node neighborhood
+      offset = ti.Vector([i, j])
+      dpos = (offset.cast(float) - fx) * dx
+      weight = w[i][0] * w[j][1]
+      grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)
+      grid_m[base + offset] += weight * p_mass
+  for i, j in grid_m:
+    if grid_m[i, j] > 0: # No need for epsilon here
+      grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] # Momentum to velocity
+      grid_v[i, j][1] -= dt * 50 # gravity
+      if i < 3 and grid_v[i, j][0] < 0:          grid_v[i, j][0] = 0 # Boundary conditions
+      if i > n_grid - 3 and grid_v[i, j][0] > 0: grid_v[i, j][0] = 0
+      if j < 3 and grid_v[i, j][1] < 0:          grid_v[i, j][1] = 0
+      if j > n_grid - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0
+  for p in x: # grid to particle (G2P)
+    base = (x[p] * inv_dx - 0.5).cast(int)
+    fx = x[p] * inv_dx - base.cast(float)
+    w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1.0), 0.5 * ti.sqr(fx - 0.5)]
+    new_v = ti.Vector.zero(ti.f32, 2)
+    new_C = ti.Matrix.zero(ti.f32, 2, 2)
+    for i, j in ti.static(ti.ndrange(3, 3)): # loop over 3x3 grid node neighborhood
+      dpos = ti.Vector([i, j]).cast(float) - fx
+      g_v = grid_v[base + ti.Vector([i, j])]
+      weight = w[i][0] * w[j][1]
+      new_v += weight * g_v
+      new_C += 4 * inv_dx * weight * ti.outer_product(g_v, dpos)
+    v[p], C[p] = new_v, new_C
+    x[p] += dt * v[p] # advection
+
+group_size = n_particles // 3
+
+@ti.kernel
+def initialize():
+  for i in range(n_particles):
+    x[i] = [ti.random() * 0.2 + 0.3 + 0.10 * (i // group_size), ti.random() * 0.2 + 0.05 + 0.32 * (i // group_size)]
+    material[i] = i // group_size # 0: fluid 1: jelly 2: snow
+    v[i] = [0, 0]
+    F[i] = [[1, 0], [0, 1]]
+    Jp[i] = 1
+
+initialize()
+
+gui = ti.GUI("Taichi MLS-MPM-99", res=512, background_color=0x112F41)
+for frame in range(20000):
+  for s in range(int(2e-3 // dt)):
+    substep()
+  colors = np.array([0x068587, 0xED553B, 0xEEEEF0], dtype=np.uint32)
+  gui.circles(x.to_numpy(), radius=1.5, color=colors[material.to_numpy()])
+  gui.show() # Change to gui.show(f'{frame:06d}.png') to write images to disk

examples/opengl_nbody.py

@@ -0,0 +1,54 @@
+import taichi as ti
+
+ti.init(arch=ti.opengl)
+
+N = 8000
+pos = ti.var(dt=ti.f32, shape=(N, 2))
+vel = ti.var(dt=ti.f32, shape=(N, 2))
+bounce = 0.8
+
+@ti.kernel
+def initialize():
+  for i in range(N):
+    pos[i, 0] = ti.random() * 0.5 - 0.25 + 0.5
+    pos[i, 1] = ti.random() * 0.5 - 0.25 + 0.5
+    vel[i, 0] = ti.random() * 2 - 1
+    vel[i, 1] = ti.random() * 2 - 1
+
+@ti.kernel
+def advance(dt: ti.f32):
+  for i in range(N):
+    if pos[i, 0] < 0 and vel[i, 0] < 0 or pos[i, 0] > 1 and vel[i, 0] > 0:
+      vel[i, 0] = -bounce * vel[i, 0]
+    if pos[i, 1] < 0 and vel[i, 1] < 0 or pos[i, 1] > 1 and vel[i, 1] > 0:
+      vel[i, 1] = -bounce * vel[i, 1]
+
+    px = pos[i, 0] - 0.5
+    py = pos[i, 1] - 0.5
+    pas = -80.0 * ti.sqrt(px ** 2 + py ** 2)
+    a_x = px * pas
+    a_y = py * pas
+
+    for j in range(N):
+      if i != j:
+        dx = pos[i, 0] - pos[j, 0]
+        dy = pos[i, 1] - pos[j, 1]
+        d2 = dx ** 2 + dy ** 2
+        das = 233.0 * 0.001 ** 2 / (0.001 ** 2 + d2)
+        a_x = a_x + dx * das
+        a_y = a_y + dy * das
+
+    pos[i, 0] = pos[i, 0] + vel[i, 0] * dt
+    pos[i, 1] = pos[i, 1] + vel[i, 1] * dt
+
+    vel[i, 0] = vel[i, 0] + a_x * dt
+    vel[i, 1] = vel[i, 1] + a_y * dt
+
+gui = ti.GUI("n-body", res=(400, 400))
+
+initialize()
+while not gui.has_key_event() or gui.get_key_event().key == ti.GUI.MOTION:
+  _pos = pos.to_numpy()
+  gui.circles(_pos, radius=1, color=0x66ccff)
+  gui.show()
+  advance(0.005)

examples/opengl_random.py

@@ -0,0 +1,20 @@
+import taichi as ti
+
+ti.init(arch=ti.opengl)
+
+n = 640
+pixels = ti.var(dt=ti.f32, shape=(n, n))
+
+@ti.kernel
+def paint():
+  for i, j in pixels:
+    pixels[i, j] = ti.random()
+
+gui = ti.GUI("Noise", res=(n, n))
+
+for i in range(1000000):
+  paint()
+  gui.set_image(pixels)
+  gui.show()
+  if gui.has_key_event() and gui.get_key_event().key != ti.GUI.MOTION:
+    break

examples/opengl_range_for.py

@@ -0,0 +1,13 @@
+import taichi as ti
+
+ti.init(arch=ti.opengl)
+
+x = ti.var(ti.i32, shape=(4, 4))
+
+@ti.kernel
+def func():
+  for i, j in x:
+    x[i, j] = 200 + 10 * i + j
+
+func()
+print(x.to_numpy())

python/taichi/__init__.py

@@ -7,6 +7,7 @@
 from taichi.misc import *
 from taichi.misc.gui import GUI
 from taichi.misc.task import Task
+from taichi.misc.test import approx
 from taichi.misc import settings as settings
 from taichi.misc.gui import rgb_to_hex
 from taichi.misc.settings import *

python/taichi/lang/__init__.py

@@ -29,6 +29,7 @@
 x64 = core.x64
 cuda = core.cuda
 metal = core.metal
+opengl = core.opengl
 profiler_print = lambda: core.get_current_program().profiler_print()
 profiler_clear = lambda: core.get_current_program().profiler_clear()
 profiler_start = lambda n: core.get_current_program().profiler_start(n)
@@ -184,6 +185,8 @@ def supported_archs():
     archs.append(cuda)
   if ti.core.with_metal():
     archs.append(metal)
+  if ti.core.with_opengl():
+    archs.append(opengl)
   return archs
 
 class _ArchCheckers(object):

python/taichi/misc/test.py

@@ -0,0 +1,9 @@
+import taichi as ti
+
+def approx(*args, **kwargs):
+  from pytest import approx
+  if ti.cfg.arch == ti.opengl:
+    kwargs['rel'] = max(kwargs.get('rel', 1e6), 1e4)
+    return approx(*args, **kwargs)
+  else:
+    return approx(*args, **kwargs)