1 files changed, 337 insertions, 0 deletions
diff --git a/hacks/euler2d.tex b/hacks/euler2d.tex
new file mode 100644
index 0000000..700ad31
--- /dev/null
+++ b/hacks/euler2d.tex
@@ -0,0 +1,337 @@
+\documentclass[12pt]{article}
+
+%\usepackage{fullpage}
+\usepackage{amsmath,amssymb}
+
+\begin{document}
+
+\title{Two Dimensional Euler Simulation}
+
+\author{
+S.J. Montgomery-Smith\\
+Department of Mathematics\\
+University of Missouri\\
+Columbia, MO 65211, U.S.A.\\
+stephen@math.missouri.edu\\
+http://www.math.missouri.edu/\~{}stephen}
+
+\date{September 10, 2000}
+
+\maketitle
+
+This document describes a program I wrote to simulate the 
+two dimensional Euler Equation --- a program that is part
+of the {\tt xlock} screensaver as the {\tt euler2d}
+mode.  A similar explanation may also be found in the
+book by Chorin \cite{C}.
+
+\section{The Euler Equation}
+
+The Euler Equation describes the motion of an incompressible
+fluid that has no viscosity.  If the fluid is contained
+in a domain $\Omega$ with boundary $\partial \Omega$, then
+the equation is in the vector field $u$ (the velocity)
+and the
+scalar field $p$ (the pressure):
+\begin{eqnarray*}
+\frac{\partial}{\partial t} u &=& -u \cdot \nabla u + \nabla p \\
+\nabla \cdot u &=& 0 \\
+u \cdot n &=& 0 \quad \text{on $\partial \Omega$}
+\end{eqnarray*}
+where $n$ is the unit normal to $\partial \Omega$.
+
+\section{Vorticity}
+
+It turns out that it can be easier write these equations
+in terms of the vorticity.  In two dimensions the vorticity
+is the scalar $w = \partial u_2/\partial x - \partial u_1/\partial y$.
+The equation for vorticity becomes
+\[ \frac{\partial}{\partial t} w = -u \cdot \nabla w .\]
+A solution to this equation can be written as follows.  The velocity
+$u$ causes a flow, that is, a function $\varphi(t,x)$ that tells where
+the particle initially at $x$ ends up at time $t$, that is
+\[
+\frac\partial{\partial t} \varphi(t,x)
+= u(t,\varphi(t,x)) .\]
+Then the equation
+for $w$ tells us that the vorticity is ``pushed'' around by the flow,
+that is, $w(t,\varphi(t,x)) = w(0,x)$.
+
+\section{The Biot-Savart Kernel}
+
+Now, once we have the vorticity, we can recover the velocity $u$ by
+solving the equation
+\begin{eqnarray*}
+\partial u_2/\partial x - \partial u_1/\partial y &=& w \\
+\nabla \cdot u &=& 0 \\
+u \cdot n &=& 0 \quad \text{on $\partial \Omega$}.
+\end{eqnarray*}
+This equation is solved by using a Biot-Savart kernel $K(x,y)$:
+$$ u(x) = \int_\Omega K(x,y) w(y) \, dy .$$
+The function $K$ depends upon the choice of domain.  First let us consider
+the case when $\Omega$ is the whole plane (in which case the boundary
+condition $u \cdot n = 0$ is replaced by saying that $u$ decays at infinity).
+Then
+\begin{equation*}
+K(x,y) = K_1(x,y) = c \frac{(x-y)^\perp}{|x-y|^2} .
+\end{equation*}
+Here $x^\perp = (-x_2,x_1)$, and $c$ is a constant, probably something
+like $1/2\pi$.  In any case we will set it to be one, which in effect
+is rescaling the time variable, so we don't need to worry about it.
+
+We can use this as a basis to find $K$ on the unit disk
+$\Omega = \Delta = \{x:|x|<1\}$.  It turns out to be
+\begin{equation*}
+K_2(x,y) = K_1(x,y) - K_1(x,y^*) ,
+\end{equation*}
+where $y^* = y/|y|^2$ is called the reflection of $y$ about the
+boundary of the unit disk.
+
+Another example is if we have a bijective analytic function
+$p:\Delta \to {\mathbb C}$, and we let $\Omega = p(\Delta)$.
+(Here we think of $\Delta$ as a subset of $\mathbb C$, that is,
+we are identifying the plane with the set of complex numbers.)
+In that case we get
+\[ K_p(p(x),p(y)) = K_2(x,y)/|p'(x)|^2 .\]
+Our simulation considers the last case.  Examples of such
+analytic functions include series 
+$p(x) = x + \sum_{n=2}^\infty c_n x^n$, where
+$\sum_{n=2}^\infty n |c_n| \le 1$.
+(Thanks to David Ullrich for pointing this out to me.)
+
+\section{The Simulation}
+
+Now let's get to decribing the simulation.  We assume a rather
+unusual initial distribution for the vorticity --- that the
+vorticity is a finite sum of dirac delta masses.
+\[ w(0,x) = \sum_{k=1}^N w_k \delta(x-x_k(0)) .\]
+Here $x_k(0)$ is the initial place where the points
+of vorticity are concentrated, with values $w_k$.  
+Then at time $t$, the vorticity becomes
+\[ w(t,x) = \sum_{k=1}^N w_k \delta(x-x_k(t)) .\]
+The points of fluid $x_k(t)$ are pushed by the
+flow, that is, $x_k(t) = \varphi(t,x_k(0))$, or
+\[ \frac{\partial}{\partial t} x_k(t) = u(t,x_k(t)) .\]
+Putting this all together, we finally obtain the equations
+\[ \frac{\partial}{\partial t} x_k = \alpha_k \]
+where
+\[ \alpha_k   = \sum_{l=1}^N w_l K(x_k,x_l) .\]
+This is the equation that our simulation solves.
+
+In fact, in our case, where the domain is $p(\Delta)$,
+the points are described by points
+$\tilde x_k$, where $x_k = p(\tilde x_k)$.  Then
+the equations become
+\begin{eqnarray}
+\label{tildex-p1}
+\frac{\partial}{\partial t} \tilde x_k &=& \tilde\alpha_k \\
+\label{tildex-p2}
+\tilde\alpha_k &=& \frac1{|p'(\tilde x_k)|^2}
+     \sum_{l=1}^N w_l K_2(\tilde x_k,\tilde x_l) .
+\end{eqnarray}
+
+We solve this $2N$ system of equations using standard
+numerical methods, in our case, using the second order midpoint method
+for the first step, and thereafter using the second order Adams-Bashforth 
+method.  (See for example the book
+by Burden and Faires \cite{BF}).
+
+\section{The Program - Data Structures}
+
+The computer program solves equation (\ref{tildex-p1}), and displays
+the results on the screen, with a boundary.  All the information
+for solving the equation and displaying the output is countained
+in the structure {\tt euler2dstruct}.  Let us describe some of
+the fields in {\tt euler2dstruct}.  
+The points $\tilde x_k$ are contained 
+in {\tt double *x}: with the coordinates of
+$\tilde x_k$ being the two numbers
+{\tt x[2*k+0]}, {\tt x[2*k+1]}.  The values $w_k$ are contained
+in {\tt double *w}.  The total number of points is
+{\tt int N}.  (But only the first {\tt int Nvortex} points
+have $w_k \ne 0$.)  The coefficients of the analytic function
+(in our case a polynomial) $p$
+are contained in {\tt double p\_coef[2*(deg\_p-1)]} --- here
+{\tt deg\_p} is the degree of $p$, and the real and imaginary
+parts of the coefficient
+$c_n$ is contained in {\tt p\_coef[2*(n-2)+0]} and {\tt p\_coef[2*(n-2)+1]}.
+
+\section{Data Initialization}
+
+The program starts in the function {\tt init\_euler2d}.  After allocating
+the memory for the data, and initialising some of the temporary variables
+required for the numerical solving program, it randomly assigns the
+coefficients of $p$, making sure that $\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$.
+Then the program figures out how to draw the boundary, and what rescaling
+of the data is required to draw it on the screen.  (This uses the
+function {\tt calc\_p} which calculates $p(x)$.)
+
+Next, it randomly assigns the initial values of $\tilde x_k$.  We want
+to do this in such a way so that the points are uniformly spread over the
+domain.  Let us first consider the case when the domain is the unit circle
+$\Delta$.  In that case the proportion of points that we would expect
+inside the circle of radius $r$ would be proportional to $r^2$.  So
+we do it as follows:
+\[ r = \sqrt{R_{0,1}},\quad \theta = R_{-\pi,\pi}, \quad
+   \tilde x_k = r (\cos \theta, \sin \theta) .\]
+Here, and in the rest of this discussion, $R_{a,b}$ is a function
+that returns a random variable uniformly distributed over the interval
+$[a,b]$.
+
+This works fine for $\Delta$, but for $p(\Delta)$, the points 
+$p(\tilde x_k)$ are not uniformly distributed over $p(\Delta)$,
+but are distributed with a density proportional to
+$1/|p'(\tilde x_k)|^2$.  So to restore the uniform density we need
+to reject this value of $\tilde x_k$ with probability proportional
+to $|p'(\tilde x_k)|^2$.  Noticing that the condition 
+$\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$ implies that 
+$|p'(\tilde x_k)| \le 2$, we
+do this by rejecting if $|p'(\tilde x_k)|^2 < R_{0,4}$.
+(This makes use of the function {\tt calc\_mod\_dp2} which calculates
+$|p'(x)|^2$.)
+
+\section{Solving the Equation}
+
+The main loop of the program is in the function {\tt draw\_euler2d}.
+Most of the drawing operations are contained in this function, and
+the numerical aspects are sent to the function {\tt ode\_solve}.
+But there is an aspect of this that I would like
+to discuss in the next section, and so we will look at a simple method for 
+numerically solving differential equations.
+
+The Euler Method
+(nothing to do with the Euler Equation), is as
+follows.  Pick a small number $h$ --- the time step (in
+the program call {\tt delta\_t}).  Then we approximate
+the solution of the equation:
+\begin{equation}
+\label{method-simple}
+\tilde x_k(t+h) = \tilde x_k(t) + h \tilde\alpha_k(t) .
+\end{equation}
+The more sophisticated methods we use are variations of 
+the Euler Method, and so the discussion in the following section
+still applies.
+
+In the program, the quantities $\tilde\alpha_k$, given by
+equations (\ref{tildex-p2}) are calculated by the function
+{\tt derivs}
+(which in turns calls {\tt calc\_all\_mod\_dp2} to
+calculate $|p'(\tilde x_k)|^2$ at all the points).
+
+
+\section{Subtle Perturbation}
+
+Added later: the scheme described here seems to not be that effective,
+so now it is not used.
+
+One problem using a numerical scheme such as the Euler Method occurs
+when the points $\tilde x_k$ get close to the boundary
+of $\Delta$.  In that case, it is possible that the new
+points will be pushed outside of the boundary.  Even if they 
+are not pushed out of the boundary, they may be much closer
+or farther from the boundary than they should be.  
+Our system of equations is very sensitive to how close points
+are to the boundary --- points with non-zero vorticity
+(``vortex points'') that are close to the boundary travel
+at great speed alongside the boundary, with speed that is
+inversely proportional to the distance from the boundary.
+
+A way to try to mitigate this problem is something that I call
+``subtle perturbation.''
+We map the points in 
+the unit disk to points in the plane using the map
+\begin{equation*}
+F(x) = f(|x|) \frac x{|x|} ,
+\end{equation*}
+where $f:[0,1]\to[0,\infty]$ is an increasing continuous
+bijection.  It turns out that a good choice is
+\begin{equation*}
+f(t) = -\log(1-t) .
+\end{equation*}
+(The reason for this is that points close to each other 
+that are a distance
+about $r$ from the boundary will be pushed around so that
+their distance from each other is about multiplied by the
+derivative of $\log r$, that is, $1/r$.)
+Note that the inverse of this function is given by
+\begin{equation*}
+F^{-1}(x) = f^{-1}(|x|) \frac x{|x|} ,
+\end{equation*}
+where 
+\begin{equation*}
+f^{-1}(t) = 1-e^{-t} .
+\end{equation*}
+
+So what we could do is the following: instead of working with
+the points $\tilde x_k$, we could work instead with the points
+$y_k = F(\tilde x_k)$.  In effect this is what we do.
+Instead of performing the computation (\ref{method-simple}),
+we do the calculation
+\begin{equation*}
+y_k = F(\tilde x_k(t)) + h {\cal A}(\tilde x_k) \tilde\alpha_k(t) 
+\end{equation*}
+where
+${\cal A}(x)$ is the matrix of partial derivatives of $F$:
+\begin{equation*}
+{\cal A}(x) = 
+\frac{f(|x|)}{|x|}
+\left[
+\begin{matrix}
+1 & 0\\
+0 & 1
+\end{matrix}
+\right]
++ \frac1{|x|}
+  \left(\frac{f'(|x|)}{|x|} - \frac{f(|x|)}{|x|^2}\right)
+\left[
+\begin{matrix}
+x_{1}^2   & x_{1} x_{2}\\
+x_{1} x_{2} & x_{2}^2
+\end{matrix}
+\right],
+\end{equation*}
+and then compute
+\begin{equation*}
+\tilde x_k(t+h) = F^{-1}(y_k).
+\end{equation*}
+These calculations are done in the function {\tt perturb}, if
+the quantity {\tt SUBTLE\_PERTURB} is set.
+
+\section{Drawing the Points}
+
+As we stated earlier, most of the drawing functions are contained
+in the function {\tt draw\_euler2d}.  If the variable 
+{\tt hide\_vortex} is set (and the function {\tt init\_euler2d}
+will set this with probability $3/4$), then we only display
+the points $\tilde x_k$ for ${\tt Nvortex} < k \le N$.  If 
+{\tt hide\_vortex} is not set, then the ``vortex points''
+$\tilde x_k$ ($1 \le k \le {\tt Nvortex}$) are displayed in white.
+In fact the points $p(\tilde x_k)$ are what are put onto the screen,
+and for this we make use of the function {\tt calc\_all\_p}.
+
+\section{Addition to Program: Changing the Power Law}
+
+A later addition to the program adds an option {\tt eulerpower},
+which allows one to change the power law that describes how
+the vortex points influence other points.  In effect, if this
+option is set with the value $m$, then the Biot-Savart Kernel
+is replace by
+$$ K_1(x,y) = \frac{(x-y)^\perp}{|x-y|^{m+1}}, $$
+and
+$$ K_2(x,y) = K_1(x,y) - |y|^{1-m} K_1(x,y) .$$
+So for example, setting $m=2$ corresponds to the 
+quasi-geostrophic equation.  (I haven't yet figured out
+what $K_p$ should be, so if $m \ne 1$ we use the unit circle
+as the boundary.)
+
+\begin{thebibliography}{9}
+
+\bibitem{BF} Richard L. Burden, J. Douglas Faires, Numerical Analysis,
+sixth edition, Brooks/Cole, 1996.
+
+\bibitem{C} Alexandre J. Chorin, Vorticity and Turbulence,
+Applied Mathematical Sciences, Vol 103, Springer Verlag, 1994.
+
+\end{thebibliography}
+
+\end{document}