\documentclass[reqno]{amsart}
\usepackage{graphicx}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
Sixth Mississippi State Conference on Differential Equations and
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conference 15 (2007),  pp. 193--210.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2007 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{193}
\title[\hfilneg EJDE-2006/Conf/15\hfil An augmented IIM-level set method]
{An augmented IIM-level set method for Stokes
equations with discontinuous viscosity}

\author[Z. Li, S. R. Lubkin, X. Wan\hfil EJDE/Conf/15 \hfilneg]
{Zhilin Li, Sharon R. Lubkin, Xiaohai Wan}  % in alphabetical order

\address{Zhilin Li \newline
Center for Research in
Scientific Computation \& Department of Mathematics, North
Carolina State University, Raleigh, NC 27695-8205, USA}
\email{zhilin@math.ncsu.edu}

\address{Sharon R. Lubkin \newline
Center for Research in Scientific Computation \& Department of Mathematics,
North Carolina State University, Raleigh, NC 27695-8205, USA}
\email{lubkin@math.ncsu.edu}

\address{Xiaohai Wan\newline
Biomathematics Graduate Program \& Department of Mathematics,
North Carolina State University, Raleigh, NC 27695-8203, USA}
\email{xwan@unity.ncsu.edu}

\thanks{Published February 28, 2007.}
\subjclass[2000]{65M06, 65M12, 76T05}
\keywords{Incompressible Stokes equations; discontinuous viscosity;
\hfill\break\indent
 interface problem; immersed interface method; level set method;
 fast Poisson solver, GMRES}

\begin{abstract}
 A finite difference method is proposed for
 incompressible Stokes equations with discontinuous viscosity.
 The method combines the augmented immersed interface method with a
 level set representation of the interface on a uniform Cartesian grid.
 In the proposed method, augmented interface variables are introduced
 along the interface and be solved  by the GMRES iterative method.
 The augmented approach allows fast Poisson
 solvers to be used and therefore is very efficient.  Numerical examples
 including two moving interface problems are presented.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}


\section{Introduction}

In this paper, we present a numerical method for solving
the stationary incompressible Stokes equations for two-phase Newtonian
flows with an arbitrary interface. The governing equations have the
following form
\begin{equation}\label{eq:1}
\nabla p = \nabla \cdot \mu \left( \nabla \mathbf{u} +
(\nabla \mathbf{u})^T \right) + \mathbf{F} (\mathbf{x}) + \mathbf g(\mathbf{x}),
\end{equation}
with the incompressibility condition
\begin{equation} \label{incomp}
\nabla \cdot \mathbf{u} = 0,
\end{equation}
where $\mathbf{u} =(u,v)^T$ is the fluid velocity, $p$ is the fluid pressure,
$\mu$ is the fluid viscosity, $\mathbf{x} = (x,y)$ is the Cartesian coordinate
variable, $\mathbf g = (g_1,g_2)^T$ is a body force
such as gravity, and $\mathbf{F}$ is a source which can have a Dirac delta
function singularity,
\begin{equation}
\mathbf{F} (\mathbf{x}) = \int_{\Gamma} \mathbf f(s)
\delta(\mathbf{x}-\mathbf{X}(s))ds
\end{equation}
where $\Gamma$ is a smooth interface in the solution domain with the
arc-length parameter $s$,  $\mathbf f = (f_1,f_2)^T$ is the source density,
and $\delta(\cdot)$ is the
Dirac delta function defined in the distribution sense. The body force
term $\mathbf g$ may also have a finite
jump across the interface $\Gamma$ as well.

The interface $\Gamma$ divides the solution domain $\Omega$ into two parts
$\Omega^+$ and $\Omega^-$ with $\Omega = \Omega^+ \cup \Omega^-$. If the
interface is closed, we use $\Omega^+$ to express the exterior domain of
the interface.

We assume that the viscosity coefficient $\mu$ is
piecewise constant across the interface:
\begin{equation} \label{eqmu}
\mu (\mathbf{x})= \begin {cases}
\mu^+  &\mbox{if }\mathbf{x}\in \Omega^+, \\
 \mu^- & \mbox{if }\mathbf{x}\in \Omega^-,
\end {cases}
\end{equation}
where $\mu^+$ and $\mu^-$ are two positive constants.

In this paper, we propose a finite difference method
coupled with a level set function representation of the interface
on a uniform Cartesian grid to solve the problem. Finite element and
boundary integral methods for this problem can be found in the literature.
These methods may be more efficient if only the solution on the interface
is solved. In our method, we solve for the velocity in the entire domain.

There are at least two difficulties in solving \eqref{eq:1}-\eqref{eqmu}
numerically using finite difference methods. The first one is to deal with
the singular force term. A
simple way is to use Peskin's discrete delta function approach \cite{peskin1} to distribute
the singular force to nearby grid points. Such a discretization is
typically first
order accurate and will smooth out the solution. The second difficulty is how
to deal with the discontinuous viscosity. A simple smoothing method  may
introduce large errors; see \cite{li:review} for a one-dimensional example.

In the case of a continuous viscosity with a Dirac delta function
source distribution along an interface, various methods have been developed. We refer the readers to
\cite{cortez-stokes,fo-pe,rjl-li:stokes,mayo-pe,tu-peskin} for various
methods and the references therein. The difficulty with a discontinuous
viscosity is that the jump conditions for the pressure and velocity are
coupled together (section \ref{sec:4}),
 which makes it difficult to
discretize the system accurately. Approximations of the jump conditions have been used to decouple the system
\cite{li-lubkin} and have been productively used to simulate problems in biological fluid dynamics \cite{lubkin-li}.

 The augmented immersed interface
method for incompressible 2D Stokes flows with discontinuous viscosity
has been developed \cite{li-ito-lai} with the interface represented by a
cubic spline interpolation, a variation of the particle approach.
It is easy with splines to obtain the surface derivatives of jumps of the material
variables (pressure and velocity) since they can be expressed as the functions
of the arc-length. However, the spline approach is difficult to use
for multi-connected domains, moving interface problems with topological
changes such as merging and splitting, and for three dimensional (3D)
problems. In addition, the level set method usually has better stability. For a
spline approach, re-parameterization, filtering, and re-griding may be
needed at every or every other time steps. All these reasons favor a level set approach over a spline approach.
Nevertheless, it is not trivial
to compute surface derivatives of an interface quantity from the level set
method since the interface is explicitly defined only at grid points.


\section{The Stokes solver: Three Poisson equations approach} \label{sec:3}

Instead of writing the singular force as a source term  for the Stokes
equations (\ref{eq:1}), we can solve the Stokes equations
in the form
\begin{equation}\label{stokes1}
\nabla p = \nabla \cdot \mu^+\left( \nabla \mathbf{u} +
(\nabla \mathbf{u})^T\right) + \mathbf g(\mathbf{x}),
\quad \mathbf{x} \in \Omega^+
\end{equation}
and
\begin{equation}\label{stokes2}
\nabla p = \nabla \cdot \mu^- \left( \nabla \mathbf{u} +
(\nabla \mathbf{u})^T \right) + \mathbf g(\mathbf{x}),
\quad \mathbf{x} \in \Omega^-.
\end{equation}
We define viscosity ratio $ \lambda = {\mu^-}/{\mu^+}$.
If $\Gamma$ is closed and $\lambda > 1$, then
the fluid phase inside the interface is more viscous and is called a drop;
if $\Gamma$ is closed and $\lambda < 1$, then the fluid phase inside the
interface is less viscous than the ambient fluid and is called a bubble
(see for example \cite{POZRIKIDIS-eds}); if $\lambda = 1$ the fluids inside
and outside have equal viscosity.

The solutions to (\ref{stokes1})-(\ref{stokes2}) are coupled by the jump
conditions which will be stated
in Section \ref{sec:4}.
If we apply the divergence operator to both sides of the momentum equation
(\ref{stokes1}, \ref{stokes2}) and make use of the incompressibility
condition(\ref{incomp}) in each fluid, we  get
\begin{equation} \label{p1}
\Delta p = \nabla \cdot \mathbf g,
\end{equation}
where $\Delta$ is the Laplacian.
This is a Poisson equation for the pressure $p$. Once $p$ is solved from the
equation above, we can rewrite the Stokes equation
(\ref{stokes1})-(\ref{stokes2}) in the form \cite{chorin:proj}
\begin{gather} \label{p2}
\Delta \tilde{u} = \frac{\partial p}{\partial x} - g_1, \\
 \label{p3}
