Taichi - Physics Engines

• Games 201 - ADVANCED PHYSICS ENGINES

  Lecturer: Dr. Yuanming Hu

• Taichi Graphics Course S1 (for completion and new feature updates)

  Lecturer: Dr. Tiantian Liu

Some supplement contents about Taichi

If overlapping contents: just add on the original GAMES201 notes. The new contents: add at the last of the notes individually.

• Exporting Results (Lecture 2)
• OOP and Metaprogramming (Lecture 3)
• Diff. Programming (Lecture 3)
• Visualization (Lecture 3)

Lecture 1 Introduction

Keyword: Discretization / Efficient solvers / Productivity / Performance / Hardware architecture / Data structures / Parallelization / Differentiability (DiffTaichi)

Taichi Programming Language

Initialization

Required every time using taichi

Use ti.init, for spec. cpu or gpu: use ti.cpu (default) or ti.gpu

ti.init(arch=ti.cuda) # spec run in any hardware

Data

Data Types

signed integers ti.i8 (~ 64, default i32) / unsigned integers ti.u8 (~64) / float-point numbers ti.f32(default) (~ 64)

Modify Default

The default can be changed via ti.init

• ti.init(default_fp=ti.f32)
• ti.init(default_ip=ti.i64)

Type Promotions

• i32 + f32 = f32
• i32 + i64 = i64

Switch to high percision automatically

Type Casts

• Implicit casts: Static types within the Taichi Scope

  xxxxxxxxxx
  def foo():  # Directly inside Python scope
      a = 1
      a = 2.7 # Python can re-def types automatically
      print(a)    # 2.7
  

  xxxxxxxxxx
  @ti.kernel  # Inside Taichi scope
  def foo():
      a = 1   # already def as a int type
      a = 2.7 # re-def failed (in Taichi)
      print(a)    # 2
  

• Explicit casts: variable = ti.casts(variable, type)

  xxxxxxxxxx
  @ti.kernel
  def foo():
      a = 1.8
      b = ti.cast(a, ti.i32)  # switch to int;    b = 1
      c = ti.cast(b, ti.f32)  # switch to floating;   c = 1.0
  

Tensor

Multi-dimensional arrays（高维数组）

• Self-defines: Use ti.types to create compound types including vector / matrix / struct

  ​x
  import taichi as ti
  ti.init(arch = ti.cpu)
  
  # Define your own types of data
  vec3f = ti.types.vector(3, ti.f32)  # 3-dim
  mat2f = ti.types.matrix(2, 2, ti.f32)   # 2x2
  ray = ti.types.struct(ro = vec3f, rd = vec3f, l = ti.f32)
  
  @ti.kernel
  def foo():
      a = vec3f(0.0)
      print(a)    # [0.0, 0.0, 0.0]
      d = vec3f(0.0, 1.0, 0.0)
      print(d)    # [0.0, 1.0, 0.0]
      B = mat2f([[1.5, 1.4], [1.3, 1.2]]) 
      print("B = ", B)    # B = [[1.5, 1.4], [1.3, 1.2]]
      r = ray(ro = a, rd = d, l = 1)
      print("r.ro = ", r.ro)  # r.ro = [0.0, 0.0, 0.0]
      print("r.rd = ", r.rd)  # r.rd = [0.0, 1.0, 0.0]
      
  foo()
  

• Pre-defines: An element of the tensors can be either a scalar (var), a vector (ti.Vector) or a matrix (ti.Matrix) (ti.Struct)

  Accessed via a[i, j, k]syntax (no pointers)

  xxxxxxxxxx
  import taichi as ti
  ti.init()
  a = ti.var(dt=ti.f32, shape=(42, 63)) # A tensor of 42x63 scalars
  b = ti.Vector(3, dt=ti.f32, shape=4) # A tensor of 4x 3D vectors (3 - elements in the vector, shape - shape of the tensor, composed by 4 3D vectors)
  C = ti.Matrix(2, 2, dt=ti.f32, shape=(3, 5)) # A tensor of 3x5 2x2 matrices
  loss = ti.var(dt=ti.f32, shape=()) # A 0-D tensor of a single scalar (1 element)
  
  a[3, 4] = 1
  print('a[3, 4] = ', a[3, 4]) # a[3, 4] = 1.000000
  b[2] = [6, 7, 8] 
  print('b[0] =', b[0][0], b[0][1], b[0][2]) # b[0] not yet supported (working)
  loss[None] = 3
  print(loss[None]) # 3
  

Field

ti.field: A global N-d array of elements

heat_field = ti.field(dtype=ti.f32, shape=(256, 256))

• global: a field can be read / written from both Taichi scope and Python scope

• N-d: Scalar (N = 0); Vector (N = 1); Matrix (N = 2); Tensor (N = 3, 4, 5, …)

• elements: scalar, vector, matrix, struct

• access elements in a field using [i, j, k, …] indexing

  xxxxxxxxxx
  import taichi as ti
  ti.init(arch=ti.cpu)
  
  pixels = ti.field(dtype=float, shape=(16, 8))
  pixels[1,2] = 42.0  # index the (1,2) pixel on the screen
  

  xxxxxxxxxx
  import taichi as ti
  ti.init(arch=ti.cpu)
  
  vf = ti.Vector.field(3, ti.f32, shape=4)    # 4x1 vectors, every vector is 3x1
  
  @ti.kernel
  def foo():
      v = ti.Vector([1, 2, 3])
      vf[0] = v
  

• Special Case: access a 0-D field using [None]

  xxxxxxxxxx
  zero_d_scalar = ti.field(ti.f32, shape=())
  zero_d_scalar[None] = 1.5   # Scalar in the scalar field
  
  zero_d_vector = ti.Vector.field(2, ti.f32, shape=())
  zero_d_vector[None] = ti.Vector([2.5, 2.6])
  

Other Examples:

• 3D gravitational field in a 256x256x128 room

  gravitational_field = ti.Vector.field(n=3, dtype=ti.f32, shape=(256, 256, 128))

• 2D strain-tensor field in a 64x64 grid

  strain_tensor_field = ti.Matrix.field(n = 2, m = 2, dtype = ti.f32, shape = (64, 64))

• a global scalar that want to access in a taichi kernel

  global_scalar = ti.field(dtype = ti.f32, shape=())

Kernels

Must be decorated with @ti.kernel (Compiled, statically-typed, lexically-scoped, parallel and differentiable - faster)

