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.

Sign up for GitHub

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

Jump to bottom

[bug] Fix stable fluid BFECC advection #1764

Merged
merged 8 commits into from
Sep 26, 2020
Merged

[bug] Fix stable fluid BFECC advection #1764

Changes from 3 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
3e0da34
apply known fixings, fix the bfecc process, the correction symbol, an…
Eydcao Aug 24, 2020
34de871
[skip ci] add standard test; fix some mistakes
Eydcao Aug 25, 2020
808c52c
[skip ci] remove the abundont '+0.5'; adjust the boundary to non-boun…
Eydcao Aug 25, 2020
55c02f7
Apply suggestions from code review
Eydcao Aug 26, 2020
0fd590f
[skip ci] apply suggestions from Yubin concerning floor and int
Eydcao Aug 26, 2020
ede5c88
[skip ci] clean
Eydcao Aug 26, 2020
d076a0c
[skip ci]Apply suggestions from code review
Eydcao Aug 26, 2020
ae0c782
[skip ci] enforce code format
taichi-gardener Aug 26, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
84 changes: 55 additions & 29 deletions examples/stable_fluid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,14 @@
ti.init(arch=ti.gpu)

_velocities = ti.Vector.field(2, float, shape=(res, res))
_intermedia_velocities = ti.Vector.field(2, float, shape=(res, res))
_new_velocities = ti.Vector.field(2, float, shape=(res, res))
velocity_divs = ti.field(float, shape=(res, res))
velocity_curls = ti.field(float, shape=(res, res))
_pressures = ti.field(float, shape=(res, res))
_new_pressures = ti.field(float, shape=(res, res))
_dye_buffer = ti.Vector.field(3, float, shape=(res, res))
_intermedia_dye_buffer = ti.Vector.field(3, float, shape=(res, res))
_new_dye_buffer = ti.Vector.field(3, float, shape=(res, res))


Expand Down Expand Up @@ -65,10 +67,10 @@ def bilerp(vf, p):
iu, iv = int(s), int(t)
# fract
fu, fv = s - iu, t - iv
a = sample(vf, iu + 0.5, iv + 0.5)
b = sample(vf, iu + 1.5, iv + 0.5)
c = sample(vf, iu + 0.5, iv + 1.5)
d = sample(vf, iu + 1.5, iv + 1.5)
a = sample(vf, iu, iv)
b = sample(vf, iu+1, iv)
c = sample(vf, iu, iv+1)
d = sample(vf, iu+1, iv+1)
return lerp(lerp(a, b, fu), lerp(c, d, fu), fv)


Expand All @@ -78,10 +80,10 @@ def sample_minmax(vf, p):
s, t = u - 0.5, v - 0.5
# floor
iu, iv = int(s), int(t)
a = sample(vf, iu + 0.5, iv + 0.5)
b = sample(vf, iu + 1.5, iv + 0.5)
c = sample(vf, iu + 0.5, iv + 1.5)
d = sample(vf, iu + 1.5, iv + 1.5)
a = sample(vf, iu, iv)
b = sample(vf, iu+1, iv)
c = sample(vf, iu, iv+1)
d = sample(vf, iu+1, iv+1)
return min(a, b, c, d), max(a, b, c, d)


Expand Down Expand Up @@ -114,7 +116,7 @@ def backtrace_rk3(vf: ti.template(), p, dt: ti.template()):