\Delta \tilde{v} = \frac{\partial p}{\partial y} - g_2,
\end{gather}
where
\begin{equation}
\tilde{u} ( {\mathbf{x}}) = \begin{cases} \mu^+ u( {\mathbf{x}})
\quad {\mathbf{x}} \in \Omega^+  \\
\mu^- u( {\mathbf{x}}) \quad {\mathbf{x}} \in \Omega^-, \end{cases}  \quad
\tilde{v}( {\mathbf{x}}) = \begin{cases} \mu^+ v ( {\mathbf{x}})\quad \mathbf{x}
\in \Omega^+  \\
\mu^- v( {\mathbf{x}}) \quad \mathbf{x} \in \Omega^-. \end{cases}
\end{equation}
We obtain $\mathbf{u} = (u,v)$ by solving $\tilde{\mathbf{u}} = (\tilde{u},
\tilde{v})$ from the Poisson equations above.
The three decoupled Poisson equations (\ref{p1})-(\ref{p3}) each can be
solved efficiently using a  fast Poisson solver if the solutions are smooth
across the interface.

\section{Jump conditions} \label{sec:4}

It is well known that if there are singular
interface forces, then the pressure $p$ is discontinuous, and the velocity
${\mathbf{u}}$ is non-smooth across the interface. If the viscosity is discontinuous, that is
$\lambda \neq 1$, then $\tilde{\mathbf{u}}$ is also discontinuous.
Due to these discontinuities, the direct application of three Poisson
solvers without modifications will not give the correct solution.
 It is critical to incorporate the jump conditions in the
 finite difference schemes so that the solutions in different phases
can be coupled together correctly across the interface. The jump condition on
the pressure $p$ at point $\mathbf{x}^*$ on the interface $\Gamma$ is defined as
\begin{equation}
[p]({\mathbf{x}^{*}}) = \Big(\lim_{\mathbf{x} \to
(\mathbf{x}^{*})^+}{p}({\mathbf{x}})\Big) -
 \Big(\lim_{\mathbf{x} \to (\mathbf{x}^{*})^-}{p}({\mathbf{x}}) \Big).