For loops at the outermost scope in a Taichi kernel is automatically parallelized (if inside - serial)

If a outermost scope is not wanted to parallelize for and an inside scope is, write the unwanted one in the python scope and call the other one as the outermost in the kernel.

Arguments

• At most 8 parameters

• Pass from Python scope to the Taichi scope

• Must be type-hinted

  xxxxxxxxxx
  @ti.kernel
  def my_kernel(x: ti.i32, y: ti.f32):    # explicit input variables with types
      print(x + y)
      
  my_kernel(2, 3.3)   # 5.3    
  

• Scalar only (if vector needs to input separately)

• Pass by value

  Actually copied from the var. in the Python scope and if the values of some var. x is modified in the Taichi kernel, it won’t change in the Python scope.

Return Value

• May or may not return

• Returns one single scalar value only

• Must be type-hinted

  xxxxxxxxxx
  @ti.kernel
  def my_kernel() -> ti.i32:  # returns i32
      return 233.666
  
  print(my_kernel())  # 233   (casted to int32)
  

Functions

Decorated with @ti.func, usually for high freq. used func.

Taichi's function can be called in taichi's kernel but can't be called by python (not global)

Function 可以被 kernel 调用但不能被 py 调用，kernel 也不能调用 kernel，但是 function可以调用 function

• Taichi functions can be nested (function in function)
• Taichi functions are force-inlined (cannot iterate)

Arguments

• Don’t need to be type-hinted

• Pass by value (for the force-inlined @ti.func)

  if want to pass outside, use return

Scalar Math

Similar to python. But Taichi also support chaining comparisons (a < b <= c != d)

Specifically, a / b outputs floating point numbers, a // b ouputs integer number (as in py3)

Matrices and Linear Algebra

ti.Matrix is for small matrices only. For large matrices, 2D tensor of scalars should be considered.

Note that element-wise product in taichi is * and mathematically defined matrix product is @ (like MATLAB)

Common Operations:（返回矩阵 A 经过变换的结果而非变换 A 本身）

xxxxxxxxxx
A.transpose() # 转置
A.inverse() # 逆矩阵
A.trace() # 迹
A.determinant(type)
A.normalized() # 向量除以自己模长
A.norm() # 返回向量长度
A.cast(type) # 转换为其他数据类型
ti.sin(A)/cos(A) # element wise（给 A 中所有元素都做 sin 或者 cos 运算）

Parallel for-loops

(Automatically parallalized)

Range-for loops

Same as python for loops

Struct-for loops

Iterates over (sparse) tensor elements. Only lives at the outermost scope looping over a ti.field

for i,j in x: can be used

Example:

xxxxxxxxxx
import taichi as ti

ti.init(arch=ti.gpu)  # initiallize every time

n = 320
pixels = ti.var(dt=ti.f32, shape=(n*2, n))  # every element in this tensor is 32 bit float point number，shape 640x320

@ti.kernel
def paint(t: ti.f32):
    for i, j in pixels:  # parallized over all pixel
        pixels[i, j] = i * 0.001 + j * 0.002 + t 
        # This struct-for loops iterate over all tensor coordinates. i.e. (0,0), (0,1), (0,2), ..., (0,319), (1,0), ..., (639,319)
        # Loop for all elements in the tensor
        # For this dense tensor, every element is active. But for sparse tensor, struct-for loop will only work for active elements
paint(0.3)

• break is NOT supported in the parallel for-loops

Atomic Operations

In Taichi, augmented assignments (e.g. x[i] += 1) are automatically atomic.

(+=: the value on the RHS is directly summed on the current value of the var. and the referncing var.(array))

(Atomic operation: an operation that will always be executed without any other process being able to read or change state that is read or changed during the operation. It is effectively executed as a single step, and is an important quality in a number of algorithms that deal with multiple independent processes, both in synchronization and algorithms that update shared data without requiring synchronization.)

When modifying global var. in parallel, make sure to use atomic operations.

xxxxxxxxxx
@ti.kernel
def sum():  # to sum up all elements in x
    for i in x:
        # Approach 1: OK
        total[None] += x[i]  # total: 0-D tensor => total[None]
        
        # Approach 2: OK
        ti.atomic_add(total[None], x[i])  # will return value before the atomic addition
        
        # Approach 3: Wrong Result (Not atomic)
        total[None] = total[None] + x[i]    # other thread may have summed

Taichi-scope vs. Python-scope

Taichi-scope: Everything decorated with ti.kernel and ti.func

Code in Taichi-scope will be compiled by the Taichi compiler and run on parallel devices (Attention that in parallel devices the order of print may not be guaranteed)

• Static data type in the Taichi scope

  The type won’t change even if actually defines some other values / fields (error)

• Static lexical scope in the Taichi scope

  xxxxxxxxxx
  @ti.kernel
  def err_out_of_scope(x:float):
      if x < 0:   # abs value of x
          y = -x
      else:
          y = x
      print(y)    # y generated in the 'if' and not pass outside, error occurs
  

• Compiled JIT (just in time) (cannot see Python scope var. at run time)

  xxxxxxxxxx
  a = 42
  
  @ti.kernel
  def print_a():
      print('a = ', a)
  
  print_a()   # 'a = 42'
  a = 53
  print('a = ', a)    # 'a = 53'
  print_a()   # still 'a = 42'
  

  Another demo:

  xxxxxxxxxx
  d = 1
  
  @ti.kernel
  def foo():
      print('d in Taichi scope = ', d)
      
  d += 1  # d = 2
  foo()   # d in Taichi scope = 2 (but after this call, ti kernel will regard d as a constant)
  d += 1  # d = 3
  foo()   # d in Taichi scope = 2 (d not changed in Ti-scope but changed in Py-scope)
  

  If want real global: use ti.field

  xxxxxxxxxx
  a = ti.field(ti.i32, shape=())
  
  @ti.kernel
  def print_a():
      print('a=', a[None])
  
  a[None] = 42
  print_a() # "a= 42"
  a[None] = 53
  print_a() # "a= 53"
  

Python-scope: Code outside the Taichi-scope

Code in Python-scope is simple py code and will be executed by the py interpreter

Phases

1. Initialization: ti.init()

2. Tensor allocation: ti.var, ti.Vector, ti.Matrix

   Only define tensors in this allocation phase and never define in the computation phase

3. Computation (launch kernel, access tensors in Py-scope)