@ti.kernel
def advect_semilag(vf: ti.template(), qf: ti.template(),
new_qf: ti.template()):
new_qf: ti.template(), intermedia_qf: ti.template()):
ti.cache_read_only(qf, vf)
for i, j in vf:
p = ti.Vector([i, j]) + 0.5
Expand All @@ -123,21 +125,29 @@ def advect_semilag(vf: ti.template(), qf: ti.template(),


@ti.kernel
def advect_bfecc(vf: ti.template(), qf: ti.template(), new_qf: ti.template()):
def advect_bfecc(vf: ti.template(), qf: ti.template(), new_qf: ti.template(), intermedia_qf: ti.template()):
ti.cache_read_only(qf, vf)
for i, j in vf:
p = ti.Vector([i, j]) + 0.5
p_mid = backtrace(vf, p, dt)
q_mid = bilerp(qf, p_mid)
p_fin = backtrace(vf, p_mid, -dt)
q_fin = bilerp(qf, p_fin)
new_qf[i, j] = q_mid + 0.5 * (q_fin - qf[i, j])
p = backtrace(vf, p, dt)
intermedia_qf[i, j] = bilerp(qf, p)

ti.cache_read_only(intermedia_qf, qf, vf)
for i, j in vf:
Copy link
Collaborator

Choose a reason for hiding this comment

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

Is the key problem about separating bfecc to two offload?

Copy link
Contributor Author

@Eydcao Eydcao Aug 26, 2020 •
edited

Choose a reason for hiding this comment

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

Yes. To be more precise, we need a forward advection (backtrack + interpolate) to get an inter-media quantity field first, then a backward advection to (backtrace in -dt direction + interpolate) to get the field at the original time stamp (assuming there is no error). Finally, we got the error and fix the inter-media field.

Since for the second one (backward advection) we need the whole inter-media field ready, it's necessary to run it in a separate another loop and store it already.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

But to be noticed, jacobbi_dual caused serious artifacts with the new bfecc.

I will post detailed comparisons between combination of different methods when I get a minute.

p = ti.Vector([i, j]) + 0.5
p_two_star = backtrace(vf, p, -dt)
Eydcao marked this conversation as resolved.
Show resolved Hide resolved
p_star = backtrace(vf, p, dt)
q_star = intermedia_qf[i, j]
new_qf[i, j] = bilerp(intermedia_qf, p_two_star)

min_val, max_val = sample_minmax(qf, p_mid)
# new
Eydcao marked this conversation as resolved.
Show resolved Hide resolved
new_qf[i, j] = q_star + 0.5 * (qf[i, j]-new_qf[i, j])

min_val, max_val = sample_minmax(qf, p_star)
cond = min_val < new_qf[i, j] < max_val
for k in ti.static(range(cond.n)):
if not cond[k]:
new_qf[i, j][k] = q_mid[k]
new_qf[i, j][k] = q_star[k]


advect = advect_bfecc
Expand Down Expand Up @@ -175,13 +185,13 @@ def divergence(vf: ti.template()):
vt = sample(vf, i, j + 1).y
vc = sample(vf, i, j)
if i == 0:
vl = -vc.x
vl = 0
if i == res - 1:
vr = -vc.x
vr = 0
if j == 0:
vb = -vc.y
vb = 0
if j == res - 1:
vt = -vc.y
vt = 0
velocity_divs[i, j] = (vr - vl + vt - vb) * 0.5


Expand Down Expand Up @@ -232,7 +242,7 @@ def pressure_jacobi_dual(pf: ti.template(), new_pf: ti.template()):
(plt + prt + prb + plb) * 2 + pcc * 4) * 0.0625


pressure_jacobi = pressure_jacobi_dual
pressure_jacobi = pressure_jacobi_single

if pressure_jacobi == pressure_jacobi_dual:
p_jacobi_iters //= 2
Expand Down Expand Up @@ -266,8 +276,10 @@ def enhance_vorticity(vf: ti.template(), cf: ti.template()):


def step(mouse_data):
advect(velocities_pair.cur, velocities_pair.cur, velocities_pair.nxt)
advect(velocities_pair.cur, dyes_pair.cur, dyes_pair.nxt)
advect(velocities_pair.cur, velocities_pair.cur,
velocities_pair.nxt, _intermedia_velocities)
advect(velocities_pair.cur, dyes_pair.cur,
dyes_pair.nxt, _intermedia_dye_buffer)
velocities_pair.swap()
dyes_pair.swap()

Expand Down Expand Up @@ -296,17 +308,23 @@ def __init__(self):
self.prev_mouse = None
self.prev_color = None

def __call__(self, gui):
def __call__(self, gui, frame):
# [0:2]: normalized delta direction
# [2:4]: current mouse xy
# [4:7]: color
mouse_data = np.zeros(8, dtype=np.float32)
if gui.is_pressed(ti.GUI.LMB):
mxy = np.array(gui.get_cursor_pos(), dtype=np.float32) * res
if gui.is_pressed(ti.GUI.LMB) or frame < 20:
if frame < 20:
Eydcao marked this conversation as resolved.
Show resolved Hide resolved
# print(frame)
mxy = np.array([0 + frame*0.01, 0 + frame*0.01],
dtype=np.float32) * res
else:
mxy = np.array(gui.get_cursor_pos(), dtype=np.float32) * res
if self.prev_mouse is None:
self.prev_mouse = mxy
# Set lower bound to 0.3 to prevent too dark colors
self.prev_color = (np.random.rand(3) * 0.7) + 0.3
# self.prev_color = (np.random.rand(3) * 0.7) + 0.3
self.prev_color = np.array([1, 1, 1], dtype=np.float32)
else:
mdir = mxy - self.prev_mouse
mdir = mdir / (np.linalg.norm(mdir) + 1e-5)
Expand All @@ -328,6 +346,7 @@ def reset():

gui = ti.GUI('Stable Fluid', (res, res))
md_gen = MouseDataGen()
frame = 0
Eydcao marked this conversation as resolved.
Show resolved Hide resolved
while gui.running:
if gui.get_event(ti.GUI.PRESS):
e = gui.event
Expand All @@ -341,8 +360,15 @@ def reset():
elif e.key == 'd':
debug = not debug

frame += 1
Eydcao marked this conversation as resolved.
Show resolved Hide resolved

if frame > 4500:
Eydcao marked this conversation as resolved.
Show resolved Hide resolved
break

if not paused:
mouse_data = md_gen(gui)
mouse_data = md_gen(gui, frame)
Eydcao marked this conversation as resolved.
Show resolved Hide resolved
# if frame < 40:
# print(mouse_data)
step(mouse_data)

gui.set_image(dyes_pair.cur)
Expand Down