Graphics Reference
In-Depth Information
We can apply that reasoning on our grid to get the semi-Lagrangian
method introduced to graphics by Stam [Stam 99]. We want to figure out
the new value of
q
at a grid point, and to do that in a Lagrangian way we
need to figure out the old value of
q
that the particle that ends up at the
grid point possesses. The particle is moving through the velocity field
u
,
and we know where it ends up—so to find where it started we simply run
backwards through the velocity field from the grid point. We can grab the
old value of
q
from this start point, which will be the new value of
q
at the
gridpoint!Butwait,yousay,whatifthatstartpointwasn'tonthegrid?
In that case we simply interpolate the old value of
q
from the old values on
the grid, and we're done.
Let's go through that again, slowly and with formulas. We'll say that
the location in space of the grid point we're looking at is
x
G
.Wewant
to find the new value of
q
at that point, which we'll call
q
n
+
G
.Weknow
from our understanding of advection that if a hypothetical particle with
old value
q
P
ends up at
x
G
, when it moves through the velocity field for
thetimestepΔ
t
,then
q
n
+1
G
=
q
P
. So the question is, how do we figure out
q
P
?
The first step is figuring out where this imaginary particle would have
started from, a position we'll call
x
P
. The particle moves according to the
simple ordinary differential equation
dx
dt
=
u
(
x
)
and ends up at
x
G
after time Δ
t
. If we now run time backwards, we can
go in reverse from
x
G
to the start point of the particle—i.e., finding where
a particle would end up under the reverse velocity field
u
“starting” from
x
G
. The simplest possible way to estimate this is to use one step of forward
Euler:
−
x
P
=
x
G
−
Δ
tu
(
x
G
)
,
where we use the velocity evaluated at the grid point to take a Δ
t
-step
backwards through the flow field. It turns out forward Euler is sometimes
adequate, but significantly better results can be obtained using a slightly
more sophisticated technique such as a higher-order Runge-Kutta method.
See Appendix A to review time integration methods. In particular, at least
a second-order Runge-Kutta method is recommended as a default, such as
1
2
Δ
tu
(
x
G
)
,
x
mid
=
x
G
−
x
P
=
x
G
−
Δ
tu
(
x
mid
)
.