4. Optional: Restart (clear memory, destory var. and kernels): ti.reset()

Practical Example (fractal)

xxxxxxxxxx
import taichi as ti

ti.init(arch=ti.gpu)

n = 320
pixels = ti.var(dt=ti.f32, shape=(n * 2, n))

@ti.func
def complex_sqr(z):
    return ti.Vector([z[0]**2 - z[1]**2, z[1] * z[0] * 2])  # calculate the square of a complex number
    # In this example, use a list and put in [] to give values for ti.Vector 

@ti.kernel
def paint(t: ti.f32):   # time t - float point 32 bit
    for i, j in pixels:  # parallized over all pixels
        c = ti.Vector([-0.8, ti.cos(t) * 0.2])
        z = ti.Vector([i / n - 1, j / n - 0.5]) * 2  # Julia set formula
        iterations = 0   # iterate for all pixels
        while z.norm() < 20 and iterations < 50:
            z = complex_sqr(z) + c  # user-defined in @ti.func
            iterations += 1
        pixels[i, j] = 1 - iterations * 0.02
        
gui = ti.GUI("Julia Set", res = (n*2, n))   # ti's gui

for i in range(1000000):
    paint(i * 0.03)
    gui.set_image(pixels)
    gui.show

Debug Mode

Lecture 2 Lagrangian View

Two Views (Def)

Lagrangian View

Sensors that move passively with the simulated material（节点跟随介质移动）

粒子会被额外记录位置，会随（被模拟的）介质不断移动

Euler View

Still sensors that never moves（穿过的介质的速度）

网格中每个点的位置不会被特别记录

Mass-Spring Systems

弹簧 - 质点模型 (Cloth / Elastic objects / ...)

Hooke's Law

(f - force, k - Spring stifness, L_i,j - Spring rest length between particle i and j, ^ - normalization, ${\widehat{{\mathbf{x}}_i-\mathbf{x}}}_j$$\widehat{\mathbf{x}_{i}-\mathbf{x}}_{j}$​ - direction vector of paticle i to j)

Newton's Second Law of Motion

Time Integration

Common Types of Integrators

• Forward Euler (explict)

  ${\mathbf{v}}_{i+1}={\mathbf{v}}_t+\mathrm{\Delta }t\frac{{\mathbf{f}}_t}{m}$$\mathbf{v}_{i+1} = \mathbf{v}_t + \Delta t \frac{\mathbf{f}_t}{m}$, ${\mathbf{x}}_{t+1}={\mathbf{x}}_t+\mathrm{\Delta }t{\mathbf{v}}_t$$\mathbf{x}_{t+1}=\mathbf{x}_{t}+\Delta t \mathbf{v}_{t}$

• Semi-implicit Euler (aka. Symplectic Euler, Explicit)

  ${\mathbf{v}}_{i+1}={\mathbf{v}}_t+\mathrm{\Delta }t\frac{{\mathbf{f}}_t}{m}$$\mathbf{v}_{i+1} = \mathbf{v}_t + \Delta t \frac{\mathbf{f}_t}{m}$, ${\mathbf{x}}_{t+1}={\mathbf{x}}_t+\mathrm{\Delta }t{\mathbf{v}}_{t+1}$$\mathbf{x}_{t+1}=\mathbf{x}_{t}+\Delta t \mathbf{v}_{t+1}$

  （准确性上提升，以 $t+1$$t+1$​​​ 步的 $\mathbf{v}$$\mathbf{v}$​​​ 替代第 $t$$t$​​ 步的使用）

