EWaveAlgorithm.pdf

(148 KB) Pobierz
eWave: Using an Exponential Solver on the
iWave Problem
Jerry Tessendorf
March 16, 2014
1
iWave Dynamics
The iWave equations of motion for the surface displacement
h
and velocity
potential
φ
are:
∂h(x, t)
∂t
∂φ(x, t)
∂t
=
−
2
φ(x, t)
(1)
(2)
(3)
=
−gh(x,
t)
To see how the exponential solver approach helps, it would be good to set
up a different way to express these equations in terms of a couplet
W
(x,
t)
=
The couplet has the equation of motion
∂W
(x,
t)
=
M
W
(x,
t)
∂t
where
M
is the matrix
M
=
0
−g
√
−
0
2
h(x, t)
φ(x, t)
(4)
(5)
(6)
2
eWave: Exponential Solution
W
(x,
t)
= exp
{Mt}
W
(x, 0)
(7)
Equation 5 has the exact exponential solution
1
This can be built more explicitly, because if you work it out, you see the following
identities:
M
2
M
3
M
2n
M
2n+1
=
=
=
=
−
1 0
0 1
g
−
2
(8)
(9)
−M
g
−
2
1 0
−g
0 1
M
−g
−
2
−
n
2
n
(10)
(11)
So, if we expand the exponential into a Taylor series:
exp
{Mt}
=
(M)
n
t
n
n!
n=0
∞
(12)
and separate the powers into even and odd sets
(M)
2n
t
2n
exp
{Mt}
=
+
(2n)!
n=0
Using the identities above,
exp
{Mt}
=
1
0
0
1
√
(−g
−
2
t
2
)
n
+
Mt
(2n)!
n=0
∞
∞
(M)
2n+1
t
2n+1
(2n + 1)!
n=0
∞
(13)
√
(−g
−
2
t
2
)
n
(14)
(2n + 1)!
n=0
∞
The first infinite series is the series for the cosine, and the second is the series
√
for the sine. So if we define
ω
≡
g
−
2
ˆ
exp
{Mt}
=
=
1
0
0
1
cos (ˆ
t)
+
ω
√
−
ω
ˆ
2
M
sin (ˆ
t)
ω
ω
ˆ
(15)
(16)
cos (ˆ
t)
ω
−g
ω
ω
sin (ˆ
t)
ˆ
sin (ˆ
t)
ω
cos (ˆ
t)
ω
The nice thing about this is that, because
g
and
2
are not time dependent
quantities, this solution is exact for any value of time
t.
This is the exact solution
for any time period:
√
−
2
sin (ˆ ∆t)
φ(x, t)
(17)
ω
h(x, t
+ ∆t) = cos (ˆ ∆t)
h(x, t)
+
ω
ω
ˆ
g
φ(x, t
+ ∆t) = cos (ˆ ∆t)
φ(x, t)
−
ω
sin (ˆ ∆t)
h(x, t)
ω
(18)
ω
ˆ
In real space, the tricky part of this is the complex way that
2
is in the
solution. Using convolution on this would be similar to the old iWave approach,
but still better because the dynamics has been integrated already.
2
3
FFT Form
In Fourier space, the appearance of
2
simplifies considerable, because now the
√
operator
ω
becomes the dispersion relation
ω
→
ω(k)
=
gk,
where
k
is the
ˆ
ˆ
absolute magnitude of the Fourier vector. For the initial Fourier amplitudes
˜
˜
h(k)
and
φ(k)
˜
h(k, t
+ ∆t)
˜
φ(k, t
+ ∆t)
k
sin (ω(k)∆t)
ω(k)
g
˜
sin (ω(k)∆t)
cos (ω(k)∆t)
φ(k, t)
−
ω(k)
˜
cos (ω(k)∆t)
h(k, t)
+
˜
φ(k, t)
(19)
˜
h(k, t)
(20)
=
=
Equations 19 and 20 are what is in the FFTDynamics code.
4
Convolution Form
If we want to evaluate the dynamics as a convolution, the FFT form can be
used to construction the convolutions. In convolution form, there are three
convolution kernels that apply as
h(x, t)
φ(x, t)
=
=
C(x)
⊗
h(x)
+
S(x)
⊗
φ(x)
C(x)
⊗
φ(x)
−
T
(x)
⊗
h(x)
(21)
(22)
where
⊗
denotes spatial convolution and the convolution kernels are:
C(x)
S(x)
T
(x)
=
=
=
d
2
k
i
k
·
x
e
cos (ω(k)∆t)
(2π)
2
d
2
k
i
k
·
x
k
e
sin (ω(k)∆t)
(2π)
2
ω(k)
d
2
k
i
k
·
x
g
e
sin (ω(k)∆t)
(2π)
2
ω(k)
(23)
(24)
(25)
Explicitly in terms of integration, these convolutions are:
h(x, t)
φ(x, t)
=
=
d
2
y C(x
−
y)
h(y)
+
S(x
−
y)
φ(y)
d
2
y C(x
−
y)
φ(y)
−
T
(x
−
y)
h(y)
(26)
(27)
The practical implementation of convolution is as a moving window filter on
the 2D gridded data for
h
and
φ.
Imagining these quantities have values
h
ij
,
φ
ij
on the rectangular grid, the time updates are
N/2
h
ij
(t + ∆t)
=
k,l=−N/2
C
kl
h
i+k
j+l
(t)
+
S
kl
φ
i+k
j+l
(t)
(28)
3
N/2
φ
ij
(t + ∆t)
=
k,l=−N/2
C
kl
φ
i+k
j+l
(t)
−
T
kl
h
i+k
j+l
(t)
(29)
Where
N
is the size of the square moving window, and
C
kl
,
S
kl
, and
T
kl
are the
values our the kernels at discrete grid intervals. Strategies for obtaining these
discrete kernels are discussed below.
Theoretically, if we take the limit that
N
is the size of the full simulation
grid and apply periodic boundary conditions, the moving window filter produces
identical results to the FFT approach, although the FFT approach is dramati-
cally faster at it. But the moving window filter opens up many other strategies
that can be applied when simple periodic boundary conditions are no longer fea-
sible, and assuming that a small-enough window size
N
produces results that
look good. In iWave, the smallest window size recommended was 13×13, but
here a size that looks good will probably be larger, perhaps 30×30 or larger
depending on the circumstances and desired visual qualities. In general, the
larger the value of
N
, the better the quality of the propagation motion.
5
Boundary Conditions
The convolution approach breaks away from the constraint that the fields
h
and
φ
be periodic. That freedom also imposes a requirement that boundary
conditions be specified when applying the convolution. Here are a few different
strategies. Note that all of these strategies can be applied in a single simulation
by choosing different ones on each boundary.
5.1
Periodic
Of course, periodic boundary conditions can still be imposed, simply by wrap-
ping the indexing of
i
+
k
and
j
+
l
in the convolution sums.
5.2
Fixed Ghost Values
When the values of
i
+
k
or
j
+
l
extend beyond the bound of the grid, ghost
values for
h
and
φ
can be prescribed.
5.3
Tiled Simulation Grids
If you are running independent eWave simulations on multiple grids that butt
up agains each other, the simultions can be coupled via the boundary conditions.
When the values of
i
+
k
or
j
+
l
extend beyond the bound of one grid, values for
h
and
φ
can be retrieved from the adjoining grid. This will naturally propagate
waves and momentum between simulation grids without requiring extending
each grid with gridpoints for ghost values (although ghost values is one approach
to implementing the approach here).
4
6
Convolution Kernel Construction
The approach below is to find expressions for the convolution kernels
C, S,
and
T
. There are two methods we will look at. The first is a highly theoretical
way that produces practically poor performance. The second is a brute-force
method that is fast, flexible, and works.
6.1
The Impractical Theoretical Way
For the deep water case, using the symmetries of the integrands, we can define
the dimensionless variable
ξ
as
ξ
2
=
gt
2
|x|
(30)
Using this variable, these convolution kernels can be written as
C(x)
=
S(x)
=
1
1
C(ξ)
|x|
2
S(ξ)
T
(ξ)
(31)
(32)
(33)
√
|x|
5/2
g
√
g
|x|
3/2
T
(x) =
where
∞
C(ξ)
=
0
∞
√
du u J
0
(u) cos(
u ξ)
(34)
(35)
(36)
S(ξ)
=
0
∞
√
du u
3/2
J
0
(u) sin(
u ξ)
√
du u
1/2
J
0
(u) sin(
u ξ)
T
(ξ) =
0
While this approach to generating the convolution kernels may prove in-
teresting in the future, trying to apply these formula has not yet produced a
practically useful approach.
6.2
The Practical Brute-Force Way
The simplest way to quickly and accurately construct
C, S,
and
T
as
N
×
N
moving window filters is to use FFTs to directly evaluate the integrals in
equations 23 – 25. The procedure is as follows:
1. Create three grids, one each for
C, S,
and
T
. The dimensions of these grids
should
not
be
N
×
N
. They should be the dimensions of the simulation
grid for
h
and
φ.
This is to ensure that spatial scales properly contribute
to the filter values.
5

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika:

Zgłoś jeśli naruszono regulamin