\end{equation}
Other jump conditions are defined similarly. Jump conditions on each variable
can be derived either from the physics of the problem itself or from the
partial differential equations.  The derivations of the jump conditions for
the Stokes equation are described elsewhere \cite{li-ito-lai}. Here we just list
the results. Define
\begin{equation}
\mathbf q = \left( q_1, q_2 \right)^T = \left(  \left[ \tilde{u} \right] ,  \left[ \tilde{v} \right] \right)^T.
\end{equation}
Then the solutions of the Stokes flow are the following three Poisson
equations:
\begin{equation} %\label{S2}
\left\{ \begin{aligned}
   & \Delta p = \nabla \cdot {\bf g},  \\
   & [p] = \hat{f}_1 - 2 \frac{\partial}{\partial {\mathbf{\tau}}}
\left( \mathbf{q} \cdot \mathbf{\mathbf{\tau}}\left( \mathbf{x}^{*}\right) \right), \\
 &\big[\frac{\partial p}{\partial {\mathbf{n}}}\big] = \frac{\partial \hat{f}_2}{\partial{\mathbf{\tau}}} +
      2\frac{\partial^2}{\partial {\mathbf{\tau}}^2} \left( \mathbf{q} \cdot
      {\mathbf{n}}\left( \mathbf{x}^{*} \right) \right) -
4 \kappa \frac{\partial}{\partial {\mathbf{\tau}}}( {\bf q}\cdot {\mathbf{\tau}}\left( \mathbf{x}^{*} \right))
 +  [ \mathbf g \cdot \mathbf{n}  ],
\end{aligned} \right.\label{Sab}
\end{equation}
\begin{equation}
\left\{\begin{aligned}
& \Delta \tilde{u} = p_x - g_1, \\
& [\tilde{u}] =q_1, \\
& \big[\frac{\partial \tilde{u}}{\partial {\mathbf{n}}}\big ]
= \Big( \hat{f}_2 +  \frac{\partial }{\partial {\mathbf{\tau}}}
          \left( {\mathbf{q}} \cdot {\mathbf{n}}\left( \mathbf{x}^{*} \right) \right) \Big)
           \sin \theta -
          \Big(  \frac{\partial }{\partial {\mathbf{\tau}}} \left( {\mathbf{q}} \cdot
         {\mathbf{\tau}}\left( \mathbf{x}^{*}  \right) \right) \Big) \cos \theta,  \\
    \end{aligned} \right. \label{Sbb}
\end{equation}
\begin{equation}
\left\{ \begin{aligned}
 & \Delta \tilde{v} = p_y-g_2, \\
 & [\tilde{v}] =q_2, \\
 & \big[\frac{\partial \tilde{v} }{\partial {\mathbf{n}}}
\big] = -\Big( \hat{f}_2 +   \frac{\partial }{\partial {\mathbf{\tau}}}
          \left( {\mathbf{q}} \cdot {\mathbf{n}}\left( \mathbf{x}^{*} \right) \right) \Big)
           \cos \theta -
          \Big(  \frac{\partial }{\partial {\mathbf{\tau}}} \left( {\mathbf{q}} \cdot
         {\mathbf{\mathbf{\tau}}}\left( \mathbf{x}^{*} \right) \right) \Big)
          \sin \theta,
     \end{aligned} \right. \label{Scb}
\end{equation}
\begin{equation}
 \big[ \frac{\tilde{u}}{\mu} \big] = 0, \quad
 \big[ \frac{\tilde{v}}{\mu} \big] = 0. \label{Scd}
\end{equation}
In the jump conditions above, $\mathbf{n} = (\cos\theta,\sin\theta)^T$ and
$\mathbf{\tau} = (-\sin\theta,\cos\theta)^T$ are unit outward normal and
tangential
directions respectively, $\theta$ is the angle between the $x$-axis and the
normal directions, $\hat{f_{1}} = \mathbf f \cdot \mathbf{n}$ and
$\hat{f_{2}} = \mathbf f \cdot \mathbf \mathbf{\tau}$ are the force components in
the normal and tangential directions. $\kappa$ is interface curvature. Unit circle
is defined to have constant curvature $-1$.
The vector $\mathbf q$ is called
the augmented variable which is
unknown unless we know the solution $\left( u,v \right)^T$. However, if we know
$\mathbf q$, then
the jump conditions for the pressure and the velocity are decoupled in
terms of $\mathbf q$ and its surface derivatives.
In this paper,
we assume that we have a periodic boundary for all the variables.

The main idea of the
augmented approach is to solve $\mathbf q$ such
that $(\tilde{u},\tilde{v})^T$
and $p$ satisfy the Stokes equations, the jump conditions above, and
the kinematic condition
\begin{equation}
[u] = [v] = 0.
\end{equation}
Then $(u,v)^T $ and $p$ must be the solution to the original problem.

The jump conditions above enable us to use the immersed interface
method (IIM) to solve the three Poisson equations (\ref{p1})-(\ref{p3}) if
a guess $\mathbf q$ is given. The generalized minimum residual (GMRES)
method \cite{saad} is used to solve for
$\mathbf q$ so that the kinematic condition can be satisfied.

\section{The numerical algorithm} \label{sec:5}

\subsection{Level set representation of interfaces and related quantities}

We use the zero level set of a two dimensional function $\varphi(x,y)$ to
represent the interface. That is, $\varphi(x,y)$ is such a function that
\begin{equation}
\Gamma = \left \{ (x,y): \varphi(x,y) = 0 \right \}.
\end{equation}
Generally $\varphi(x,y)$ is chosen as the signed distance function from the
interface.
We denote $\Omega^+ = \{(x,y), \, |\, \varphi(x,y) \geq 0\}$ and
 $\Omega^- = \{(x,y), \, |\, \varphi(x,y) <0\}$. For simplicity, we assume
 that the domain $\Omega$ is a rectangle $[a,\;b]\times [c,\;d]$.
The solution domain is discretized using a uniform Cartesian grid
\begin{equation}
\begin{aligned}
x_i = a + ih_x, \quad i = 0,1,\dots,m, \quad h_x = \frac{b-a}{m}; \\
y_j = c + jh_y, \quad j = 0,1,\dots,n, \quad h_y = \frac{d-c}{n}.
\end{aligned}
\end{equation}
In the rest of the paper,  we take
$h_x = h_y = h$ to simplify the notation. The level set function
$\varphi_{i,j}=\varphi(x_i,y_j)$ is defined at grid points.
A grid point $(x_i,y_j)$ is an \emph{irregular inner grid point} if
\begin{equation}
\varphi_{i,j} <0
\end{equation}
and any of the following four inequalities is true:
\begin{equation}
\begin{gathered}
\varphi_{i,j}\cdot \varphi_{i-1,j} \leq 0, \quad
\varphi_{i,j}\cdot \varphi_{i+1,j} \leq 0, \\
\varphi_{i,j}\cdot \varphi_{i,j-1} \leq 0, \quad
\varphi_{i,j}\cdot \varphi_{i,j+1} \leq 0.
\end{gathered}
\end{equation}
Similarly, a grid point $(x_i,y_j)$ is an \emph{irregular outer grid point} if
\begin{equation}
\varphi_{i,j} \geq 0
\end{equation}
and any of the following four inequalities is true:
\begin{equation}
\begin{gathered}
\varphi_{i,j}\cdot \varphi_{i-1,j} < 0, \quad
\varphi_{i,j}\cdot \varphi_{i+1,j} < 0, \\
\varphi_{i,j}\cdot \varphi_{i,j-1} < 0, \quad
\varphi_{i,j}\cdot \varphi_{i,j+1} < 0.
\end{gathered}
\end{equation}
For a given irregular grid point $\mathbf{x} = (x_i,y_j)$ near a
sufficiently flat section of interface,
there is a corresponding orthogonal projection point
$\mathbf{x}^* = (x^*,y^*)$ on the interface satisfying
\begin{equation}
\mathbf{x}^* = \mathbf{x} + a \mathbf{n}(\mathbf{x})
\end{equation}
where $a$ is a real unknown variable and
$\mathbf{n} = {\nabla \varphi}/{|\nabla \varphi|}$ is the unit outward normal
direction at $\mathbf{x}$. See Fig. \ref{f1} for an illustration.

\begin{figure}[hpbt]
\begin{center}
\includegraphics[width=0.6\textwidth]{figures/x1}
\end{center}
\caption{An irregular grid point $\mathbf{x} = (x_i,y_j)$ and
its normal projection point $\mathbf{x}^* = (x^*,y^*)$ on
the interface $\Gamma$. $(\xi,\eta)$ is the local coordinate system at
$\mathbf{x}^*$.}
\label{f1}
\end{figure}

Since $\varphi(\mathbf{x}^*) = 0$,  we have
\begin{equation}
0 = \varphi(\mathbf{x} + a \mathbf{n}) \approx \varphi(\mathbf{x}) + |\nabla \varphi|a + \frac{1}{2}(\mathbf{n}^{\mathsf T}\mathbf H(\varphi)\mathbf{n})a^2
\end{equation}
where the Hessian matrix $\mathbf H$ is defined as
\begin{equation}
\mathbf H(\varphi) =   \begin{pmatrix}
\frac{\partial ^2 \varphi}{\partial x^2} &&
\frac{\partial ^2 \varphi}{\partial x \partial y} \\
 \frac{\partial ^2 \varphi}{\partial y \partial x} &&
\frac{\partial ^2 \varphi}{\partial y^2}\end{pmatrix}.
\end{equation}
The derivatives are approximated using the central finite difference schemes:
\begin{gather} \label{d1}
\frac{\partial \varphi}{\partial x} \approx \frac{\varphi_{i+1,j}-\varphi_{i-1,j}}{2h},
\\ \label{d2}
\frac{\partial \varphi}{\partial y} \approx
\frac{\varphi_{i,j+1}-\varphi_{i,j-1}}{2h},
\\ \label{d3}
\frac{\partial^2 \varphi}{\partial x^2} \approx
\frac{\varphi_{i+1,j}-2\varphi_{i,j}+\varphi_{i-1,j}}{h^2},
\\ \label{d4}
\frac{\partial^2 \varphi}{\partial y^2} \approx
\frac{\varphi_{i,j+1}-2\varphi_{i,j}+\varphi_{i,j-1}}{h^2},
\\ \label{d5}
\frac{\partial^2 \varphi}{\partial x \partial y} =
\frac{\partial^2 \varphi}{\partial y \partial x} \approx
\frac{\varphi_{i+1,j+1}+\varphi_{i-1,j-1}-\varphi_{i+1,j-1}-
\varphi_{i-1,j+1}}{4h^2}.
\end{gather}

At each orthogonal projection $\mathbf{x}^*$, the geometrical information
needed includes the curvature $\kappa$ and the angle $\theta$ between the
normal direction and the $x$-axis at $\mathbf{x}^*$.
The curvature is defined as usual
\begin{equation}
\kappa =  -\frac{\frac{\partial^2 \varphi}{\partial x^2}(\frac{\partial \varphi}{\partial y})^2 - 2\frac{\partial^2 \varphi}{\partial x \partial y}\frac{\partial \varphi}{\partial x}\frac{\partial \varphi}{\partial y} + \frac{\partial^2 \varphi}{\partial y^2}(\frac{\partial \varphi}{\partial x})^2}
{\left( (\frac{\partial \varphi}{\partial x})^2 + ( \frac{\partial \varphi}{\partial y})^2 \right) ^{\frac{3}{2}}}.
\end{equation}
The unit normal direction is determined by
\begin{equation}
\cos\theta = \frac{\frac{\partial \varphi}{\partial x}}{\sqrt{(\frac{\partial
\varphi}{\partial x})^2+(\frac{\partial \varphi}{\partial y})^2}}, \quad
\sin\theta = \frac{\frac{\partial \varphi}{\partial y}}{\sqrt{(\frac{\partial
\varphi}{\partial x})^2+(\frac{\partial \varphi}{\partial y})^2}}.
\end{equation}
The derivatives at $\mathbf{x}^*$ can be approximated using the bilinear interpolation formula from the derivatives already calculated at
grid points using (\ref{d1})-(\ref{d5}) as following. Given an interface point $\mathbf{x}^*$ and a grid
function $f(\cdot)$, we label the four grid points surrounding $\mathbf{x}^*$
 as
\begin{equation}
(x_i,y_j), (x_{i+1},y_j), (x_{i},y_{j+1}), (x_{i+1},y_{j+1}).
\end{equation}
See Fig. \ref{f1} for an example.
If we define
\begin{equation}
a = \frac{2(x^*-x_i)}{h} - 1, \quad
b = \frac{2(y^*-y_j)}{h} - 1,
\end{equation}
then we can use the  bilinear interpolation scheme
\begin{equation} \label{bilinear}
\begin{aligned}
& \big\{(1-a)(1-b)f_{i,j} + (1+a)(1-b)f_{i+1,j}\\
&+ (1-a)(1+b)f_{i,j+1} + (1+a)(1+b)f_{i+1,j+1}\big\}/4
\end{aligned}
\end{equation}
to approximate $f(\mathbf{x}^*)$. We can calculate $\kappa$,  $\cos \theta$
and $\sin \theta$
at the orthogonal projection point $\mathbf{x}^*$ if we set $f=\kappa(\mathbf{x})$,
 $f=\cos \theta$, and $f=\sin \theta$ respectively using the interpolated
formula above. Note that these quantities are poorly conditioned if $|\nabla \varphi| \sim 0$.
However our choice of $\varphi$ as the signed distance function avoids this difficulty.

\subsection{The immersed interface method for elliptic interface problems}

Given an intermediate approximation of $\mathbf{q}$, the solution to the
Stokes equations can be obtained by solving the three Poisson equations
(\ref{p1})-(\ref{p3}) with
known jump conditions in the solution and the flux:
\begin{equation} \label{poisson}
\begin{gathered}
 (\Psi_x)_x + (\Psi_y)_y  = R(x,y) \\
 [\Psi] = w_1, \quad [\Psi_n]=w_2.
\end{gathered}
\end{equation}
using the immersed interface method \cite{rjl-li:sinum,li:fast,li-ito}.
The finite difference scheme using the immersed interface method
can be simply written as
\begin{equation}
\frac{\Psi_{i-1,j} + \Psi_{i+1,j} + \Psi_{i,j-1} + \Psi_{i,j+1} -
4 \Psi_{ij}}{h^2} = R_{ij} + C_{ij}
\end{equation}
where $C_{ij}=0$ at regular grid points. At an irregular grid point,
$C_{ij}$ depends on the interface curvature, first and second order surface
derivatives of $w_1$ and $w_2$ (see \cite{li:thesis} for its formulation and derivation). The immersed interface method
not only gives a second order accurate solution, but also enables us to
use widely available fast Poisson solvers to solve the discrete systems of
equations. The fast Poisson solver, for example the one from Fishpack
\cite{fish},
only requires $O(N \log(N))$ operations, where $N$ is the total number of
grid points. Therefore the Stokes equations with given ${\bf q}$ can be
solved in $O(3N \log(N))$ operations, which corresponds to the cost of
one GMRES iteration.

\subsection{Computing an interface quantity and its tangential derivatives
using the least squares interpolation scheme}

To determine the correction term $C_{ij}$ we need to know the jumps
$[\Psi] = w_1(s)$, $[\frac{\partial \Psi}{\partial \mathbf{n}}] = w_2(s)$
and the tangential derivatives $w_1'$, $w_1''$ and $w_2'$,
where $s$ is the arc-length parameterization of the interface \cite{li:thesis}.
The problem may be formulated thus: given the values of an interface
function $\phi(s)$ at a set of interface points ${s=s_1,\dots,s=s_n}$,
find the values  of $\phi(s^*)$, $\phi'(s^*)$ and $\phi''(s^*)$.
It is critical to choose an interpolation scheme for the immersed
interface-level set method to work well.

We assume that the values of the interface function $\phi(s)$ are known at
all the projection points from the inner irregular grid points and we want
to interpolate $\phi(s)$ at the projection points from the outer irregular
grid points, and $\phi'(s)$ and $\phi''(s)$ at the projection points
from both inner and outer irregular grid points. See Fig. \ref{f2} for an illustration.

\begin{figure}[hpbt]
\begin{center}
\includegraphics[width=0.5\textwidth]{figures/x2}
\end{center}
\caption{On the interface $\Gamma$, the values of the interface function $\phi(s)$
at the projection points from the inner irregular grid points (dark circles) can be used
to interpolate both the values of $\phi(s)$ at the projection points from the outer irregular
grid points (white circles) and the values of $\phi'(s)$ and
$\phi''(s)$ at all the projection points.}
\label{f2}
\end{figure}


The least squares interpolation scheme for approximating $\phi(s^*)$
(or $\phi'(s^*)$ and $\phi''(s^*)$) can be written as
\begin{equation} \label{se3}
\sum_{i} \gamma_{i}\phi(s_i),
\end{equation}
where $\gamma_{i}$'s are the coefficients that need to be determined from
the Taylor expansion at the interpolation point $s^*$
\begin{equation} \label{se2}
\phi(s_i) = \phi(s^*) + \phi'(s^*)(s_i-s^*) +
\frac{1}{2}\phi''(s^*)(s_i-s^*)^2 + \dots.
\end{equation}

We define
\begin{equation} \label{least1}
a_1 =  \sum_{i}\gamma_{i}, \quad
a_2 =   \sum_{i}(s_i-s^*)\gamma_{i} \quad
a_3 =  \sum_{i}\frac{1}{2}(s_i-s^*)^2\gamma_{i},
\end{equation}
where $(s_i-s^*)$ is the signed arc length from the interface point
at $s = s_i$ to the point at $s = s^*$(see \cite{li:void} for its formula
using the Hermite cubic spline approximation method). Substituting (\ref{se2})
into (\ref{se3}), we obtain the coefficients $\gamma_{i}$ from setting
\begin{equation} \label{se1}
a_1 = 1, \quad  a_2 = 0,  \quad a_3 = 0.
\end{equation}
If we choose more than three interface points for the interpolation, then
the linear system of equations (\ref{se1}) is under-determined, that is, there
are an infinite number of solutions.  We use the singular value
decomposition (SVD) method to solve the system of equations (\ref{se1})
to get the unique solution that has the least 2-norm among all feasible
solutions.

For a given $s = s^*$, the interpolation points are chosen as a few
orthogonal projections of inner grid points ($\varphi_{i,j}< 0$) that are
the first few closest to the point $(X(s^*),Y(s^*))$. The criterion is
the arc length distance such that $|s-s^*| \leq \alpha$, where
$\alpha$ is a parameter. We choose $\alpha=6.5h$.

\subsection{Augmented immersed interface method for the Stokes equations}

The augmented unknown $\mathbf q$ needs to be solved only at orthogonal
projection points from the inner irregular
grid points because we can use interface interpolation to get the values
at outer projection points.

Given a discrete approximation of $(q_1,q_2)$ at $\{\mathbf{X}_k\}$,
we can solve the first three equations \eqref{Sab}-\eqref{Scb} to get an
approximate
solution:  the pressure $\mathbf{P}({\bf Q})$, the scaled velocity
$\tilde{\mathbf{U}}({\bf Q})$
and $\tilde{\mathbf{V}}({\bf Q}).$ Generally the computed velocity
$(\tilde{\mathbf{U}}({\bf Q}),\,\tilde{\mathbf{V}}({\bf Q}))$ does
not satisfy the two augmented equations in \eqref{Scd}, that is,
$(\mathbf{U},\mathbf{V})=(\tilde{\mathbf{U}}/\mu,\tilde{\mathbf{V}}/\mu)$ may not
be continuous across the interface.

Let us concatenate the discrete solutions $\{P_{ij}\}$, $\{U_{ij}\}$, and
$\{V_{ij}\}$ as a vector $\mathbf{\tilde{\mathcal{U}}} $ whose dimension is
$O(3MN)$, where $M$ and $N$ are the number of grid lines in the $x$ and $y$
direction respectively. We denote also the vector of the discrete values of
$(q_1,q_2)$ at the control points  $\{\mathbf{X}_k\}$ (consisted of projection points
only from the inner irregular grid points) by
$\mathbf{Q}$ whose dimension is $O(2 N_b)$. Then the discrete
solution of \eqref{Sab}-\eqref{Scb} given ${\bf Q}$ can be written as
\begin{equation} \label{bigeq1}
A \, \mathbf{\tilde{\mathcal{U}}}  + B \mathbf{Q} = \mathbf{F}_1
\end{equation}
for some vector $\mathbf{F}_1$ and sparse matrices $A$ and $B$. It
requires solving three Poisson equations with different source
terms and jump conditions to get $\mathbf{\tilde{\mathcal{U}}} $.

Once we know the solution $\mathbf{\tilde{\mathcal{U}}} $ given $\mathbf{Q}$, we can use
$(\tilde{\mathbf{U}},\,\tilde{\mathbf{V}})$ and the jump conditions
$[{\tilde{\mathbf{U}}}_{\mathbf{n}}]$ and
$[{ \tilde{\mathbf{V}}}_{\mathbf{n}}]$ which also depend on $\mathbf{Q}$,
 to get
$[\mathbf{U}(\mathbf{Q})]=[\tilde{\mathbf{U}}(\mathbf{Q})/\mu]$ and
$[\mathbf{V}(\mathbf{Q})]=[\tilde{\mathbf{V}}(\mathbf{Q})/\mu]$
at those control points $\{\mathbf{X}_k\}$, $1\le k \le N_b$. If both
$\|\, [\mathbf{U}(\mathbf{Q})]\, \|$ and
$\|\,[\mathbf{V}(\mathbf{Q})]\,\|$
are smaller than a given tolerance, then the method has already
converged and the set of $\mathbf{Q}$, $\tilde{\mathbf{U}}/\mu$,
$\tilde{\mathbf{V}}/\mu$ is an approximate
solution. The interpolation scheme to get
$[\tilde{\mathbf{U}}(\mathbf{Q})/\mu]$ and
$[\tilde{\mathbf{V}}(\mathbf{Q})/\mu]$, which will be explained in
the next sub-section, depends on
$\mathbf{\tilde{\mathcal{U}}} ,\mathbf{Q}$  \emph{linearly}.  Therefore we can write
\begin{equation} \label{stinterp1}
\left( \big[{\tilde{\mathbf{U}}}/{\mu} \big], \,
\big[{\tilde{\mathbf{V}}}/{\mu}\big] \right)^T =
S \,\mathbf{\tilde{\mathcal{U}}}  + E \mathbf{Q} - \mathbf{F}_2,
\end{equation}
where ${S}$ and $E$ are two sparse matrices, and $\mathbf{F}_2$ is a vector.
The matrices depend on the interpolation scheme but do not need to be
actually constructed in the
algorithm. We need to choose such a vector $\mathbf{Q}$ that the continuity
condition for the velocity is satisfied along the interface $\Gamma$.
If we put the two
matrix-vector equations \eqref{bigeq1} and \eqref{stinterp1} together we get
\begin{equation}\label{MatVec1}
\begin{bmatrix}
  A & B \\
  S & E \end{bmatrix}
\begin{bmatrix}
  \mathbf{\mathbf{\tilde{\mathcal{U}}} } \\
 \mathbf{Q}  \end{bmatrix}
=\begin{bmatrix}
  \mathbf{F}_1 \\
  \mathbf{F}_2 \end{bmatrix}.
\end{equation}
Note that $ \mathbf{Q}$ is defined only on a set of points $\{\mathbf{X}_k\}$
on the interface while $\mathbf{\mathbf{\tilde{\mathcal{U}}} }$ is defined at all grid points.
The Schur complement for $\mathbf{Q}$ is
\begin{equation} \label{smalleq}
(E-SA^{-1} B) \mathbf{Q} = \mathbf{F}_2 -SA^{-1} \mathbf{F}_1 =
\bar{\bf F}.
\end{equation}
If we can solve the system above to get $\mathbf{Q}$, then we can get \
$\mathbf{\mathbf{\tilde{\mathcal{U}}} }$ easily. Because the dimension of $\mathbf{Q}$ is much smaller
than that of $\mathbf{\mathbf{\tilde{\mathcal{U}}} }$, we  expect to get a reasonably fast algorithm for
the two-phase Stokes equations if we can solve \eqref{smalleq}
efficiently.

In implementation, the GMRES iterative method \cite{saad} is used to solve
\eqref{smalleq}. The GMRES method only requires matrix-vector
multiplication. We explain below how to evaluate the right hand side
$\bar{\bf F}$ of the Schur complement, and how to evaluate the matrix-vector
multiplication needed by the GMRES iteration. We can see
why we do not need to form the coefficient matrix $E-SA^{-1} B$ explicitly.

\subsubsection{Evaluation of the right hand side of the Schur complement}

First we set $\mathbf{Q}={\bf 0}$ and solve the de-coupled system
\eqref{Sab}-\eqref{Scb},
or \eqref{bigeq1} in the discrete form, to get $\mathbf{\mathbf{\tilde{\mathcal{U}}} \!({\bf 0})}$ which
is $A^{-1}\mathbf{F}_1$ from \eqref{bigeq1}. From the interpolation
\eqref{stinterp1}, we also have
\[
\left( \big[  \mathbf{U}(\mathbf{0}) \big] ,
\big[  \mathbf{V}(\mathbf{0}) \big] \right)^T=
S \,\mathbf{\tilde{\mathcal{U}}} \!(\mathbf{0}) + E \,\mathbf{0} - {\bf F}_2=S \,\mathbf{\tilde{\mathcal{U}}} (\mathbf{0})
 - {\bf F}_2.
\]
Note that the residual
of the Schur complement for $\mathbf{Q}={\bf 0}$ is
\begin{equation} \label{resid2b}
\begin{aligned}
R(\mathbf{0}) & =  (E-SA^{-1}B) \, \mathbf{0}  - \bar{\mathbf{F}}
= - \bar{\mathbf{F}} \\
& =  - \left( {\bf F}_2  -   SA^{-1} {\bf F}_1 \right)
= - {\bf F}_2 + S \, \mathbf{\mathbf{\tilde{\mathcal{U}}} (0)} \\
&= \left( \left [  \mathbf{U}(\mathbf{0}) \right ] , \left [  \mathbf{V}(\mathbf{0}) \right ] \right)^T,
\end{aligned}
\end{equation}
which gives the right hand side of the Schur complement system with an
opposite sign.

\subsubsection{Evaluation of the matrix-vector multiplication}
The matrix-vector multiplication of the Schur complement system given
$\mathbf{Q}$ is obtained from the following two steps:
\begin{description}

\item[Step 1] Solve the coupled system \eqref{Sab}-\eqref{Scb}, or \eqref{bigeq1}
in the discrete form,  to get $\mathbf{\mathbf{\tilde{\mathcal{U}}} \!(Q)}$.

\item[Step 2] Interpolate $\mathbf{\mathbf{\tilde{\mathcal{U}}} \!(Q)}$ using \eqref{stinterp1} to get
$\left( \left [  \mathbf{U}(\mathbf{0}) \right ] , \left [  \mathbf{V}(\mathbf{0}) \right ] \right)^T$.
Then the matrix vector multiplication is
\begin{equation}
\left( E-S A^{-1} B\right)  \mathbf{Q}=\left( \left [  \mathbf{U}(\mathbf{Q}) \right ] , \left [  \mathbf{V}(\mathbf{Q}) \right ] \right)^T-
\left( \left [  \mathbf{U}(\mathbf{0}) \right ] , \left [  \mathbf{V}(\mathbf{0}) \right ] \right)^T.
\end{equation}
This is because
\begin{align*}
\left( E-S A^{-1}B\right)  \mathbf{Q}
& = E \mathbf{Q}-{S}A^{-1}B \mathbf{Q}\\
& = E \mathbf{Q} - S \left(  A^{-1}{\bf F}_1 - \mathbf{\mathbf{\tilde{\mathcal{U}}} \!(Q)} \right)
\quad  \text{(from \eqref{bigeq1}) }, \\
& = E \mathbf{Q} + S \,\mathbf{\mathbf{\tilde{\mathcal{U}}} \!(Q)} - {\bf F}_2 +
{\bf F}_2 - S  A^{-1}{\bf F}_1  \\
& = \left( \big[  \mathbf{U}(\mathbf{Q}) \big] ,
\big[  \mathbf{V}(\mathbf{Q}) \big] \right)^T +
\bar{\mathbf{F}} \quad \mbox{(from \eqref{stinterp1})}, \\
& = \left( \big[  \mathbf{U}(\mathbf{Q}) \big] ,
 \big[  \mathbf{V}(\mathbf{Q}) \big ] \right)^T-
\left( \big[  \mathbf{U}(\mathbf{0}) \big] ,
\big[  \mathbf{V}(\mathbf{0}) \big] \right)^T,
 \quad  \mbox{(from \eqref{resid2b})}.
\end{align*}
\end{description}
Now we can see that a matrix-vector multiplication is equivalent to solving
the coupled system \eqref{Sab}-\eqref{Scb}, or \eqref{bigeq1} in the
discrete form, to get
$\mathbf{\mathbf{\tilde{\mathcal{U}}} }\!,$ and using an interpolation scheme \eqref{stinterp1} to get
$\left( [  \mathbf{U}] , [  \mathbf{V}] \right)^T$ at the control points.

Since we know the right hand side of the linear system of equations and the
matrix-vector multiplication of the coefficient matrix, it is straightforward
to use GMRES or other iterative methods.

In summary, the system of equations for ${\bf Q}$ can be written as
\begin{equation}
G \mathbf Q = \mathbf b,
\end{equation}
where $\mathbf b = -\begin{pmatrix} [\mathbf U(\mathbf 0)] \\
[\mathbf V(\mathbf 0)] \end{pmatrix}$. The coefficient matrix $G$ does not need
to be formed explicitly.
Given a guess of $\mathbf Q^{k}$, the matrix vector multiplication that is
needed for the GMRES iterative method is simply
\begin{equation} \label{mv}
G\mathbf{Q^{k}} = \begin{pmatrix} [\mathbf U(\mathbf Q^{k} )]
\\ [\mathbf V(\mathbf Q^{k})] \end{pmatrix}-\begin{pmatrix}
[\mathbf U(\mathbf 0)] \\ [\mathbf V(\mathbf 0)] \end{pmatrix}.
\end{equation}

\subsection{Computing the residual vector using the least squares
interpolation}
The matrix-vector multiplication (\ref{mv}) only involves the velocity
jumps at the inner projection points. However, the velocity solved from the
Poisson equations (\ref{p1}, \ref{p2}, \ref{p3}) is defined only at  grid points. We use a least squares
interpolation scheme again to get the velocity at the interface points. The
interpolation scheme is critical for the accuracy
and convergence speed of the  GMRES method and thus is very important for
the new method.

For an inner projection point $\mathbf{x}^* = (x^*,y^*)$ and its neighboring
grid points $\{\mathbf{x} = (x_m,y_n)|m,n \in \mathbb N\}$, we can
find
\begin{equation}
(\tilde{u}(\mathbf{x}^*))^+ =
\lim_{\mathbf{x} \to (\mathbf{x}^*)^+}{\tilde{u}(\mathbf{x})}
\end{equation}
using the following interpolation
\begin{equation} \label{se4}
\sum_{m,n} \beta_{m,n}\tilde{u}(x_m,y_n) + C
\end{equation}
where the coefficients $\beta_{m,n}$ and the correction term $C$ are
to be determined. It is more convenient to use  the local coordinate
system in the tangential and normal directions (Fig. \ref{f1}) defined as
\begin{equation}
\begin{gathered}
\xi_m = (x_m - x^*)\cos\theta + (y_n - y^*)\sin\theta, \\
\eta_n = -(x_m - x^*)\sin\theta + (y_n - y^*)\cos\theta.
\end{gathered}
\end{equation}
The Taylor expansion of $\tilde{u}(x_m,y_n)$ at $\left( \xi,\eta\right) = (0,0)$ is
\begin{equation} \label{taylor2}
\begin{aligned}
&\tilde{u}(\xi_m,\eta_n) \\
&=  \tilde{u}^{\pm} +
\frac{\partial \tilde{u}^{\pm}}{\partial \xi}\, \xi_m +
\frac{\partial \tilde{u}^{\pm}}{\partial \eta} \,\eta_n +
\frac{1}{2} \, \frac{\partial ^2 \tilde{u}^{\pm}}{\partial \xi^2}\,\xi_m^2
 + \frac{1}{2}\,
\frac{\partial ^2 \tilde{u}^{\pm}}{\partial \eta^2}\,\eta_n^2 +
\frac{\partial ^2 \tilde{u}^{\pm}}{\partial \xi \partial \eta}\,\xi_m \eta_n
+ \dots.
\end{aligned}
\end{equation}
The ${\pm}$  sign is determined according to whether
$\mathbf{x} \in \Omega ^{\pm}$.
Let us define
\begin{equation} \label{aks}
\begin{gathered}
 a_1^{+} = \sum_{(x_m,y_n) \in \Omega^{+}} \beta_{m,n}, \quad
 a_1^{-} = \sum_{(x_m,y_n) \in \Omega^{-}} \beta_{m,n}, \\
 a_2^{+} = \sum_{(x_m,y_n) \in \Omega^{+}} \beta_{m,n}\xi_m, \quad
 a_2^{-} = \sum_{(x_m,y_n) \in \Omega^{-}} \beta_{m,n} \xi_m, \\
 a_3^{+} = \sum_{(x_m,y_n) \in \Omega^{+}} \beta_{m,n}\eta_n, \quad
 a_3^{-} = \sum_{(x_m,y_n) \in \Omega^{-}} \beta_{m,n} \eta_n, \\
 a_4^{+} = \sum_{(x_m,y_n) \in \Omega^{+}} \beta_{m,n}\frac{1}{2}\xi_m^2, \quad
 a_4^{-} = \sum_{(x_m,y_n) \in \Omega^{-}} \beta_{m,n} \frac{1}{2}\xi_m^2, \\
 a_5^{+} = \sum_{(x_m,y_n) \in \Omega^{+}} \beta_{m,n}\frac{1}{2}\eta_n^2, \quad
 a_5^{-} = \sum_{(x_m,y_n) \in \Omega^{-}} \beta_{m,n} \frac{1}{2}\eta_n^2, \\
 a_6^{+} = \sum_{(x_m,y_n) \in \Omega^{+}} \beta_{m,n}\xi_m\eta_n, \quad
 a_6^{-} = \sum_{(x_m,y_n) \in \Omega^{-}} \beta_{m,n} \xi_m\eta_n.
\end{gathered}
\end{equation}
Substituting (\ref{taylor2}) into (\ref{se4}), the coefficients $\beta_{m,n}$ of the interpolation scheme for $\tilde{u}^{+}$ are seen to be the
solution of the
following linear system of equations
\begin{equation} \label{lineq2}
\begin{gathered}
  a_1^{+}+a_1^{-} = 1, \quad a_2^{+}+a_2^{-} = 0,
\quad a_3^{+}+a_3^{-}  = 0, \\
 a_4^{+}+a_4^{-} = 0, \quad a_5^{+}+a_5^{-} = 0, \quad a_6^{+}+a_6^{-} = 0.
\end{gathered}
\end{equation}
The correction term is
\begin{equation} \label{correc2}
C =  a_1^{-}[\tilde{u}] + a_2^{-}\big[\frac{\partial
 \tilde{u}}{\partial \xi} \big] + a_3^{-}
 \big[\frac{\partial \tilde{u}}{\partial \eta}\big]
 + a_4^{-}\big[\frac{\partial^2 \tilde{u}}{\partial \xi^2} \big]
+ a_5^{-} \big[\frac{\partial^2 \tilde{u}}
 {\partial \eta^2} \big]  + a_6^{-} \big[\frac{\partial^2 \tilde{u}}
 {\partial \xi \partial \eta} \big].
\end{equation}
Again, we choose more than six grid points for the interpolation to get an
under-determined system of equations. The linear system of equations
is solved by the SVD method. The jump conditions in the correction term
can be represented by the jump in the solution and the flux (see \cite{li-ito-lai}).

Similarly, we can get an interpolation scheme for $\tilde{u}^{-}$. Once we
have $\tilde{u}^{\pm}$ at all orthogonal projections of inner irregular
grid points, it is also straightforward to get
$\begin{pmatrix} [\mathbf U] \\ [\mathbf V] \end{pmatrix}$
which is needed in the matrix-vector multiplication.

To get the interpolation grid points for a given inner projection point
$\mathbf{x}^* = (x^*,y^*)$, we choose at least 6 closest grid points to
$(x^*,y^*)$ as the interpolation stencil so that we can have a third order accurate
interpolation scheme. To avoid excess points we use $|\left( \xi_m,\eta_n\right)| \leq 2.5h$.

\subsection{The level set method for moving interfaces}
We use the level set method first introduced by
Osher \& Sethian \cite{osher-sethian} to solve the moving interface problems
modeled by the incompressible Stokes equations. The level set function
is evolved according to the level set function (a Hamilton-Jacobi equation)
\begin{equation} \label{level1}
\frac{\partial \varphi}{\partial t} + V_n |\nabla \varphi|= 0
\end{equation}
where $V_n = \mathbf{u} \cdot \mathbf{n}$ is the component of the velocity in the
level sets normal direction and $\varphi = \varphi(\mathbf{x},t)$. The level set function $\varphi$ is often chosen as the
signed distance function which satisfies the Eikonal equation
\begin{equation}
|\nabla \varphi| = 1
\end{equation}
at least in a neighborhood of the interface. To ensure that the level set
function $\varphi$  remains the signed distance
function, we solve the following re-initialization equation at each time step after interface advection
until $|\nabla \varphi| = 1 + \mathrm{O}(h^2)$:
\begin{equation} \label{level2}
\begin{gathered}
\varphi(\mathbf{x},t=0) = \varphi_0, \\
\frac{\partial \varphi}{\partial t} =\text{sign}(\varphi_0)(1-|\nabla \varphi|).
\end{gathered}
\end{equation}

Numerical methods \cite{osher:book, sethian:book} for Hamilton-Jacobi (HJ)
equations can be used to solve the level set equations
(\ref{level1})-(\ref{level2}).
In our implementation, we use the explicit forward Euler scheme for the time
derivative and the third order WENO (Weighted Essentially Non-Oscillatory)
scheme for the spatial derivatives.


The standard level set methods do not preserve volume well. Sussman et. al. \cite{sussman:re}
incorporated the volume conservation condition into the re-initialization
equation (\ref{level2}) and a hybrid particle
level set approach has also been proposed \cite{ron-hybrid}. We use a simple
approach described in \cite{g1} to preserve the area.

\section{Numerical examples} \label{sec:6}

In this section we  present one example with fixed interface, one with a
moving interface and one with four moving interfaces.
All the computations are done on a P4 2.80GHz DELL desktop PC with
512MB memory. The computational
domain is $\Omega = [-2,2] \times [-2,2]$. We choose to fix
$\mu^+ = 1.0$, so $\lambda =  \frac{\mu^-}{\mu^+} = \mu^-$.

\subsection{A general fixed interface test problem}

We first present some results to a general fixed interface test problem.
The interface is chosen to be a unit circle
\begin{equation}
\varphi(x,y) = \sqrt{x^2+y^2} - 1.0.
\end{equation}
The initial guess of GMRES for the augmented variable $\mathbf Q$ takes zero values.
The non-physical exact test solutions of $(\mathbf{u},p)$ are constructed as \cite{li-ito-lai}
\begin{equation}
\begin{gathered}
u = \begin{cases}
 \frac{y}{4}& \text{if $x^2+y^2 < 1$}, \\
 \frac{y}{4}(x^2+y^2)& \text{if $x^2+y^2 \ge 1$};
\end{cases}
\\
v = \begin{cases}
 -\frac{x}{4}(1-x^2)& \text{if $x^2+y^2 < 1$}, \\
 -\frac{xy^2}{4}& \text{if $x^2+y^2 \ge 1$}; \\
\end{cases}
\\
p = \begin {cases}
 (-\frac{3}{4}x^3+\frac{3}{8}x)y& \text{if $x^2+y^2 < 1$}, \\
0& \text{if $x^2+y^2 \ge 1$}. \\
\end {cases}
\end{gathered}
\end{equation}
The body force term $\mathbf g = (g_1,g_2)^T$ is chosen to be
\begin{equation}
\begin{gathered}
g_1 = \begin {cases}
 \left( -\frac{9}{4}x^2+\frac{3}{8} \right) y& \text{if $x^2+y^2 < 1$}, \\
-2\mu^+y& \text{if $x^2+y^2 \ge 1$};
\end{cases}
\\
g_2 = \begin{cases}
 -\frac{3}{4}x^3+\frac{3}{8}x-\frac{3\mu^-}{2}x& \text{if $x^2+y^2 < 1$}, \\
 \frac{\mu^+}{2}x& \text{if $x^2+y^2 \ge 1$}.
\end {cases}
\end{gathered}
\end{equation}
And the singular interface force terms are

\begin{equation}
\begin{gathered}
\hat{f}_1 = (\frac{3}{4}\cos^3\theta-\frac{3}{8}\cos\theta)\sin\theta-\frac{3}{2}\cos^3\theta\sin\theta[\mu],\\
\hat{f}_2 = \frac{1}{2}\mu^++\frac{3}{4}[\mu]\cos^2\theta(1-2\cos^2\theta).
\end{gathered}
\end{equation}
To accommodate the exact solutions, the jump condition for
$\frac{\partial p}{\partial n}$ is modified by
adding a correction term,
which does not change the difficulty of the problem.
Note that this is a very general fixed interface problem with non-homogeneous jumps and
tangential derivatives in all quantities.

\begin{figure}[hptb]

\begin{center}
\begin{tabular}{cc}
(a) & (b)\\
\includegraphics[width=0.47\textwidth]{figures/fig2}
& \includegraphics[width=0.47\textwidth]{figures/fig1}\\
(c) & (d)\\
 \includegraphics[width=0.47\textwidth]{figures/fig3}
& \includegraphics[width=0.47\textwidth]{figures/fig4}
\end{tabular}
\end{center}
\caption{Convergence analysis in the $\mathcal{L}^{\infty}$ norm for four cases in
$\log$-$\log$ scale. The slope of the linear regression line is the convergence rate which
is close to number $2.0$ for all cases. (a) $\mu^- = 0.001$ and $\mu^+ = 1.0$.
(b) $\mu^- = 0.1$ and $\mu^+ = 1.0$. (c) $\mu^- = 10$ and $\mu^+ = 1.0$.
(d) $\mu^- = 1000$ and $\mu^+ = 1.0$. Linear regressions $R^2>0.98$ for all cases.}
\label{fig1}
\end{figure}

\begin{figure}[hpbt]
\begin{center}
\includegraphics[width=0.8\textwidth]{figures/fig5}
\end{center}
\caption{Algorithm efficiency analysis for 32 by 32 through 896 by 896 grids.
 Left axis represents
number of GMRES iterations and right axis represents computational
CPU cost with unit second using log scale. GMRES iterations (thick curves) remain relatively constant with
grid size, and CPU time (thin curves) $\sim$ grid number$^2$.}
\label{fig2}
\end{figure}

We use linear regression analysis to find the convergence order of all
quantities. For this purpose, we choose $m = n = 32 + 16k(k = 0, \dots, 54)$
and run the convergence analysis for
$\lambda = 0.001, 0.1, 10$ and
$1000$.  Fig. \ref{fig1} is the  grid refinement analysis of the test problem.
The error is the $\mathcal{L}^{\infty}$ error at all grid points. The slope of the linear
regression line is the order of the convergence.
It is clear from these results that the algorithm is
second order accurate in velocity $\mathbf{u} = (u,v)^T$, as expected, and
almost second order accurate in pressure $p$.
For this problem, it also seems that with larger $\lambda$, the numerical solution for
$u$ (or $v$) tends to be more accurate than the solution for $p$.

Fig. \ref{fig2} shows the algorithm efficiency analysis.
The number of GMRES iterations seems to
be independent of grid size for fixed $\lambda$, and increases slowly as
$\lambda$ increases. Most simulations finish within two minutes.

\begin{figure}[hptb]
\begin{center}
\includegraphics[width=0.6\textwidth]{figures/05}
\end{center}
\caption{A moving interface simulation with $\lambda = 0.5$,
at time $t = 0, 0.86, 1.90$
and $6.75$. Starting from the ellipse, the
interface evolves to a circle. The inner reference circle (dash curve) helps to show the
symmetry of the equilibrium interface.
The final area inside the interface is $3.14\approx \pi$ so our algorithm preserves area well.}
\label{fig3}
\end{figure}

\begin{figure}[hpbt]
\begin{center}
\includegraphics[width=0.6\textwidth]{figures/fig7}
\end{center}
\caption{Comparison of viscosity effects on interface motion. Three curves (thin curves) at time
 $t = 1.5$, for $\lambda = 0.5,1.0,2.0$. $\lambda = 0.5$ approaches equilibrium
(thick curve) fastest.}
\label{fig4}
\end{figure}

\begin{figure}[hptb]
\begin{center}
\begin{tabular}{cc}
(a) & (b) \\
\includegraphics[width=0.47\textwidth]{figures/x3}
& \includegraphics[width=0.47\textwidth]{figures/x4}
\end{tabular}
\end{center}
\caption{A moving interface simulation starting with four elliptic interfaces (white boundary).
Streamlines are drawn and symmetry of streamlines are observed as expected. We choose $\lambda=1$ here.
(a) $t=5$. (b) $t=50$.}
\label{fig5}
\end{figure}


\subsection{A moving interface test problem}

We present a moving interface test problem where the surface tension is the only interface force, i.e.
the normal force $\hat{f}_1 = \alpha \kappa$ where $\alpha > 0$ and $\kappa$ is the curvature as
defined before; we set the tangential force $\hat{f}_2 = 0.0$. We choose $\alpha = 1.0$ and $m = n = 64$
in our simulations. It is known that starting from an arbitrary interface shape without any body forces,
i.e. $\mathbf{g = 0}$, the interface will evolve to the equilibrium circle shape \cite{rjl-li:stokes}.
The initial interface is chosen to be an ellipse
\begin{equation}
\varphi(x,y) = \sqrt{\frac{(x+y)^2}{2a^2}+\frac{(-x+y)^2}{2b^2}} - 1.0
\end{equation}
where $a = 2.0$ and $b = 0.5$. We use periodic boundary conditions for both $\mathbf{u}$ and $p$.
See Fig. \ref{fig3} for an example of the interface evolution. Also see Fig. \ref{fig4} for a
comparison of the evolving interfaces with different $\lambda$'s at a fixed time. It is clear from
Fig. \ref{fig4} that the more viscous fluid moves slower. For all the moving interface simulations, the
area conservation algorithm \cite{g1} preserves the initial area at each time step exactly.

\subsection{A moving multiple-interface problem}

We present another moving interface test problem with a
multiply-connected interface consisting initially of four ellipses. As
in the other test problems, we use periodic boundary conditions on the
edges of the simulation rectangle. The two fluids have equal viscosity,
the surface tension is given as $\hat{f}_1 = \alpha \kappa$ where $\alpha = 0.1$,
and we set $\hat{f}_2 = 0$. The initial interfaces are easily constructed using
the union of individual level set functions \cite{osher:book, sethian:book}. As the four
symmetrically placed ellipses relax to circles, we observe preservation
of the symmetries in the flow. Equilibrium is reached when the pressure
is uniformly distributed.

\section*{Discussion} \label{sec:7}

In this paper we proposed a new numerical method for Stokes equations
with interfaces, by combining the augmented interface method with the
level set method. Our method handles moving interfaces and multiply
connected interfaces more easily than a spline approach, and could more
easily be extended to 3D. Our algorithm is shown to be 2nd order
accurate in both velocity and pressure, and conserves volume well. Our
method should prove useful for interface problems in fluid flow and
biological models.

One of the remaining open questions is how to use a preconditioning technique
for the Schur complement system since the coefficient matrix for the Schur
complement system is not explicitly formed. Even the simplest diagonal
preconditioning technique needs to use the fast Stokes solver $N_b$ times
which would cost more than that of the entire solution process for the
testing examples in this paper. On the other hand, as discussed in
\cite{li:fast}, a good interpolation
scheme for the residual that takes advantages of the interface/jump conditions
can reduce the number of iteration significantly. Therefore, it is more
practical to tune-in the interpolation scheme as a preconditioner than to
use preconditioning techniques that require structures of the coefficient
matrix. Fortunately, for all the test examples in this paper, the number
of GMRES iterations is modest (fewer than $55$).


\subsection*{Acknowledgments}
This work was partially supported by grants: 0201094 from
NSF/NIH, 43751-MA from ARO, and DMS-0412654 from NSF.
The first author would like to thank Prof. John Bishir
 for helpful discussions and proofreading the paper.

\begin{thebibliography}{00}

\bibitem {chorin:proj} A. J. Chorin,
\emph{Numerical solution of the Navier-Stokes equations},
Math. Comput. 22:745-762, 1968.

\bibitem {cortez-stokes} R. Cortez,
\emph{The Method of Regularized Stokeslets},
sisc, 23: 1204-1225, 2001.

\bibitem {ron-hybrid} D. Enright, R. Fedkiw, J. Ferziger and I. Mitchell,
\emph{A hybrid particle level set method for improved
interface capturing}, J. Comp. Phys., 183:83-116, 2002.

\bibitem {fauci1} L. J. Fauci and A. McDonald,
\emph{Sperm motility in the presence of boundaries}, B. Math. Biol.,
57:679-699, 1994.

\bibitem {fo-pe} A. L. Fogelson and C. S. Peskin;
\emph{Numerical solution of the three dimensional Stokes equations in
	 the presence of suspended particles},
Proc. SIAM Conf. Multi-phase Flow, SIAM. June 1986. 	 

\bibitem {givelberg1} E. Givelberg and J. Bunn,
\emph{A comprehensive three-dimensional model of the
cochlea}, J. Comput. Phys., 191:377-391, 2003.

\bibitem {g1} S. Gro\ss, V. Reichelt and A. Reusken,
\emph{A finite element based level set method for two-phase incompressible flows},
IGPM-Report 243, RWTH Aachen, 2004.

\bibitem {jadhav1} S. Jadhav, C. D. Eggleton and K. Konstantopoulos,
\emph{A 3-D computational model predicts that
cell deformation affects selectin-mediated leukocyte rolling},
 Biophys. J., 88:96-104, 2005.

\bibitem {rjl-li:sinum} R. J. LeVeque and Z. Li,
\emph{The immersed interface method for elliptic equations with discontinuous coefficients
and singular sources}, SIAM J. Numer. Anal., 31:1019, 1994.

\bibitem {rjl-li:stokes} R. J. LeVeque and Z. Li,
\emph{Immersed interface methods for Stokes flow with elastic boundaries or
surface tension}, SIAM J. Sci. Comput., 18:709-735, 1997.

\bibitem {li:thesis} Z. Li,
 \emph{The immersed interface method - a numerical approach for partial differential equations
with interfaces}, PhD thesis, University of Washington, 1994.

\bibitem {li:review} Z. Li,
 \emph{An overview of the immersed interface method and its applications},
Taiwanese J. Mathematics, 7:1-49, 2003.

\bibitem {li:fast} Z. Li,
\emph{A Fast Iterative Algorithm for Elliptic Interface Problems}, SIAM J. Numer. Anal.,
35:230-254, 1998.

\bibitem {li-ito} Z. Li and K. Ito,
\emph{Maximum Principle Preserving Schemes for Interface Problems with
Discontinuous Coefficients}, SIAM J. Sci. Comput., 23:1225-1242, 2001.

\bibitem{li-ito-lai} Z. Li and K. Ito and M-C. Lai;
\emph{An augmented approach for Stokes equations with a
         discontinuous viscosity and singular forces},
NCSU-CRSC Tech. Report: CRSC-TR04-23, North Carolina 
                State University, 2004.

\bibitem {li-lubkin} Z. Li and S. R. Lubkin,
\emph{Numerical Analysis of Interfacial 2D Stokes Flow with
Discontinuous Viscosity and Variable Surface Tension},
Numerical Methods in Fluids, 37:525-540, 2001.

\bibitem {li:void} Z. Li, H. Zhao and H. Gao,
\emph{A numerical study of electro-migration voiding by evolving level
set functions on a fixed Cartesian grid}, J. Comput. Phys., 152:281-304, 1999.

\bibitem {lubkin-li} S. R. Lubkin and Z. Li,
\emph{Force and Deformation on Branching Rudiments: Cleaving Between Hypotheses},
 Biomechanics and Modeling in Mechanobiology, 1(1):5-16, 2002.

\bibitem {mayo-pe} A.  Mayo and C. S. Peskin,
\emph{An implicit numerical method for fluid dynamics problems with
  immersed elastic boundaries}, Contemp. Math., 141:261-277. 1991.

\bibitem {mcqueen1} D. M. McQueen and C. S. Peskin,
\emph{Heart simulation by an immersed boundary method with formal second-order
accuracy and reduced numerical viscosity},
In Mechanics for a New Millennium, Proceedings of the International
Conference on Theoretical and Applied Mechanics
(H. Aref and J. W. Phillips, eds.), Kluwer Academic Publishers, 2001.

\bibitem {osher:book} S. J. Osher and R. P. Fedkiw,
\emph{Level set methods and dynamic implicit surfaces}, Springer-Verlag, 2002.

\bibitem {osher-sethian} S. J. Osher and J. A. Sethian,
\emph{Fronts propagating with curvature-dependent speed: algorithms based on
Hamilton-Jacobi formulations}, J. Comput. Phys., 79:12-49, 1988.

\bibitem {fish} P. Swarztrauber and R. Sweet,
\emph{Efficient FORTRAN subprograms for the solution of elliptic
partial differential equations}, NCAR Technical Note-TN/IA-109, 1975.

\bibitem {peskin1} C. S. Peskin,
\emph{The immersed boundary method}, Acta Numerica, 1-39, 2002.

\bibitem {peskin2} C. S. Peskin and D. M. McQueen,
\emph{Fluid dynamics of the heart and its valves}, In
Case Studies in Methematical Modeling: Ecology, Physiology,
 and Cell Biology (H. G. Othmer, F. R. Adler, M. A. Lewis
and J. C. Dallon, eds.), Prentice-Hall, Englewood Cliffs, NJ, 309-337, 1996.

\bibitem {POZRIKIDIS-eds} C. Pozrikidis eds.,
\emph{Modeling and simulation of capsules and biological cells}, Chapman \& Hall, London, 2003.
\bibitem {saad1} Y. Saad and M. Schultz, \emph{GMRES: a generalized
minimal residual algorithm for solving nonsymmetric linear systems},
SIAM J. Sci. Statist. Comput., 7:856-869, 1986.

\bibitem{saad} Y. Saad,
\emph{GMRES: A Generalized minimal residual algorithm for solving
  nonsymmetric linear systems}, sissc, 7:856-869, 1986.

\bibitem {sethian1} J. A. Sethian and P. Smereka,
\emph{Level set methods for fluid interfaces},
Annu. Rev. Fluid Mech., 35:341-72, 2003.

\bibitem {sethian:book} J.~A.~Sethian,
\emph{Level Set Methods and Fast Marching methods},
Cambridge University Press, nd edition, 1999.

\bibitem {sussman:re} M. Sussman and E. Fatemi,
\emph{An efficient, interface-preserving level set redistancing algorithm
and its application to interfacial incompressible fluid flow},
 SIAM J. Sci. Comput., 20:1165-1191, 1999.

\bibitem{tu-peskin} C. Tu and C. S. Peskin,
\emph{Stability and instability in the computation of flows with moving
  immersed boundaries: a comparison of three methods},
 sissc, 13:1361-1376, 1992.

\bibitem {wang1} X. Wang and W. K. Liu,
\emph{Extended immersed boundary method using FEM and RKPM}, Comput. Methods
Appl. Mech. Eng., 193:1305-1321, 2004.

\bibitem {wiegmann1} A. Wiegmann and K. P. Bube,
\emph{The explicit-jump immersed interface method: finite difference
methods for PDEs with piecewise smooth solutions},
SIAM J. Numer. Anal., 37:827-862, 2000.

\bibitem {zhang1} L. Zhang, A. Gerstenberger, X. Wang and W. K. Liu,
\emph{Immersed finite element method}, Comput.
Methods Appl. Mech. Engrg., 193:2051-2067, 2004.

\end{thebibliography}

\end{document}