From http://www.jwz.org/xscreensaver/xscreensaver-5.30.tar.gz
[xscreensaver] / hacks / euler2d.tex
1 \documentclass[12pt]{article}
2
3 %\usepackage{fullpage}
4 \usepackage{amsmath,amssymb}
5
6 \begin{document}
7
8 \title{Two Dimensional Euler Simulation}
9
10 \author{
11 S.J. Montgomery-Smith\\
12 Department of Mathematics\\
13 University of Missouri\\
14 Columbia, MO 65211, U.S.A.\\
15 stephen@math.missouri.edu\\
16 http://www.math.missouri.edu/\~{}stephen}
17
18 \date{September 10, 2000}
19
20 \maketitle
21
22 This document describes a program I wrote to simulate the 
23 two dimensional Euler Equation --- a program that is part
24 of the {\tt xlock} screensaver as the {\tt euler2d}
25 mode.  A similar explanation may also be found in the
26 book by Chorin \cite{C}.
27
28 \section{The Euler Equation}
29
30 The Euler Equation describes the motion of an incompressible
31 fluid that has no viscosity.  If the fluid is contained
32 in a domain $\Omega$ with boundary $\partial \Omega$, then
33 the equation is in the vector field $u$ (the velocity)
34 and the
35 scalar field $p$ (the pressure):
36 \begin{eqnarray*}
37 \frac{\partial}{\partial t} u &=& -u \cdot \nabla u + \nabla p \\
38 \nabla \cdot u &=& 0 \\
39 u \cdot n &=& 0 \quad \text{on $\partial \Omega$}
40 \end{eqnarray*}
41 where $n$ is the unit normal to $\partial \Omega$.
42
43 \section{Vorticity}
44
45 It turns out that it can be easier write these equations
46 in terms of the vorticity.  In two dimensions the vorticity
47 is the scalar $w = \partial u_2/\partial x - \partial u_1/\partial y$.
48 The equation for vorticity becomes
49 \[ \frac{\partial}{\partial t} w = -u \cdot \nabla w .\]
50 A solution to this equation can be written as follows.  The velocity
51 $u$ causes a flow, that is, a function $\varphi(t,x)$ that tells where
52 the particle initially at $x$ ends up at time $t$, that is
53 \[
54 \frac\partial{\partial t} \varphi(t,x)
55 = u(t,\varphi(t,x)) .\]
56 Then the equation
57 for $w$ tells us that the vorticity is ``pushed'' around by the flow,
58 that is, $w(t,\varphi(t,x)) = w(0,x)$.
59
60 \section{The Biot-Savart Kernel}
61
62 Now, once we have the vorticity, we can recover the velocity $u$ by
63 solving the equation
64 \begin{eqnarray*}
65 \partial u_2/\partial x - \partial u_1/\partial y &=& w \\
66 \nabla \cdot u &=& 0 \\
67 u \cdot n &=& 0 \quad \text{on $\partial \Omega$}.
68 \end{eqnarray*}
69 This equation is solved by using a Biot-Savart kernel $K(x,y)$:
70 $$ u(x) = \int_\Omega K(x,y) w(y) \, dy .$$
71 The function $K$ depends upon the choice of domain.  First let us consider
72 the case when $\Omega$ is the whole plane (in which case the boundary
73 condition $u \cdot n = 0$ is replaced by saying that $u$ decays at infinity).
74 Then
75 \begin{equation*}
76 K(x,y) = K_1(x,y) = c \frac{(x-y)^\perp}{|x-y|^2} .
77 \end{equation*}
78 Here $x^\perp = (-x_2,x_1)$, and $c$ is a constant, probably something
79 like $1/2\pi$.  In any case we will set it to be one, which in effect
80 is rescaling the time variable, so we don't need to worry about it.
81
82 We can use this as a basis to find $K$ on the unit disk
83 $\Omega = \Delta = \{x:|x|<1\}$.  It turns out to be
84 \begin{equation*}
85 K_2(x,y) = K_1(x,y) - K_1(x,y^*) ,
86 \end{equation*}
87 where $y^* = y/|y|^2$ is called the reflection of $y$ about the
88 boundary of the unit disk.
89
90 Another example is if we have a bijective analytic function
91 $p:\Delta \to {\mathbb C}$, and we let $\Omega = p(\Delta)$.
92 (Here we think of $\Delta$ as a subset of $\mathbb C$, that is,
93 we are identifying the plane with the set of complex numbers.)
94 In that case we get
95 \[ K_p(p(x),p(y)) = K_2(x,y)/|p'(x)|^2 .\]
96 Our simulation considers the last case.  Examples of such
97 analytic functions include series 
98 $p(x) = x + \sum_{n=2}^\infty c_n x^n$, where
99 $\sum_{n=2}^\infty n |c_n| \le 1$.
100 (Thanks to David Ullrich for pointing this out to me.)
101
102 \section{The Simulation}
103
104 Now let's get to decribing the simulation.  We assume a rather
105 unusual initial distribution for the vorticity --- that the
106 vorticity is a finite sum of dirac delta masses.
107 \[ w(0,x) = \sum_{k=1}^N w_k \delta(x-x_k(0)) .\]
108 Here $x_k(0)$ is the initial place where the points
109 of vorticity are concentrated, with values $w_k$.  
110 Then at time $t$, the vorticity becomes
111 \[ w(t,x) = \sum_{k=1}^N w_k \delta(x-x_k(t)) .\]
112 The points of fluid $x_k(t)$ are pushed by the
113 flow, that is, $x_k(t) = \varphi(t,x_k(0))$, or
114 \[ \frac{\partial}{\partial t} x_k(t) = u(t,x_k(t)) .\]
115 Putting this all together, we finally obtain the equations
116 \[ \frac{\partial}{\partial t} x_k = \alpha_k \]
117 where
118 \[ \alpha_k   = \sum_{l=1}^N w_l K(x_k,x_l) .\]
119 This is the equation that our simulation solves.
120
121 In fact, in our case, where the domain is $p(\Delta)$,
122 the points are described by points
123 $\tilde x_k$, where $x_k = p(\tilde x_k)$.  Then
124 the equations become
125 \begin{eqnarray}
126 \label{tildex-p1}
127 \frac{\partial}{\partial t} \tilde x_k &=& \tilde\alpha_k \\
128 \label{tildex-p2}
129 \tilde\alpha_k &=& \frac1{|p'(\tilde x_k)|^2}
130      \sum_{l=1}^N w_l K_2(\tilde x_k,\tilde x_l) .
131 \end{eqnarray}
132
133 We solve this $2N$ system of equations using standard
134 numerical methods, in our case, using the second order midpoint method
135 for the first step, and thereafter using the second order Adams-Bashforth 
136 method.  (See for example the book
137 by Burden and Faires \cite{BF}).
138
139 \section{The Program - Data Structures}
140
141 The computer program solves equation (\ref{tildex-p1}), and displays
142 the results on the screen, with a boundary.  All the information
143 for solving the equation and displaying the output is countained
144 in the structure {\tt euler2dstruct}.  Let us describe some of
145 the fields in {\tt euler2dstruct}.  
146 The points $\tilde x_k$ are contained 
147 in {\tt double *x}: with the coordinates of
148 $\tilde x_k$ being the two numbers
149 {\tt x[2*k+0]}, {\tt x[2*k+1]}.  The values $w_k$ are contained
150 in {\tt double *w}.  The total number of points is
151 {\tt int N}.  (But only the first {\tt int Nvortex} points
152 have $w_k \ne 0$.)  The coefficients of the analytic function
153 (in our case a polynomial) $p$
154 are contained in {\tt double p\_coef[2*(deg\_p-1)]} --- here
155 {\tt deg\_p} is the degree of $p$, and the real and imaginary
156 parts of the coefficient
157 $c_n$ is contained in {\tt p\_coef[2*(n-2)+0]} and {\tt p\_coef[2*(n-2)+1]}.
158
159 \section{Data Initialization}
160
161 The program starts in the function {\tt init\_euler2d}.  After allocating
162 the memory for the data, and initialising some of the temporary variables
163 required for the numerical solving program, it randomly assigns the
164 coefficients of $p$, making sure that $\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$.
165 Then the program figures out how to draw the boundary, and what rescaling
166 of the data is required to draw it on the screen.  (This uses the
167 function {\tt calc\_p} which calculates $p(x)$.)
168
169 Next, it randomly assigns the initial values of $\tilde x_k$.  We want
170 to do this in such a way so that the points are uniformly spread over the
171 domain.  Let us first consider the case when the domain is the unit circle
172 $\Delta$.  In that case the proportion of points that we would expect
173 inside the circle of radius $r$ would be proportional to $r^2$.  So
174 we do it as follows:
175 \[ r = \sqrt{R_{0,1}},\quad \theta = R_{-\pi,\pi}, \quad
176    \tilde x_k = r (\cos \theta, \sin \theta) .\]
177 Here, and in the rest of this discussion, $R_{a,b}$ is a function
178 that returns a random variable uniformly distributed over the interval
179 $[a,b]$.
180
181 This works fine for $\Delta$, but for $p(\Delta)$, the points 
182 $p(\tilde x_k)$ are not uniformly distributed over $p(\Delta)$,
183 but are distributed with a density proportional to
184 $1/|p'(\tilde x_k)|^2$.  So to restore the uniform density we need
185 to reject this value of $\tilde x_k$ with probability proportional
186 to $|p'(\tilde x_k)|^2$.  Noticing that the condition 
187 $\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$ implies that 
188 $|p'(\tilde x_k)| \le 2$, we
189 do this by rejecting if $|p'(\tilde x_k)|^2 < R_{0,4}$.
190 (This makes use of the function {\tt calc\_mod\_dp2} which calculates
191 $|p'(x)|^2$.)
192
193 \section{Solving the Equation}
194
195 The main loop of the program is in the function {\tt draw\_euler2d}.
196 Most of the drawing operations are contained in this function, and
197 the numerical aspects are sent to the function {\tt ode\_solve}.
198 But there is an aspect of this that I would like
199 to discuss in the next section, and so we will look at a simple method for 
200 numerically solving differential equations.
201
202 The Euler Method
203 (nothing to do with the Euler Equation), is as
204 follows.  Pick a small number $h$ --- the time step (in
205 the program call {\tt delta\_t}).  Then we approximate
206 the solution of the equation:
207 \begin{equation}
208 \label{method-simple}
209 \tilde x_k(t+h) = \tilde x_k(t) + h \tilde\alpha_k(t) .
210 \end{equation}
211 The more sophisticated methods we use are variations of 
212 the Euler Method, and so the discussion in the following section
213 still applies.
214
215 In the program, the quantities $\tilde\alpha_k$, given by
216 equations (\ref{tildex-p2}) are calculated by the function
217 {\tt derivs}
218 (which in turns calls {\tt calc\_all\_mod\_dp2} to
219 calculate $|p'(\tilde x_k)|^2$ at all the points).
220
221
222 \section{Subtle Perturbation}
223
224 Added later: the scheme described here seems to not be that effective,
225 so now it is not used.
226
227 One problem using a numerical scheme such as the Euler Method occurs
228 when the points $\tilde x_k$ get close to the boundary
229 of $\Delta$.  In that case, it is possible that the new
230 points will be pushed outside of the boundary.  Even if they 
231 are not pushed out of the boundary, they may be much closer
232 or farther from the boundary than they should be.  
233 Our system of equations is very sensitive to how close points
234 are to the boundary --- points with non-zero vorticity
235 (``vortex points'') that are close to the boundary travel
236 at great speed alongside the boundary, with speed that is
237 inversely proportional to the distance from the boundary.
238
239 A way to try to mitigate this problem is something that I call
240 ``subtle perturbation.''
241 We map the points in 
242 the unit disk to points in the plane using the map
243 \begin{equation*}
244 F(x) = f(|x|) \frac x{|x|} ,
245 \end{equation*}
246 where $f:[0,1]\to[0,\infty]$ is an increasing continuous
247 bijection.  It turns out that a good choice is
248 \begin{equation*}
249 f(t) = -\log(1-t) .
250 \end{equation*}
251 (The reason for this is that points close to each other 
252 that are a distance
253 about $r$ from the boundary will be pushed around so that
254 their distance from each other is about multiplied by the
255 derivative of $\log r$, that is, $1/r$.)
256 Note that the inverse of this function is given by
257 \begin{equation*}
258 F^{-1}(x) = f^{-1}(|x|) \frac x{|x|} ,
259 \end{equation*}
260 where 
261 \begin{equation*}
262 f^{-1}(t) = 1-e^{-t} .
263 \end{equation*}
264
265 So what we could do is the following: instead of working with
266 the points $\tilde x_k$, we could work instead with the points
267 $y_k = F(\tilde x_k)$.  In effect this is what we do.
268 Instead of performing the computation (\ref{method-simple}),
269 we do the calculation
270 \begin{equation*}
271 y_k = F(\tilde x_k(t)) + h {\cal A}(\tilde x_k) \tilde\alpha_k(t) 
272 \end{equation*}
273 where
274 ${\cal A}(x)$ is the matrix of partial derivatives of $F$:
275 \begin{equation*}
276 {\cal A}(x) = 
277 \frac{f(|x|)}{|x|}
278 \left[
279 \begin{matrix}
280 1 & 0\\
281 0 & 1
282 \end{matrix}
283 \right]
284 + \frac1{|x|}
285   \left(\frac{f'(|x|)}{|x|} - \frac{f(|x|)}{|x|^2}\right)
286 \left[
287 \begin{matrix}
288 x_{1}^2   & x_{1} x_{2}\\
289 x_{1} x_{2} & x_{2}^2
290 \end{matrix}
291 \right],
292 \end{equation*}
293 and then compute
294 \begin{equation*}
295 \tilde x_k(t+h) = F^{-1}(y_k).
296 \end{equation*}
297 These calculations are done in the function {\tt perturb}, if
298 the quantity {\tt SUBTLE\_PERTURB} is set.
299
300 \section{Drawing the Points}
301
302 As we stated earlier, most of the drawing functions are contained
303 in the function {\tt draw\_euler2d}.  If the variable 
304 {\tt hide\_vortex} is set (and the function {\tt init\_euler2d}
305 will set this with probability $3/4$), then we only display
306 the points $\tilde x_k$ for ${\tt Nvortex} < k \le N$.  If 
307 {\tt hide\_vortex} is not set, then the ``vortex points''
308 $\tilde x_k$ ($1 \le k \le {\tt Nvortex}$) are displayed in white.
309 In fact the points $p(\tilde x_k)$ are what are put onto the screen,
310 and for this we make use of the function {\tt calc\_all\_p}.
311
312 \section{Addition to Program: Changing the Power Law}
313
314 A later addition to the program adds an option {\tt eulerpower},
315 which allows one to change the power law that describes how
316 the vortex points influence other points.  In effect, if this
317 option is set with the value $m$, then the Biot-Savart Kernel
318 is replace by
319 $$ K_1(x,y) = \frac{(x-y)^\perp}{|x-y|^{m+1}}, $$
320 and
321 $$ K_2(x,y) = K_1(x,y) - |y|^{1-m} K_1(x,y) .$$
322 So for example, setting $m=2$ corresponds to the 
323 quasi-geostrophic equation.  (I haven't yet figured out
324 what $K_p$ should be, so if $m \ne 1$ we use the unit circle
325 as the boundary.)
326
327 \begin{thebibliography}{9}
328
329 \bibitem{BF} Richard L. Burden, J. Douglas Faires, Numerical Analysis,
330 sixth edition, Brooks/Cole, 1996.
331
332 \bibitem{C} Alexandre J. Chorin, Vorticity and Turbulence,
333 Applied Mathematical Sciences, Vol 103, Springer Verlag, 1994.
334
335 \end{thebibliography}
336
337 \end{document}