\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}