• Backward Euler (often with Newton's Method, Implicit)

  ${\mathbf{x}}_{t+1}={\mathbf{x}}_t+\mathrm{\Delta }t{\mathbf{v}}_{t+1}$$\mathbf{x}_{t+1}=\mathbf{x}_{t}+\Delta t \mathbf{v}_{t+1}$​, ${\mathbf{v}}_{i+1}={\mathbf{v}}_t+\mathrm{\Delta }t{\mathbf{\text{M}}}^{-1}\mathbf{f}({\mathbf{x}}_{t+1})$$\mathbf{v}_{i+1} = \mathbf{v}_t + \Delta t \textbf{M}^{-1} \mathbf{f}(\mathbf{x}_{t+1})$​ （相互依赖）

  Full implement see later

Comparison

• Explicit: (forward Euler / symplectic Euler / RK (2, 3, 4) / ...)

  Future depends only on past, easy to implement, easy to explode, bad for stiff materials

  $\mathrm{\Delta }t\le c\sqrt{\frac{m}{k}}(c\sim 1)$$\Delta t \leq c\sqrt{\frac{m}{k}}\; (c\sim 1)$​（当粒子质量变大的时候允许的 $\mathrm{\Delta }t$$\Delta t$ 就会变大，$k$$k$​（硬度）越大的时候允许的也会变小 - 不适合硬度大的）

  ~ 数值稳定性，不会衰减而是随时间指数增长

• Implicit (back Euler / middle-point)

  Future depends on both future and past, harder to implement, need to solve a syustem of (linear) equation, hard to optimize, time steps become more expensive but time steps are larger, extra numerical damping and locking (overdamping) (but generally, uncoonditionally stable)

  （显式 - 容易实现，数值上较为不稳定，受 dt 步长影响较大；隐式 - 难实现，允许较长步长）

Implementing a Mass-Spring System with Symplectic Euler

Steps

• Compute new velocity using ${\mathbf{v}}_{i+1}={\mathbf{v}}_t+\mathrm{\Delta }t\frac{{\mathbf{f}}_t}{m}$$\mathbf{v}_{i+1} = \mathbf{v}_t + \Delta t \frac{\mathbf{f}_t}{m}$
• Collision with ground
• Compute new position using ${\mathbf{x}}_{t+1}={\mathbf{x}}_t+\mathrm{\Delta }t{\mathbf{v}}_{t+1}$$\mathbf{x}_{t+1}=\mathbf{x}_{t}+\Delta t \mathbf{v}_{t+1}$

Demo

xxxxxxxxxx
import taichi as ti

ti.init(debug=True)

max_num_particles = 256

dt = 1e-3

num_particles = ti.var(ti.i32, shape=())
spring_stiffness = ti.var(ti.f32, shape=())
paused = ti.var(ti.i32, shape=())
damping = ti.var(ti.f32, shape=())

particle_mass = 1
bottom_y = 0.05

x = ti.Vector(2, dt=ti.f32, shape=max_num_particles)
v = ti.Vector(2, dt=ti.f32, shape=max_num_particles)

A = ti.Matrix(2, 2, dt=ti.f32, shape=(max_num_particles, max_num_particles))
b = ti.Vector(2, dt=ti.f32, shape=max_num_particles)

# rest_length[i, j] = 0 means i and j are not connected
rest_length = ti.var(ti.f32, shape=(max_num_particles, max_num_particles))

connection_radius = 0.15

gravity = [0, -9.8]

@ti.kernel
def substep():  # 每一个新的时间步，把每一帧分为若干步（在模拟中 time_step 需要取相应较小的数值）
    
    # Compute force and new velocity
    n = num_particles[None]
    for i in range(n):  # 枚举全部 i
        v[i] *= ti.exp(-dt * damping[None]) # damping
        total_force = ti.Vector(gravity) * particle_mass  # 总受力 G = mg
        for j in range(n):  # 枚举其余所有粒子 是否有关系
            if rest_length[i, j] != 0:  # 两粒子之间有弹簧？
                x_ij = x[i] - x[j]
                total_force += -spring_stiffness[None] * (x_ij.norm() - rest_length[i, j]) * x_ij.normalized()  
                # 胡克定律公式 -k * ((xi - xj).norm - rest_length) * (xi - xj).norm
        v[i] += dt * total_force / particle_mass    # sympletic euler: 用力更新一次速度
        
    # Collide with ground （计算完力和速度之后立刻与地面进行一次碰撞）
    for i in range(n):
        if x[i].y < bottom_y:   # 一旦发现有陷入地下的趋势就把速度的 y component 设置为 0
            x[i].y = bottom_y
            v[i].y = 0

    # Compute new position
    for i in range(num_particles[None]):
        x[i] += v[i] * dt   # 新的位置的更新：把速度产生的位移累加到位置上去

        
@ti.kernel
def new_particle(pos_x: ti.f32, pos_y: ti.f32): # Taichi doesn't support using Matrices as kernel arguments yet
    new_particle_id = num_particles[None]
    x[new_particle_id] = [pos_x, pos_y]
    v[new_particle_id] = [0, 0]
    num_particles[None] += 1
    
    # Connect with existing particles
    for i in range(new_particle_id):
        dist = (x[new_particle_id] - x[i]).norm()
        if dist < connection_radius:
            rest_length[i, new_particle_id] = 0.1
            rest_length[new_particle_id, i] = 0.1
    
    
gui = ti.GUI('Mass Spring System', res=(512, 512), background_color=0xdddddd)

spring_stiffness[None] = 10000
damping[None] = 20

new_particle(0.3, 0.3)
new_particle(0.3, 0.4)
new_particle(0.4, 0.4)

while True:
    for e in gui.get_events(ti.GUI.PRESS):
        if e.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]:
            exit()
        elif e.key == gui.SPACE:
            paused[None] = not paused[None]
        elif e.key == ti.GUI.LMB:
            new_particle(e.pos[0], e.pos[1])
        elif e.key == 'c':
            num_particles[None] = 0
            rest_length.fill(0)
        elif e.key == 's':
            if gui.is_pressed('Shift'):
                spring_stiffness[None] /= 1.1
            else:
                spring_stiffness[None] *= 1.1
        elif e.key == 'd':
            if gui.is_pressed('Shift'):
                damping[None] /= 1.1
            else:
                damping[None] *= 1.1
                
    if not paused[None]:
        for step in range(10):
            substep()
    
    X = x.to_numpy()
    gui.circles(X[:num_particles[None]], color=0xffaa77, radius=5)
    
    gui.line(begin=(0.0, bottom_y), end=(1.0, bottom_y), color=0x0, radius=1)
    
    for i in range(num_particles[None]):
        for j in range(i + 1, num_particles[None]):
            if rest_length[i, j] != 0:
                gui.line(begin=X[i], end=X[j], radius=2, color=0x445566)
    gui.text(content=f'C: clear all; Space: pause', pos=(0, 0.95), color=0x0)
    gui.text(content=f'S: Spring stiffness {spring_stiffness[None]:.1f}', pos=(0, 0.9), color=0x0)
    gui.text(content=f'D: damping {damping[None]:.2f}', pos=(0, 0.85), color=0x0)
    gui.show()

Backward Euler Implicit Implement

Steps

