Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Example] Fix stable_fluid.py boundary broken: int -> ti.floor, and make mpm3d.py 3x faster on CUDA #1784

Merged
merged 3 commits into from
Aug 27, 2020
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 13 additions & 8 deletions examples/mpm3d.py
Original file line number Diff line number Diff line change
@@ -1,15 +1,18 @@
export_file = '/tmp/mpm3d.ply'
export_file = '' # use '/tmp/mpm3d.ply' for exporting result to disk

import taichi as ti
import numpy as np

ti.init(arch=ti.opengl)
ti.init(arch=ti.gpu)

#dim, n_grid, steps, dt = 2, 128, 20, 2e-4
#dim, n_grid, steps, dt = 2, 256, 32, 1e-4
dim, n_grid, steps, dt = 3, 32, 25, 4e-4
#dim, n_grid, steps, dt = 3, 64, 25, 2e-4
#dim, n_grid, steps, dt = 3, 128, 25, 8e-5

dim = 3
n_grid = 32
n_particles = n_grid**dim // 2**(dim - 1)
dx = 1 / n_grid
dt = 4e-4

p_rho = 1
p_vol = (dx * 0.5)**2
Expand All @@ -34,6 +37,7 @@ def substep():
for I in ti.grouped(grid_m):
grid_v[I] = grid_v[I] * 0
grid_m[I] = 0
ti.block_dim(n_grid)
for p in x:
Xp = x[p] / dx
base = int(Xp - 0.5)
Expand All @@ -55,6 +59,7 @@ def substep():
cond = I < bound and grid_v[I] < 0 or I > n_grid - bound and grid_v[
I] > 0
grid_v[I] = 0 if cond else grid_v[I]
ti.block_dim(n_grid)
for p in x:
Xp = x[p] / dx
base = int(Xp - 0.5)
Expand All @@ -79,7 +84,7 @@ def substep():
@ti.kernel
def init():
for i in range(n_particles):
x[i] = ti.Vector([ti.random() for i in range(dim)]) * 0.4 + 0.2
x[i] = ti.Vector([ti.random() for i in range(dim)]) * 0.4 + 0.15
J[i] = 1


Expand All @@ -101,12 +106,12 @@ def T(a):
init()
gui = ti.GUI('MPM3D', background_color=0x112F41)
while gui.running and not gui.get_event(gui.ESCAPE):
for s in range(15):
for s in range(steps):
substep()
pos = x.to_numpy()
if export_file:
writer = ti.PLYWriter(num_vertices=n_particles)
writer.add_vertex_pos(pos[:, 0], pos[:, 1], pos[:, 2])
writer.export_frame(gui.frame, export_file)
gui.circles(T(pos), radius=2, color=0x66ccff)
gui.circles(T(pos), radius=1.5, color=0x66ccff)
gui.show()
3 changes: 2 additions & 1 deletion examples/stable_fluid.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,8 @@ def bilerp(vf, p):
u, v = p
s, t = u - 0.5, v - 0.5
# floor
iu, iv = int(s), int(t)
iu, iv = ti.floor(s), ti.floor(t)
Copy link
Collaborator Author

@archibate archibate Aug 26, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

>>> import math
>>> math.floor(-1.1)
-2
>>> int(-1.1)
-1

So int is quite different from floor when things turn negative. The same to Taichi-scope.
This causes the -x and -y side having funny effects when non-zero velocity is encountered.
Hi, dear @Eydcao:
Did you assumed int to be exactly same with ti.floor, return the largest integer below the input argument during the whole debugging process?
It warns that we should never do over-assumptions when debugging. If you find all possible cause doesn't match, then you might overlooked some tiny places. Hope we all have fun and learnt a lot from this vivid lesson :)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting!

but in the old code, we used int() instead of math.floor(), so in python these two are the same thing?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Both in Python and Taichi, they are different thing :)

Copy link
Contributor

@Eydcao Eydcao Aug 26, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great! This one is good to merge.

I have got the idea why int() is irrigating here as it leads to the wrong coefficients in lerp. So when I fixed these two lines in #1764.

The strange boundary artifact is gone, great!

The new bfecc works perfec with jacobi_single now, but still has problems with jacobi_dual and mgpcg. Nonetheless, those two are the problems concerning the PPE solving, and I assume the problems with boundary artifact and bfecc are solved. So will organize a bit and make #1764 ready for review.

iu, iv = int(s + 1) - 1, int(t + 1) - 1
archibate marked this conversation as resolved.
Show resolved Hide resolved
# fract
fu, fv = s - iu, t - iv
a = sample(vf, iu + 0.5, iv + 0.5)
Expand Down