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 all 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
74 changes: 44 additions & 30 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 @@ -62,13 +64,13 @@ 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)
# 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 @@ -77,11 +79,11 @@ def sample_minmax(vf, p):
u, v = 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)
iu, iv = ti.floor(s), ti.floor(t)
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 @@ -113,8 +115,8 @@ 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()):
def advect_semilag(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
Expand All @@ -123,21 +125,31 @@ 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
# star means the temp value after a back tracing (forward advection)
# two star means the temp value after a forward tracing (reverse advection)
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)

new_qf[i, j] = q_star + 0.5 * (qf[i, j] - new_qf[i, j])

min_val, max_val = sample_minmax(qf, p_mid)
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 +187,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 +244,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 +278,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 @@ -347,9 +361,9 @@ def reset():

gui.set_image(dyes_pair.cur)
# To visualize velocity field:
#gui.set_image(velocities_pair.cur.to_numpy() * 0.01 + 0.5)
# gui.set_image(velocities_pair.cur.to_numpy() * 0.01 + 0.5)
# To visualize velocity divergence:
#divergence(velocities_pair.cur); gui.set_image(velocity_divs.to_numpy() * 0.1 + 0.5)
# divergence(velocities_pair.cur); gui.set_image(velocity_divs.to_numpy() * 0.1 + 0.5)
# To visualize velocity vorticity:
#vorticity(velocities_pair.cur); gui.set_image(velocity_curls.to_numpy() * 0.03 + 0.5)
# vorticity(velocities_pair.cur); gui.set_image(velocity_curls.to_numpy() * 0.03 + 0.5)
gui.show()