• ${\mathbf{x}}_{t+1}={\mathbf{x}}_t+\mathrm{\Delta }t{\mathbf{v}}_{t+1}$$\mathbf{x}_{t+1}=\mathbf{x}_{t}+\Delta t \mathbf{v}_{t+1}$​, ${\mathbf{v}}_{i+1}={\mathbf{v}}_t+\mathrm{\Delta }t{\mathbf{\text{M}}}^{-1}\mathbf{f}({\mathbf{x}}_{t+1})$$\mathbf{v}_{i+1} = \mathbf{v}_t + \Delta t \textbf{M}^{-1} \mathbf{f}(\mathbf{x}_{t+1})$​ （相互依赖）
• Eliminate ${\mathbf{x}}_{t+1}$$\mathbf{x}_{t+1}$: ${\mathbf{v}}_{i+1}={\mathbf{v}}_t+\mathrm{\Delta }t{\mathbf{\text{M}}}^{-1}\mathbf{f}({\mathbf{x}}_t+\mathrm{\Delta }t{\mathbf{v}}_{t+1})$$\mathbf{v}_{i+1} = \mathbf{v}_t + \Delta t \textbf{M}^{-1} \mathbf{f}(\mathbf{x}_{t}+\Delta t\mathbf{v}_{t+1})$​ （f 通常为非线性，用下述泰勒展开一次线性化）
• Linearize (Newton's): ${\mathbf{v}}_{t+1}={\mathbf{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}[\mathbf{f}\left({\mathbf{x}}_t\right)+\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\left({\mathbf{x}}_t\right)\mathrm{\Delta }t{\mathbf{v}}_{t+1}]$$\mathbf{v}_{t+1}=\mathbf{v}_{t}+\Delta t \mathbf{M}^{-1}\left[\mathbf{f}\left(\mathbf{x}_{t}\right)+\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\left(\mathbf{x}_{t}\right) \Delta t \mathbf{v}_{t+1}\right]$​
• Clean up: $[\mathbf{I}-\mathrm{\Delta }{t}^2{\mathbf{M}}^{-1}\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\left({\mathbf{x}}_t\right)]{\mathbf{v}}_{t+1}={\mathbf{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}\mathbf{f}\left({\mathbf{x}}_t\right)$$\left[\mathbf{I}-\Delta t^{2} \mathbf{M}^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{x}}\left(\mathbf{x}_{t}\right)\right] \mathbf{v}_{t+1}=\mathbf{v}_{t}+\Delta t \mathbf{M}^{-1} \mathbf{f}\left(\mathbf{x}_{t}\right)$​ (Linear System)
• To solve the equation: Jacobi / Gauss-Seidel iteration OR conjugate gradients

Solving the system:

Demo

xxxxxxxxxx
import taichi as ti
import random

ti.init()

n = 20  

A = ti.var(dt=ti.f32, shape=(n, n))  # 20 x 20 Matrix
x = ti.var(dt=ti.f32, shape=n)
new_x = ti.var(dt=ti.f32, shape=n)
b = ti.var(dt=ti.f32, shape=n)

@ti.kernel  # iteration kernel
def iterate():
    for i in range(n):
        r = b[i]
        for j in range(n):
            if i != j:
                r -= A[i, j] * x[j]
                
        new_x[i] = r / A[i, i]  # 每次都更新 x 使得 x 能满足矩阵一行或一个的线性方程组（局部更新）- 对性质好的矩阵逐渐收敛
        
    for i in range(n):
        x[i] = new_x[i]

@ti.kernel  # Compute residual b - A * x
def residual() -> ti.f32:   # residual 一开始会非常大 经过若干次迭代会降到非常低
    res = 0.0
    
    for i in range(n):
        r = b[i] * 1.0
        for j in range(n):
            r -= A[i, j] * x[j]
        res += r * r
        
    return res

for i in range(n):  # 初始化矩阵
    for j in range(n):
        A[i, j] = random.random() - 0.5

    A[i, i] += n * 0.1
    
    b[i] = random.random() * 100

for i in range(100):    # 执行迭代
    iterate()   
    print(f'iter {i}, residual={residual():0.10f}')
    

for i in range(n):
    lhs = 0.0
    for j in range(n):
        lhs += A[i, j] * x[j]
    assert abs(lhs - b[i]) < 1e-4

Unifying Explicit and Implicit Integrators

• $\beta =0$$\beta = 0$, forward / semi-implicit Euler (Explicit)
• $\beta =1/2$$\beta = 1/ 2$, middle point (Implicit)
• $\beta =1$$\beta = 1$, backward Euler (Implicit)

Solve Faster

For millions of mass points and springs

• Sparse Matrices
• Conjugate Gradients
• Preconditioning
• Use Position-based Dynamics

Lagrangian Fluid Simulation (SPH)

Smoothed Particle Hydrodynamics

Use particles carrying samples of physical quntities, and a kernel $W$$W$​, to approximate continuous fields ($A$$A$ can be almost any spatially varying phisical attributes: density / pressure / ...)

在 x 这点的物理场的值（即 $A(\mathbf{x})$$A(\mathbf{x})$​​​）等于核函数（如密度 / 压力等）加权平均了周围的粒子（本质是光滑了）的求和。核函数中间高四周低，使得更靠近该点影响更大，远离该点的影响变小（超过 support 半径 $h$$h$​ 则为 0）

不需要 mesh，适合模拟自由表面 free-surface flows 流体（如水和空气接触的表面，反之，不适合烟雾（空间需要被填满）等），可以使用 “每个粒子就是以小包水” 理解

Equation of States (EOS)

aka. Weakly Compressible SPH (WCSPH)

Momentum Equation ($\rho$$\rho$ - density (${\rho }_0$$\rho_0$ - ideal), $B$$B$ - bulk modulus, $\gamma$$\gamma$ - constant (usually ~ 7) )

(Actually the Navier-Stoke's without considering viscosity)

${\rho }_i$$\rho_i$​ 通过质量加权, $D$$D$ - Lagrange derivative

Gradient in SPH

Not really accurate but at least symmetric and momentum conserving (to add viscosity etc. Laplacian should be introduced)

SPH Simulation Cycle

• For each particle $i$$i$, compute ${\rho }_i=\sum _j{m}_{j}W(||{\mathbf{x}}_i-{\mathbf{x}}_j|{|}_2,h)$$\rho_i = \sum_j m_j W(||\mathbf{x}_i - \mathbf{x}_j||_2,h)$

• For each particle $i$$i$, compute $\mathrm{\nabla }{p}_i$$\nabla p_i$ using the gradient operator

• Symplectic Euler steps (similar as the mass-spring model)

  $${\mathbf{v}}_{t+1}={\mathbf{v}}_t+\mathrm{\Delta }t\frac{D\mathbf{v}}{Dt},{\mathbf{x}}_{t+1}={\mathbf{x}}_t+\mathrm{\Delta }t\frac{D\mathbf{v}}{Dt}$$

Variant of SPH

• Predictive-Corrective Impcompressible SPH (PCI-SPH) - 无散，更接近不可压缩
• Position-based Fluids (PBF) (Demo:ti eample pbf2d) - Position-based dynamics + SPH
• Divergence-free SPH (DFSPH) - (velocity)

Paper: SPH Fluids in Computer Graphics, Smooted Particle Hydrodynamics Techniques for Physics Based Simulation of Fluids and Solids (2019)

Courant-Friedrichs-Lewy (CFL) Conditions (Explicit)

One upper bound of time step size: (def. using velocity other than stiffness)

($C$$C$ - Courant number, CFL number, $u$$u$ - maximum velocity, $\mathrm{\Delta }t$$\Delta t$ - time step, $\mathrm{\Delta }x$$\Delta x$ - length interval (e.g. particle radius and grid size))

Application: estimating allowed time step in (explicit) time integration.

SPH: ~ 0.4; MPM (Material Point Method): 0.3 ~ 1; FLIP fluid (smoke): 1 ~ 5+

Accelerating SPH: Neighborhood Search

Build spatial data structure such as voxel grids $O(n^2)\to O(n)$$O(n^2) \rightarrow O(n)$

Neighborhood search with hashing

精确找到距离不超过 $h$$h$​​​​ 的粒子，每个 voxel 维护一个 list（所有位置在这个 voxel 里的粒子），每个需要查找的粒子先确定相应 voxel 再查找周围的 voxel（枚举所有包含的粒子，通常是常数）- SPH 的瓶颈

Ref: Compact hashing

Other Paricle-based Simulation Methods

• Discrete Element Method (DEM) - 刚体小球，如沙子模拟
• Moving Particle Semi-implicit (MPS) - 增强 fluid incompressibility
• Power Particles - Incompressible fluid solver based on power diagrams
• A peridynamics perspective on spring-mass fracture

Exporting Results

Exporting Videos

• ti.GUI.show: save screenshots / ti.imwrite(img, filename)
• ti video: creates video.mp4 (-f 24 / -f 60)
• ti gif -i input.mp4: Convert mp4 to gif

Lecture 3 Lagrangian View (2)

弹性、有限元基础、Taichi 高级特性

Elasticity

Deformation

Deformation Map $\varphi$$\phi$​​ (to describe): A (vector to vector) function that relates rest material position and deformed material position

${\mathbf{x}}_{\mathrm{deformed}}=\varphi ({\mathbf{x}}_{\mathrm{rest}})$$\mathbf{x}_{\mathrm{deformed}} = \phi (\mathbf{x}_{\mathrm{rest}})$

Deformation gradient $\mathbf{F}$$\mathbf{F}$: $\mathbf{F}=\frac{\partial {\mathbf{x}}_{\mathrm{deformed}}}{\partial {\mathbf{x}}_{\mathrm{rest}}}$$\mathbf{F} = \frac{\partial \mathbf{x}_{\mathrm{deformed}}}{\partial \mathbf{x}_{\mathrm{rest}}}$ (Translational invariant, same deformation gradients for ${\varphi }_1$$\phi_1$ and ${\varphi }_1+c$$\phi_1 + c$)

Deform / rest volume ratio: $J=\det (\mathbf{F})$$J = \mathrm{det} (\mathbf{F})$ （F[None] = [[2, 0], [0, 1]]（横向拉伸了））

Elasticity

Hyperelastic material

Hyperelastic material: Stress-strain relationship is defined by a strain energy density function $\psi =\psi (\mathbf{F})$$\psi = \psi(\mathbf{F})$

$\psi$$\psi$​ is a potential function that penalizes deformation; $\mathbf{F}$$\mathbf{F}$ represents the strain (means deformation gradient)

Stress Tensor

• The First Piola-Kirchhoff stress tensor (PK1): $\mathbf{P}\mathbf{(}\mathbf{F})=\frac{\partial \psi (\mathbf{F})}{\partial \mathbf{F}}$$\mathbf{P(F)} = \frac{\partial \psi(\mathbf{F})}{\partial \mathbf{F}}$​ (easy to compute but in rest space, asymmetric)
• Kirchhoff stress: $\tau$$\tau$
• Cauchy stress tensor: $\sigma$$\sigma$ (True stress, symmetric, conservation of angular momentum)

Relationship: $\tau =j\sigma ={\mathbf{PF}}^T$$\tau = j\sigma = \mathbf{PF}^T$; $\mathbf{P}=J\sigma {\mathbf{F}}^{-T}$$\mathbf{P} = J\sigma \mathbf{F}^{-T}$; Traction $t={\sigma }^T\mathbf{n}$$t = \sigma^{T} \mathbf{n}$ （考虑的是截面的法向量，转置为 ^-T）​

Elastic Moduli (Isotropic Materials)

• Young’s Modulus: $E=\frac{\sigma }{ϵ}$$E = \frac{\sigma }{\epsilon}$
• Bulk Modulus: $K=-V\frac{dP}{dV}$$K = -V \frac{dP}{dV}$ （常用与可压缩液体）
• Poisson’s Ratio: $\nu \in [0.0,0.5]$$\nu \in [0.0, 0.5]$ （泊松比越大，越接近于不可压缩；Auxetics 材料可以拥有负泊松比，施加拉力材料反而变粗）

Lame parameters

• Lame’s first parameter: $\mu$$\mu$
• Lame’s second parameter (shear modulus / $G$$G$): $\lambda$$\lambda$​

Conversions:（通常指定 Young’s Modulus & Poisson’s Ratio，或 $E$$E$ 和 $\nu$$\nu$）$K=\frac{E}{3(1-2\nu )}$$K = \frac{E}{3(1-2\nu)}$ ; $\lambda =\frac{E\nu }{(1+\nu )(1-2\nu )}$$\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}$ ; $\mu =\frac{E}{2(1+\nu )}$$\mu = \frac{E}{2(1+\nu)}$

Hyperelastic Material Models

Popular in graphics:

• Linear elasticity (small deformation only, not consistence in rotation) - linear equation systems

• Neo-Hookean: (Commonly used in isotropic materials)

  • $\psi (\mathbf{F})=\frac{\mu }{2}\sum _i[{\left({\mathbf{F}}^T\mathbf{F}\right)}_{ii}-1]-\mu \log \left(J\right)+\frac{\lambda }{2}{\log }^2(J)$$\psi(\mathbf{F})=\frac{\mu}{2} \sum_{i}\left[\left(\mathbf{F}^{T} \mathbf{F}\right)_{i i}-1\right]-\mu \log (J)+\frac{\lambda}{2} \log ^{2}(J)$
  • $\mathbf{P}(\mathbf{F})=\frac{\partial \psi }{\partial \mathbf{F}}=\mu (\mathbf{F}-{\mathbf{F}}^T)+\lambda \log \left(J\right){\mathbf{F}}^{-T}$$\mathbf{P}(\mathbf{F})=\frac{\partial \psi}{\partial \mathbf{F}}=\mu\left(\mathbf{F}-\mathbf{F}^{T}\right)+\lambda \log (J) \mathbf{F}^{-T}$

• (Fixed) Corotated:

  • $\psi (\mathbf{F})=\mu \sum _i{({\sigma }_i-1)}^2+\frac{\lambda }{2}(J-1{)}^2$$\psi(\mathbf{F})=\mu \sum_{i}\left(\sigma_{i}-1\right)^{2}+\frac{\lambda}{2}(J-1)^{2}$ (${\sigma }_i$$\sigma_i$ are singlular values of $\mathbf{F}$$\mathbf{F}$)
  • $\mathbf{P}(\mathbf{F})=\frac{\partial \psi }{\partial \mathbf{F}}=2\mu (\mathbf{F}-\mathbf{R})+\lambda (J-1)J{\mathbf{F}}^{-T}$$\mathbf{P}(\mathbf{F})=\frac{\partial \psi}{\partial \mathbf{F}}=2 \mu(\mathbf{F}-\mathbf{R})+\lambda(J-1) J \mathbf{F}^{-T}$

Lagrangian Finite Elements on Linear Tetrahedral Meshes

The Finite Element Method (FEM)

Discretization sheme to build discrete equations using weak formulations of continuous PDEs

Linear Tetrahedral (Triangular) FEM

Assumption: The deformation map $\varphi$$\phi$​ is affine and thereby the deformation gradient $\mathbf{F}$$\mathbf{F}$​ is constant

${\mathbf{x}}_{\mathrm{deformed}}=\mathbf{F}({\mathbf{x}}_{\mathrm{rest}})+\mathbf{b}$$\mathbf{x}_{\mathrm{deformed}} = \mathbf{F}(\mathbf{x}_{\mathrm{rest}}) + \mathbf{b}$ （$\mathbf{b}$$\mathbf{b}$ - offset，平移产生的形不产生弹性，不影响 gradient）

For every element $e$$e$​, the elastic potential energy $U(e)={\int }_e\psi (\mathbf{F}\mathbf{(}\mathbf{x}))d\mathbf{x}=V_e\psi ({\mathbf{F}}_e)$$U(e) = \int_e\psi (\mathbf{F(x)}) d\mathbf{x} = V_e \psi (\mathbf{F} _e)$​​ （对每个元素这么设定，计算其弹性势能）

Computing ${\mathbf{F}}_e$$\mathbf{F}_e$​ (for 2D triangles)

• Vertices: ${\mathbf{a}}_{\mathrm{rest}}$$\mathbf{a}_{\mathrm{rest}}$​; Deformed positions: ${\mathbf{a}}_{\mathrm{deformed}}$$\mathbf{a}_\mathrm{deformed}$​ (same as $\mathbf{b}$$\mathbf{b}$​ and $\mathbf{c}$$\mathbf{c}$​)

• ${\mathbf{a}}_{\mathrm{deformed}}={\mathbf{Fa}}_{\mathrm{deformed}}+\mathbf{P}$$\mathbf{a}_\mathrm{deformed} = \mathbf{Fa}_\mathrm{deformed} + \mathbf{P}$ (same as $\mathbf{b}$$\mathbf{b}$ and $\mathbf{c}$$\mathbf{c}$​)

• Eliminate $\mathbf{P}$$\mathbf{P}$: $({\mathrm{a}}_{\text{deformed}}-{\mathbf{c}}_{\text{deformed}})=\mathrm{F}({\mathrm{a}}_{\text{rest}}-{\mathrm{c}}_{\text{rest}})$$\left(\mathrm{a}_{\text {deformed}}-\mathbf{c}_{\text {deformed}}\right)=\mathrm{F}\left(\mathrm{a}_{\text {rest}}-\mathrm{c}_{\text {rest}}\right)$ $({\mathrm{b}}_{\text{deformed}}-{\mathbf{c}}_{\text{deformed}})=\mathrm{F}({\mathrm{b}}_{\text{rest}}-{\mathrm{c}}_{\text{rest}})$$\left(\mathrm{b}_{\text {deformed}}-\mathbf{c}_{\text {deformed}}\right)=\mathrm{F}\left(\mathrm{b}_{\text {rest}}-\mathrm{c}_{\text {rest}}\right)$

• ${\mathbf{F}}_{(2\times 2)}$$\mathbf{F}_{(2\times 2)}$ has 4 linear constrains (equations)

  $\mathbf{B}=[{\mathbf{a}}_{\mathrm{rest}}-{\mathbf{c}}_{\mathrm{rest}}|{\mathbf{b}}_{\mathrm{rest}}-{\mathbf{c}}_{\mathrm{rest}}{]}^{-1}$$\mathbf{B} = [\mathbf{a}_\mathrm{rest} - \mathbf{c}_\mathrm{rest} | \mathbf{b}_\mathrm{rest} - \mathbf{c}_\mathrm{rest}]^{-1}$ ; $\mathbf{D}=[{\mathbf{a}}_{\mathrm{deformed}}-{\mathbf{c}}_{\mathrm{deformed}}|{\mathbf{b}}_{\mathrm{deformed}}-{\mathbf{c}}_{\mathrm{deformed}}]$$\mathbf{D} = [\mathbf{a}_\mathrm{deformed} - \mathbf{c}_\mathrm{deformed} | \mathbf{b}_\mathrm{deformed} - \mathbf{c}_\mathrm{deformed}]$ ; $\mathbf{F}\mathbf{=}\mathbf{DB}$$\mathbf{F = DB}$​

  $\mathbf{B}$$\mathbf{B}$ is constant throughout the process, should be pre-computed

Explicit Linear Triangular (FEM) Simulation

Semi-implicit Euler time integration scheme:

${\mathbf{v}}_{t+1,i}={\mathbf{v}}_{t,i}+\mathrm{\Delta }t\frac{{\mathbf{f}}_{t,i}}{m_i}{\mathbf{x}}_{t+1,i}={\mathbf{x}}_{t,i}+\mathrm{\Delta }t{\mathbf{v}}_{t+1,i}$$\mathbf{v}_{t+1, i} =\mathbf{v}_{t, i}+\Delta t \frac{\mathbf{f}_{t, i}}{m_{i}} \quad \mathbf{x}_{t+1, i} =\mathbf{x}_{t, i}+\Delta t \mathbf{v}_{t+1, i}$​

${\mathbf{x}}_{t,i}$$\mathbf{x}_{t,i}$​ and ${\mathbf{v}}_{t,i}$$\mathbf{v}_{t,i}$​​ are stored on the vertices of finite elements (triangles / tetrahedrons), $V_e$$V_e$ - constant volumes

${\mathbf{f}}_{t,i}=-\frac{\partial U}{\partial {\mathbf{x}}_i}=-\sum _e\frac{\partial U(e)}{\partial {\mathbf{x}}_i}=-\sum _e{V}_e\frac{\partial \psi \left({\mathbf{F}}_e\right)}{\partial {\mathbf{F}}_e}\frac{\partial {\mathbf{F}}_e}{{\partial }_{{\mathbf{x}}_i}}=-\sum _e{V}_e\mathbf{P}\left({\mathbf{F}}_e\right)\frac{\partial {\mathbf{F}}_e}{\partial {\mathbf{x}}_i}$$\mathbf{f}_{t, i}=-\frac{\partial U}{\partial \mathbf{x}_{i}}=-\sum_{e} \frac{\partial U(e)}{\partial \mathbf{x}_{i}}=-\sum_{e} V_{e} \frac{\partial \psi\left(\mathbf{F}_{e}\right)}{\partial \mathbf{F}_{e}} \frac{\partial \mathbf{F}_{e}}{\partial_{ \mathbf{x}_{i}}}=-\sum_{e} V_{e} \mathbf{P}\left(\mathbf{F}_{e}\right) \frac{\partial \mathbf{F}_{e}}{\partial \mathbf{x}_{i}}$

Taichi’s AutoDiff system can use to compute $\mathbf{p}({\mathbf{F}}_e)$$\mathbf{p}(\mathbf{F}_e)$​

xxxxxxxxxx
for s in range(30):
    with ti.Tape(total_energy):  # Automatically diff. and write into a 'tape' and use in semi-implicit Euler

Implicit Linear Triangular FEM Simulation

Backward Euler Time Integration: $[\mathbf{I}-\beta \mathrm{\Delta }{t}^2{\mathbf{M}}^{-1}\frac{\partial \mathbf{f}}{\partial \mathbf{x}}\left({\mathbf{x}}_t\right)]{\mathbf{v}}_{t+1}={\mathbf{v}}_t+\mathrm{\Delta }t{\mathbf{M}}^{-1}\mathbf{f}\left({\mathbf{x}}_t\right)$$\left[\mathbf{I}-\beta \Delta t^{2} \mathbf{M}^{-1} \frac{\partial \mathbf{f}}{\partial \mathbf{x}}\left(\mathbf{x}_{t}\right)\right] \mathbf{v}_{t+1}=\mathbf{v}_{t}+\Delta t \mathbf{M}^{-1} \mathbf{f}\left(\mathbf{x}_{t}\right)$

($\mathbf{M}$$\mathbf{M}$ - diagonal matrix to record the mass)

Compute force differentials: $\frac{\partial \mathbf{f}}{\partial \mathbf{x}}=-\frac{{\partial }^2\psi }{\partial {\mathbf{x}}^2}$$\frac{\partial \mathbf{f}}{\partial \mathbf{x}} = -\frac{\partial ^2 \psi}{\partial \mathbf{x}^2}$ (second order DE are not available in taichi)

Higher Level of Taichi

Objective Data-Oriented Programming

OOP -> Extensibility / Maintainability

An “object” contains its own data (py-var. / ti.field) and method (py-func. / @ti.func / @ti.kernel)

Python OOP in a Nutshell

xxxxxxxxxx
class Wheel:    # def a class of 'wheel'
    def __init__(self, radius, width, rolling_fric):    # data, for a wheel: radius / width / fric
        self.radius = radius    # convert all data as `self.` (inside member of 'self')
        self.width = width
        self.rolling_fric = rolling_fric
        
    def Roll(self): # Method of the 'wheel'
        ...
         
w1 = Wheel(5, 1, 0.1)   # Instantiated Objects (different wheels)
w2 = Wheel(5, 1, 0.1)
w3 = Wheel(6, 1.2, 0.15)
w4 = Wheel(6, 1.2, 0.15)

If want to add some features, this method can inherit all past features and add on them (easy to maintain and extent)

Using Classes in Taichi

Hybrid scheme: Objective Data-Oriented Programming (ODOP)

-> More data-oriented -> usually use more vectors / matrices in the class

Important Decorators:

• @ti.data_oriented to decorate class
• @ti.kernel to decorate class members functions that are Taichi kernels
• @ti.func to decorate class members functions that are Taichi functions

Caution: if the variable is passed from Python scope, the self.xxx will still regard as a constant

Features:

• Encapsulation: Different classes can be stored in various .py scripts and can be called using from [filename] import [classname] in the main script.
• Inheritance: A @data_oriented class can inherit from another @data_oriented class. Both the data and methods are inherited from the base class.
• Polymorphism: Define methods in the child class that have the same name as the methods in the parent class. Proper methods will be called according to the instantiated objects.

Demo: ti example odop_solar: $\mathbf{a}=GM\mathbf{r}/||\mathbf{r}|{|}_2^3$$\mathbf{a} = GM\mathbf{r}/||\mathbf{r}||^3_2$

xxxxxxxxxx
import taichi as ti

@ti.data_oriented
class SolarSystem:
    def __init__(self, n, dt):  # n - planet number; dt - time step size
        self.n = n
        self.dt = dt
        self.x = ti.Vector(2, dt=ti.f32, shape=n)
        self.v = ti.Vector(2, dt=ti.f32, shape=n)
        self.center = ti.Vector(2, dt=ti.f32, shape=())
        
    @staticmethod   # @ti.func 还可以被额外的 @staticmethod（静态函数，同 py）修饰
    @ti.func    # 生成一个随机数
    def random_around(center, radius):  # random number in [center - radius, center + radius]
        return center + radius * (ti.random() - 0.5)  * 2
    
    @ti.kernel
    def initialize(self):   # initialize all the tensors
        for i in range (self.n):
            offset = ti.Vector([0.0, self.random_around(0.3, 0.15)])
            self.x[i] = self.center[None] + offset
            self.v[i] = [-offset[1], offset[0]]
            self.v[i] *= 1.5 / offset.norm()
            
    @ti.func    # still in class
    def gravity(self, pos):
        offset = -(pos - self.center[None])
        return offset / offset.norm()**3
    
    @ti.kernel  # sympletic Euler
    def integrate(self):
        for i in range(self.n):
            self.v[i] += self.dt * self.gravity(self.x[i])
            self.x[i] += self.dt * self.v[i]
            
solar = SolarSystem(9, 0.0005)  # 9 for n (planet number); 0.0005 for dt
solar.center[None] = [0.5, 0.5]
solar.initialize()

gui = ti.GUI("Solar System", background_color = 0x25A6D9)

while True:
    if gui.get_event():
        if gui.event.key == gui.SPACE and gui.event.type == gui.PRESS:
            solar.initialize()
    for i in range(10):
        solar.integrate()
    gui.circle([0.5, 0.5], radius = 20, color = 0x8C274C)
    gui.circles(solar.x.to_numpy(), radius = 5, color = 0xFFFFFF)
    gui.show()