\section{Introduction\label{s:intro}}
\input introduction.tex

\section{Formulation\label{s:formulation}} 
\input formulation.tex

\section{Method\label{s:method}} 
\input method.tex

\section{Computing Local Averages of Pressure and Stress\label{s:aver}}
\input aver.tex

\section{Results\label{s:results}} 
\input results.tex

\section{Conclusions\label{s:conclusions}}
\input conclusions.tex


\begin{appendices}
\section{Error estimates for near-singular integration \label{A:AppendixA}} 
\input appen1.tex

\section{Jumps in pressure and stress \label{A:AppendixB}} 
\input appen2.tex

\section{Variable curvature formulation \label{A:AppendixC}}
\input appen3.tex
\end{appendices}


\bibliographystyle{plainnat} 






\subsection{Integral equation formulation}
We give the formulation for the general case in which we have viscosity
contrast between the interior and exterior fluid ($\mu_0$ is the
viscosity of the exterior fluid, $\mu_p$ is the viscosity of the
interior fluid for vesicle $p$, and $\nu_p = \mu_p/\mu_0$), and solid
boundaries with a prescribed velocity.  First we introduce $\SS_{pq}$
and ${\mathcal{D}}_{pq}$, the single- and double-layer potentials for Stokes flow,
where the constant factors are chosen so that our formulation is in
agreement with Pozrikidis~\cite[equation 2.2]{pozrikidis2001b}. The
subscripts denote the potentials induced by a hydrodynamic density on
the membrane of vesicle $q$ and evaluated on the membrane of vesicle
$p$:
\begin{align*}
  \SS_{pq}[{\mathbf{f}}]({\mathbf{x}}) &:= \frac{1}{4\pi\mu_{0}}\int_{\gamma_{q}}\left(
    -\boldsymbol{I} \log \rho
  + \frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}} \right){\mathbf{f}}\, ds_{{\mathbf{y}}},
  && {\mathbf{x}} \in \gamma_{p}, \\
  {\mathcal{D}}_{pq}[{\mathbf{u}}]({\mathbf{x}}) &:= \frac{1-\nu_{q}}{\pi}\int_{\gamma_{q}}
    \frac{{\mathbf{r}} \cdot {\mathbf{n}}}{\rho^{2}}\frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}}{\mathbf{u}}\,
  ds_{{\mathbf{y}}}, && {\mathbf{x}} \in \gamma_{p},
\end{align*}
where ${\mathbf{r}} = {\mathbf{x}} - {\mathbf{y}},$ and $\rho = \|{\mathbf{r}}\|_{2}$.      
Also, we define
\[\SS_{p} := \SS_{pp}\quad\mbox{and}\quad {\mathcal{D}}_{p} := {\mathcal{D}}_{pp},\]
to indicate vesicle self-interactions.  Next, we define
\begin{align*}
  {\mathcal{E}}_{pq}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) &= \SS_{pq}[{\mathbf{f}}]({\mathbf{x}}) + {\mathcal{D}}_{pq}[{\mathbf{u}}]({\mathbf{x}}), && {\mathbf{x}} \in \gamma_{p}, \\
  {\mathcal{E}}_{p}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) &= \sum_{q=1}^{M} {\mathcal{E}}_{pq}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}), && {\mathbf{x}} \in \gamma_{p}.
\end{align*}
${\mathcal{B}}$ is the completed double-layer operator for confined
Stokes flow with density ${\boldsymbol\eta}$
\begin{align*}
  {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}) = {\mathcal{D}}_{\Gamma}[{\boldsymbol\eta}]({\mathbf{x}}) + 
    \sum_{q=1}^{M}R[\xi_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}) + 
    \sum_{q=1}^{M}S[{\boldsymbol\lambda}_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}), 
    \quad {\mathbf{x}} \in \gamma \cup \Gamma.
\end{align*}
If ${\mathbf{x}} \in \Gamma_{0}$, the rank one modification
$\mathcal{N}_{0}[{\boldsymbol\eta}]({\mathbf{x}}) = \int_{\Gamma_{0}}({\mathbf{n}}({\mathbf{x}}) \otimes
{\mathbf{n}}({\mathbf{y}})) {\boldsymbol\eta}({\mathbf{y}})ds_{{\mathbf{y}}}$ is added to ${\mathcal{B}}$ which removes the
one-dimensional null space of the corresponding integral
equation~\cite{pozrikidis1992}.  The Stokeslets and rotlets are defined
as
\begin{align*}
  S[{\boldsymbol\lambda}_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}) = \frac{1}{4\pi\mu_{0}}
    \left(\log \rho + \frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}}\right)
    {\boldsymbol\lambda}_{q}({\boldsymbol\eta}) 
  \qquad \text{and} \qquad
  R[\xi_{q}({\boldsymbol\eta}),{\mathbf{c}}_{q}]({\mathbf{x}}) = \frac{\xi_{q}({\boldsymbol\eta})}{\mu_{0}}
    \frac{{\mathbf{r}}^{\perp}}{\rho^{2}},
\end{align*}
where ${\mathbf{c}}_{q}$ is a point inside $\omega_{q}$, $\omega_{q}$ is the
interior of vesicle $q$, ${\mathbf{r}}={\mathbf{x}} - {\mathbf{c}}_{q}$, and ${\mathbf{r}}^{\perp} =
(r_{2},-r_{1})$.  The operator ${\mathcal{B}}$ satisfies the jump condition
\begin{align*}
  \lim_{\substack{{\mathbf{x}} \rightarrow {\mathbf{x}}_{0} \\ {\mathbf{x}} \in \Omega}}
    {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}) = -\frac{1}{2}{\boldsymbol\eta}({\mathbf{x}}_{0}) + 
    {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}_{0}), \qquad {\mathbf{x}}_{0} \in \Gamma,
\end{align*}
and the size of the Stokeslets and rotlets are
\begin{align*}
  {\boldsymbol\lambda}_{q,i} = \frac{1}{2\pi} \int_{\gamma_{q}} 
    {\boldsymbol\eta}_{i}({\mathbf{y}}) ds_{{\mathbf{y}}}, \quad i=1,2
  \qquad \text{and} \qquad
  \xi_{q} = \frac{1}{2\pi}\int_{\gamma_{q}} {\mathbf{y}}^{\perp}
    {\boldsymbol\eta}({\mathbf{y}})ds_{{\mathbf{y}}}.
\end{align*}
The inextensibility constraint is written in operator form as
\begin{align*}
  \mathcal{P}[{\mathbf{u}}]({\mathbf{x}}) = {\mathbf{x}}_{s} \cdot {\mathbf{u}}_{s}.
\end{align*}

Putting everything together, the integral formulation equation of~\eqref{e:vesicles-pde} is
given by~\cite{rah:vee:bir}
\begin{subequations}
\label{e:ves:dyn}
\begin{align}
  &(1+\nu_{p}) {\mathbf{u}}({\mathbf{x}}) = {\mathcal{E}}_{p}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) + 
  {\mathcal{B}}_{p}[{\boldsymbol\eta}]({\mathbf{x}}), \quad &&{\mathbf{x}} \in \gamma_{p},
  \label{e:vesicle:dynamic} \\
  &\mathcal{P}[{\mathbf{u}}]({\mathbf{x}}) = 0, &&{\mathbf{x}} \in \gamma_{p},
  \label{e:inextensibility} \\
  &\mathbf{U}({\mathbf{x}}) = -\frac{1}{2}{\boldsymbol\eta}({\mathbf{x}}) + {\mathcal{E}}_{\Gamma}[{\mathbf{f}},{\mathbf{u}}]({\mathbf{x}}) + 
    {\mathcal{B}}[{\boldsymbol\eta}]({\mathbf{x}}), &&{\mathbf{x}} \in \Gamma,
  \label{e:BIE}
\end{align}
\end{subequations}
where ${\mathbf{f}}$ is the hydrodynamic density, which in our case is the stress
jump on the vesicle interface, which also is referred to as the traction
jump, and is given by 
\begin{align*}
  {\mathbf{f}} = -\kappa_{b}{\mathbf{x}}_{ssss} + (\sigma {\mathbf{x}}_{s})_{s}.
\end{align*}
Since ${\mathbf{u}} = d{\mathbf{x}}/dt$ and ${\mathbf{f}}$ depends on $\sigma$ and
${\mathbf{x}}$, \eqref{e:ves:dyn} is a system of integro-differential-algebraic equations
for ${\mathbf{x}},\sigma$, and ${\boldsymbol\eta}$. 




\subsection{Spatial Discretization\label{s:spaceDisc}} 
Let ${\mathbf{x}}(\alpha),$ with $\alpha \in (0,2\pi]$, be a parametrization of
the interface $\gamma_p$, $\{\alpha_k = 2k\pi/n\}_{k=1}^n$ be $n$
material points, and
\begin{align*}
  {\mathbf{x}}(\alpha) = \sum_{k = -n/2+1}^{n/2} \hat{{\mathbf{x}}}(k) e^{ik\alpha}.
\end{align*}
Using the FFT to compute $\hat{{\mathbf{x}}}$, we compute derivatives of ${\mathbf{x}}$
with spectral accuracy, assuming that $\gamma$ belongs to $C^{\infty}$.
In particular, we can compute the arclength derivative with spectral
accuracy since 
\begin{align*}
  \pderiv{}{s} = \pderiv{\alpha}{s}\pderiv{}{\alpha} = 
    \frac{1}{\|\partial{\mathbf{x}}/\partial\alpha\|}\pderiv{}{\alpha},
\end{align*}
where $s$ is the arclength parametrization.

Since the single-layer potential $\SS$ has a logarithmic singularity, we
use the hybrid Gauss-trapezoid quadrature rule given in Table 8 of
\cite{alpert1999}.  The error of this quadrature rule is $\mathcal{O}(h^{8}
\log h)$ for integrands with logarithmic singularities. Letting ${\mathbf{x}}_k =
{\mathbf{x}}(\alpha_k)$, the quadrature rule is
\begin{align*}
  \mathcal S[{\mathbf{f}}]({\mathbf{x}}) \approx \sum_{k=1}^{n+m} w_{k}
      S({\mathbf{x}},{\mathbf{x}}_{k}) {\mathbf{f}}({\mathbf{x}}_{k})\| {\mathbf{x}}_{\alpha,k}\|,
\end{align*}
where $n$ is the number of nodes, $m$ is a number of additional
quadrature nodes, $w_{k}$ are the quadrature weights, ${\mathbf{x}}_{k}$ are the
quadrature abscissae, and ${\mathbf{x}}_{\alpha}$ is $\partial {\mathbf{x}}/\partial
\alpha$.  The first $n$ abscissae are the usual equispaced quadrature
points.  Then, the other $m$, which is determined by the desired order
of convergence for the integral, are additional quadrature abscissas
clustered around the singularity.

The double-layer potential has no singularity in two dimensions since
\begin{align*}
  \lim_{\substack{{\mathbf{x}}' \to {\mathbf{x}} \\ {\mathbf{x}}' \in \gamma}} D({\mathbf{x}}',{\mathbf{x}}) = 
    \frac{\kappa({\mathbf{x}})}{2\pi}({\bf t}({\mathbf{x}}) \otimes {\bf t}({\mathbf{x}})),
  \quad {\mathbf{x}} \in \gamma,
\end{align*}
where $\kappa$ is the curvature and ${\bf t}$ is the unit tangent.
Thus, a composite trapezoid rule will give spectral accuracy since the
integrand is periodic and smooth.


\subsection{Time discretization}
Time derivatives are discretized as
\begin{align*}
  \frac{d{\mathbf{x}}}{dt} \approx \frac{\beta{\mathbf{x}}^{N+1} - {\mathbf{x}}^{0}}{\Delta t},
\end{align*}
where ${\mathbf{x}}^{0}$ is a linear combination of previous time steps.
Operators and terms that are treated explicitly are discretized at
${\mathbf{x}}^{e}$ which is also a linear combination of previous time steps.
The simplest IMEX method is IMEX Euler which is given by $\beta=1$,
${\mathbf{x}}^{0} = {\mathbf{x}}^{N}$, and ${\mathbf{x}}^{e} = {\mathbf{x}}^{N}.$ The second-order time
integrator we use is given by $\beta=3/2,$
\begin{align*}
{\mathbf{x}}^{0} = 2{\mathbf{x}}^{N} - \frac{1}{2}{\mathbf{x}}^{N-1}, \hspace{20pt}
{\mathbf{x}}^{e} = 2{\mathbf{x}}^{N} - {\mathbf{x}}^{N-1}.
\end{align*}
To avoid significant time stepping constraints, we treat the bending
term ${\mathbf{f}}_{b}$ and tension term ${\mathbf{f}}_{\sigma}$ semi-implicitly.  An
approximation for the position and tension of vesicle $p$ at time $N+1$
satisfies
\begin{subequations}
  \begin{align}
    \frac{\alpha_{p}}{\Delta t}(\beta {\mathbf{x}}_{p}^{N+1} - {\mathbf{x}}_{p}^{0}) &= 
      \SS_{p}^{e} {\mathbf{f}}_{p}^{N+1} + 
      {\mathcal{D}}_{p}^{e}{\mathbf{u}}_{p}^{N+1} + {\mathcal{B}}_{p}[{\boldsymbol\eta}^{N}] 
      + \sum_{\substack{q=1 \\ q \neq p}}^{M} {\mathcal{E}}_{pq}^{e}[{\mathbf{f}}_{q}^{N},{\mathbf{u}}_{q}^{N}], 
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:vesicle:dynamic:exp} \\
      \beta \mathcal{P}^{e}{\mathbf{x}}_{p}^{N+1} &= \mathcal{P}^{e}{\mathbf{x}}_{p}^{0},
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:inextensibility:exp} \\
      U({\mathbf{x}}) &= -\frac{1}{2}{\boldsymbol\eta}^{N}({\mathbf{x}}) + {\mathcal{E}}_{\Gamma}^{e}[{\mathbf{f}}^{N},{\mathbf{u}}^{N}]({\mathbf{x}}) + 
      {\mathcal{B}}[{\boldsymbol\eta}^{N}]({\mathbf{x}}) + \mathcal{N}_{0}[{\boldsymbol\eta}^{N}]({\mathbf{x}}), &&{\mathbf{x}} \in \Gamma, 
      \label{e:BIE:exp} \\
      {\mathbf{u}}_{p}^{N+1} &= \frac{\beta {\mathbf{x}}_{p}^{N+1}-{\mathbf{x}}_{p}^{0}}{\Delta t},
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:velocity:exp}
  \end{align}
  \label{e:explicit:method}
\end{subequations}
where $\alpha_{p} = (1+\nu_{p})/2$, and operators with a superscript $e$
are discretized at ${\mathbf{x}}^{e}$.  Note that~\eqref{e:explicit:method} has
explicit vesicle-vesicle and vesicle-boundary interactions.  That is, we
are approximating the solution of~\eqref{e:ves:dyn} with the solution of
a problem where the vesicles are fully decoupled from one another and
from the boundary.  To solve~\eqref{e:explicit:method}, we first find
${\boldsymbol\eta}^{N}$ by solving~\eqref{e:BIE:exp} and then
solve~\eqref{e:vesicle:dynamic:exp} and~\eqref{e:inextensibility:exp}
for the new position and tension of each vesicle independently of all
the others.

In high-concentration flows, a new source of stiffness arises.  As two
vesicles approach one another, ${\mathcal{E}}_{pq}$ becomes increasingly large.
A similar result holds for vesicle-boundary interactions where
${\mathcal{B}}_{p}$ introduces stiffness.  We propose the new method with
semi-implicit vesicle-vesicle and vesicle-boundary interactions (the
differences with \eqref{e:explicit:method} being highlighted in red)
\begin{subequations}
  \begin{align}
    \frac{\alpha_{p}}{\Delta t}
    (\beta {\mathbf{x}}_{p}^{N+1} - {\mathbf{x}}_{p}^{0}) &= 
      \SS_{p}^{e} {\mathbf{f}}_{p}^{N+1} + {\mathcal{D}}_{p}^{e}{\mathbf{u}}_{p}^{N+1} + 
      {\mathcal{B}}_{p}[{\boldsymbol\eta}^{N\textcolor{red}{+1}}] + \sum_{\substack{q=1 \\ q \neq p}}^{M} 
      {\mathcal{E}}_{pq}^{e}[{\mathbf{f}}_{q}^{N\textcolor{red}{+1}},
      {\mathbf{u}}_{q}^{N\textcolor{red}{+1}}],
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:vesicle:dynamic:imp} \\
      \beta \mathcal{P}^{e}{\mathbf{x}}_{p}^{N+1} &= \mathcal{P}^{e}{\mathbf{x}}_{p}^{0},
      &&{\mathbf{x}} \in \gamma_{p},
      \label{e:inextensibility:imp} \\
      U({\mathbf{x}}) &= -\frac{1}{2}{\boldsymbol\eta}^{N\textcolor{red}{+1}}({\mathbf{x}}) + 
      {\mathcal{E}}_{\Gamma}^{e}[{\mathbf{f}}^{N\textcolor{red}{+1}},
      {\mathbf{u}}^{N\textcolor{red}{+1}}]({\mathbf{x}}) + 
      {\mathcal{B}}[{\boldsymbol\eta}^{N\textcolor{red}{+1}}]({\mathbf{x}})+\mathcal{N}_{0}[{\boldsymbol\eta}^{N\textcolor{red}{+1}}]({\mathbf{x}}), &&{\mathbf{x}} \in \Gamma, 
      \label{e:BIE:imp} \\
      {\mathbf{u}}_{p}^{N+1} &= \frac{\beta {\mathbf{x}}_{p}^{N+1}-{\mathbf{x}}_{p}^{0}}{\Delta t},
      &&{\mathbf{x}} \in \gamma_{p}.
      \label{e:velocity:imp}
  \end{align}
  \label{e:implicit:method}
\end{subequations}
Unlike in~\eqref{e:explicit:method}, vesicles
in~\eqref{e:implicit:method} are implicitly coupled since the operator
${\mathcal{E}}_{pq}^{e}$ is applied to ${\mathbf{f}}_{q}^{N+1}$ and ${\mathbf{u}}_{q}^{N+1}$.
Similarly, the operator ${\mathcal{B}}_{p}$ implicitly couples the vesicles to the
solid walls.  Equation~\eqref{e:implicit:method} is still an
approximation of~\eqref{e:ves:dyn}, but its benefit is that it is
experimentally a more stable method (see Sections~\ref{s:taylorGreen}
and~\ref{s:couette}).

Since~\eqref{e:implicit:method} is fully coupled, it is too expensive to
solve without preconditioning.  Roughly speaking, the number of GMRES
iterations depends on the sizes of ${\mathcal{E}}_{pp}$ and ${\mathcal{E}}_{pq}$ and the
latter depends on the inter-vesicle proximity.  We apply preconditioned
GMRES~\cite{saad} to~\eqref{e:implicit:method} where we use the
block-diagonal preconditioner.  The number of GMRES iterations depends
on the minimum inter-vesicle distance, but this preconditioned linear
system experimentally results in mesh-independence\footnote{We
experimented with multigrid preconditioners~\cite{quaife:biros2013} but
mesh-independence was only achieved for computationally expensive
preconditioners.}.

To test the preconditioner, one step of a Taylor-Green flow is done
with the number of points per vesicle ranging from $16$ to $128$
(Table~\ref{t:GMRES:preco}).  We solve~\eqref{e:implicit:method} to a
tolerance of $10^{-8}$ with and without preconditioning.  The required
GMRES iterations for the block-diagonal solver grows almost linearly
and the results are consistent with those reported
in~\cite{rah:vee:bir}.  The implicit solver without preconditioning
requires an unacceptable number of GMRES iterations, but with
preconditioning, we see that the required number of GMRES steps is
independent of $N$ \footnote{Notice, however, that the number of GMRES
iterations will depend on how close the vesicles are to one another as
this determines the overall spectrum structure of the operator. But for
a given flow, the method is independent of $N$.}.  Using the
block-diagonal preconditioner will not result in a speedup unless we
precompute and store the block diagonal preconditioner.  We perform
this precomputation for all the experiments in Section~\ref{s:results}.

\begin{table}[htps]
\centering
\begin{tabular}{>{\centering}m{0.5cm} >{\centering}m{3cm} >{\centering}m{3cm} >{\centering}m{3cm}}
$N$ & Block- \\ Diagonal & Unpreconditioned \\ Implicit & Preconditioned
\\ Implicit \tabularnewline
\hline
16 & 46 & 315 & 18 \tabularnewline
32 & 79 & 341 & 15 \tabularnewline
64 & 164 & 882 & 15 \tabularnewline
128 & 294 & 2117 & 15 
\end{tabular}
\mcaption{The number of GMRES iterations required to take one time step
of a Taylor-Green flow using semi-implicit inter-vesicle
interactions~\eqref{e:implicit:method}.  $N$ is the number of
discretization points per vesicle, the second column is the number of
GMRES iterations to invert the block-diagonal preconditioner, the third
column is the number of unpreconditioned GMRES iterations, and the
final column is the number of preconditioned GMRES iterations.  Notice
that without preconditioning, the number of iterations quickly grows
while the preconditioned linear system is
mesh-independent.}{t:GMRES:preco}
\end{table}
We have also observed experimentally that when vesicles are close,
treating inter-vesicle interactions explicitly with second-order time
stepping is unconditionally unstable.  We have experimented with time
step sizes of varying orders of magnitude and each experiment
inevitably becomes unstable.  However, scheme~\eqref{e:implicit:method}
with second-order time stepping is stable once $\Delta t$ is small
enough for the simulation under consideration.

\subsection{Fast summation}
For high-concentration flows, the most expensive part of the simulation
is evaluating the integral operators $\SS$ and ${\mathcal{D}}$.  The single-layer
potential $\SS$ can be written as
\begin{align*}
  \frac{1}{4\pi\mu_{0}} \int_{\gamma}  \left(
    -\boldsymbol{I}\log \rho  + 
    \frac{{\mathbf{r}} \otimes {\mathbf{r}}}{\rho^{2}} \right){\mathbf{f}} ds = 
    -\frac{1}{4\pi\mu_{0}} \int_{\gamma} (\log \rho) {\mathbf{f}} ds +
    \frac{1}{4\pi\mu_{0}} \int_{\gamma}
    \left(\frac{{\mathbf{r}} \cdot {\mathbf{f}}}{\rho^{2}}\right){\mathbf{r}} ds \\
  =-\frac{1}{4\pi\mu_{0}} \int_{\gamma}(\log \rho) {\mathbf{f}} ds + 
    \frac{1}{4\pi\mu_{0}} \int_{\gamma}   
    \frac{{\mathbf{r}}}{\rho^{2}} ({\mathbf{x}} \cdot {\mathbf{f}}) ds_{{\mathbf{y}}} - 
    \frac{1}{4\pi\mu_{0}} \int_{\gamma} 
    \frac{{\mathbf{r}}}{\rho^{2}} ({\mathbf{y}} \cdot {\mathbf{f}}) ds_{{\mathbf{y}}},
\end{align*}
which can be computed with three Laplace FMM applications.  We also use
the FMM to evaluate the single-layer potential at arbitrary target
locations by assigning a charge of zero to these locations.
Unfortunately, the double-layer potential can not be decomposed into a
sequence of Laplace single- and double-layer potentials.  Future work
involves using the kernel-independent FMM~\cite{ying-biros-zorin-03} to
accelerate evaluating the double-layer potential.


\subsection{Near-singular integration}\label{s:near-singular}
Another difficulty with high-concentration simulations is the
evaluation of layer potentials at points close to a boundary.  Integral
operators ${\mathcal{D}}$ and $\SS$ can be approximated with high accuracy using
the trapezoid rule if the target point is sufficiently far from the
source points.  As the target point approaches the vesicle,
more points are required to achieve the same accuracy, and, the
required number of points for a fixed error grows without
bound (Appendix~\ref{A:AppendixA}).  A near-singular integration
strategy that guarantees a uniform accuracy for all target points is
required.  

We adopt the strategy from~\cite{ying-biros-zorin06} since it is not
depend on the nature of the integrand and it is simple to implement.
Suppose that a vesicle $\gamma$ is discretized with $N$ points and that
the resulting arclength term is of size $h$.  Let $d({\mathbf{x}},\gamma) =
\inf_{{\mathbf{y}} \in \gamma} \|{\mathbf{x}} - {\mathbf{y}}\|$ be the distance from ${\mathbf{x}}$ to
$\gamma$.  We define the far zone of $\gamma$ as $\Omega_{1}=\{{\mathbf{x}}
\:|\: d({\mathbf{x}},\gamma) \geq h\}$, and the near zone of $\gamma$ as
$\Omega_{0}=\{{\mathbf{x}} \:|\: d({\mathbf{x}},\gamma) < h\}.$ The near and far zones
can be constructed efficiently using a quadtree-structure, similar to
the near-field and far-field interactions in fast multipole codes.

For ${\mathbf{x}} \in \Omega_{1}$, the trapezoid rule with $N^{3/2}$ points
guarantees that the error for $\SS$ is $\mathcal{O}(h^{M/2-2})$ and for
${\mathcal{D}}$ is $\mathcal{O}(h^{M/2-4})$, where the density function (in our
case, either ${\mathbf{f}}$ or ${\boldsymbol\eta}$) belongs to $C^{M}$.  These error bounds
are proved in Appendix~\ref{A:AppendixA}.  For a target point ${\mathbf{x}}$ in
the near zone, we first find the closest boundary point, ${\mathbf{x}}_{0}$,
using two applications of Newton's method with the initial guess being
the closest discretization point to ${\mathbf{x}}$.  We then place $m$
interpolation points at the points ${\mathbf{x}}_{j} = {\mathbf{x}}_{0} + j\beta h
({\mathbf{x}}-{\mathbf{x}}_{0})/\|{\mathbf{x}}-{\mathbf{x}}_{0}\|$ for $j=0,\ldots,m$, where $\beta$ is
slightly larger than $1$ to guarantee that all the interpolation points
are in the far zone (left plot of Figure~\ref{f:nearsing:conv}).  The
layer potential is evaluated at ${\mathbf{x}}_{0}$ using a local interpolant with
$N_{\mathrm{int}}$ points on $\gamma$ and evaluated at ${\mathbf{x}}_{j},$
$j=1,\ldots,m$ using the $N^{3/2}$-point trapezoid rule.  Then, a
one-dimensional Lagrange interpolation assigns a value to the layer
potential at ${\mathbf{x}}$.  Note that the proposed scheme can become more
accurate by using interpolation schemes that compute derivatives of the
solution and building more accurate polynomial interpolations.  This,
however, would significantly complicate the implementation.  Pseudocode
for this algorithm is given in Algorithm~\ref{a:algorithm2}.  In
Appendix~\ref{A:AppendixA}, we show that the error of our algorithm
applied to the single-layer potential $\SS$ is 
\begin{align*}
  \mathcal{O}\left(h^{\min\left(N_{\mathrm{int}}-1,
      m,\frac{M}{2}-2\right)}\right),
\end{align*}
and applied to the double-layer potential ${\mathcal{D}}$ is
\begin{align*}
  \mathcal{O}\left(h^{\min\left(N_{\mathrm{int}}-1,
      m,\frac{M}{2}-4\right)}\right).
\end{align*}

\begin{algorithm}
\begin{algorithmic} 
\FORALL[Loop over far zone]{${\mathbf{x}} \in \Omega_{1}$} 
\STATE ${\mathbf{u}} = \mathtt{trapezoid}(\gamma,{\mathbf{x}},{\mathbf{f}})$
\COMMENT{Apply the $N^{3/2}$-point trapezoid rule}
\ENDFOR 
\FORALL[Loop over near zone]{${\mathbf{x}} \in \Omega_{0}$}
\STATE ${\mathbf{x}}_{0} = \argmin d({\mathbf{x}},\gamma)$
\COMMENT{Find the point on $\gamma$ closest to ${\mathbf{x}}$}
\STATE ${\mathbf{u}}_{0} = \mathtt{interpolate}
    ({\mathbf{x}}_{0},{\mathbf{u}}_{1},\ldots,{\mathbf{u}}_{N_{\mathrm{int}}})$
\COMMENT{Interpolate ${\mathbf{u}}$ at ${\mathbf{x}}_{0}$}
\FOR{$j=1,\ldots,m$}
\STATE ${\mathbf{x}}_{j} = {\mathbf{x}}_{0} + jh\beta({\mathbf{x}}-{\mathbf{x}}_{0})/\|{\mathbf{x}}-{\mathbf{x}}_{0}\|$
\COMMENT{Define Lagrange interpolation points}
\STATE ${\mathbf{u}}_{j} = \mathtt{trapezoid}(\gamma,{\mathbf{x}}_{j},{\mathbf{f}})$
\COMMENT{Evaluate layer potential at ${\mathbf{x}}_{j}$}
\ENDFOR 
\STATE ${\mathbf{u}} = \mathtt{interpolate}({\mathbf{x}},{\mathbf{u}}_{0},\ldots,{\mathbf{u}}_{m})$
\COMMENT{Compute ${\mathbf{u}}$ using Lagrange interpolation}
\ENDFOR 
\RETURN ${\mathbf{u}}$ \COMMENT{Return the layer potential at all target points}
\end{algorithmic}
\mcaption{LayerPotential($\gamma$,  ${\mathbf{f}}$, ${\mathbf{u}}$, $\Omega_{0}$,
$\Omega_{1}$)} {a:algorithm2}
\end{algorithm}

\begin{figure}[htps]
  \begin{center}
    \begin{minipage}{0.45\textwidth}
    \input{near-sing.tikz}
    \end{minipage}
    \begin{minipage}{0.45\textwidth}
    \input{nearSingErrors.tikz} 
    \end{minipage} 
  \end{center} 
  \mcaption{\emph{Left:} Lagrange interpolation points in $\Omega_{1}$
  and $\gamma$ are used to approximate layer potentials at ${\mathbf{x}} \in
  \Omega_{0}$.  The value at ${\mathbf{x}}_{0}$ uses an interpolant from
  neighboring points on the boundary.  All other interpolation points
  are in $\Omega_{1}$ so that the $N^{3/2}$-point trapezoid rule is
  sufficiently accurate.  \emph{Right:} Here we report convergence
  rates. For the blue and the black curves, the expected order of
  convergence is $5$. We obtain $4.7$ and $5.0$, respectively. As the
  green curve corresponds to a linear flow, Lagrange interpolation gives
  exact values.  Hence, machine accuracy is reached once $N$ is
  sufficiently large.}{f:nearsing:conv}
\end{figure}

We test our near-singular integration with three examples.  In all the
examples, 32 target points are placed in $\Omega_{0}$ for
$N=16,\ldots,512$.  The boundary is a $3:2$ ellipse, the boundary values
are interpolated with $N_{\mathrm{int}}=6$ points, and the Lagrange
interpolation uses $m=6$ points.  Since $m = N_{\mathrm{int}}$, we are
using one more point to interpolate in the direction coming off of
$\gamma$.  We do this because the interpolant that is coming off of the
vesicle is always evaluated near the first interpolation point.  Thus,
it may suffer from the Runge phenomenon (we have observed this
experimentally).  This is not the case for the local interpolation used
to find ${\mathbf{u}}({\mathbf{x}}_{0})$ since ${\mathbf{x}}_{0}$ is located near the middle of the
local interpolation points.  For the first two examples, we compute the
density function of $\SS$ for the two solutions of the Stokes equations:
\begin{align*}
  &{\mathbf{u}}({\mathbf{x}})=(x,-y), \\
  &{\mathbf{u}}({\mathbf{x}})=\left(-\log|{\mathbf{x}}| + \frac{{\mathbf{x}} \otimes {\mathbf{x}}}{|{\mathbf{x}}|^{2}}\right)
  \left[ \begin{array}{c}
    1 \\ 1 
  \end{array}\right].
\end{align*}
For the third example, we pick an arbitrary density function and compare
the results from near-singular integration with the results from an
over-refined trapezoid rule.  The error for all three examples is
reported in~Figure~\ref{f:nearsing:conv}.  For the extensional flow ${\mathbf{u}}
= (x,-y)$, the Lagrange interpolation is exact since the flow is linear
and the error decays faster than the other examples.  For the other
examples, the error is close to fifth-order which agrees with the order
of accuracy of the near-singular integration

If we have $M$ vesicles and $N$ points per vesicle, the major costs of
the near-singular integration algorithm are:
\begin{itemize}
  \item \emph{Upsampling each vesicle from $N$ to $N^{3/2}$ points:}
  Using the FFT, this requires $\mathcal{O}(MN^{3/2}\log N)$ operations.
  \item \emph{Applying the $N^{3/2}$-point quadrature rule to
  $\Omega_{1}$:} With the FMM, this requires $\mathcal{O}(MN^{3/2} + MN)$
  operations, and without the FMM, this requires $\mathcal{O}(M^{2}N^{5/2})$
  operations.  
  \item \emph{Evaluating the layer potential in $\Omega_{0}$:} Assuming
  there are $\mathcal{O}(1)$ points in the near zone, this requires
  $\mathcal{O}(MN^{3/2})$ operations.
\end{itemize}
The algorithm can be further sped up by introducing the intermediate and
far zones~\cite{ying-biros-zorin06}.  The intermediate zone is
$\Omega_{1}=\{{\mathbf{x}} \:|\: d({\mathbf{x}},\gamma) \in (\sqrt{h},h]\}$ and the far
zone is $\Omega_{2}=\{{\mathbf{x}} \:|\: d({\mathbf{x}},\gamma) \geq \sqrt{h}\}$.  Then,
the same accuracy can be achieved by using the $N^{3/2}$-point trapezoid
rule in $\Omega_{1}$ and the $N$-point trapezoid rule in $\Omega_{2}$.
While this would reduce the constant of the complexity, it does not
change the overall order and complicates the implementation.


\subsection{Collision detection}\label{s:collision}
Even with semi-implicit treatment of inter-vesicle interactions, due
to time discretization errors, the vesicles can come so close that they
collide or cross a solid wall.\footnote{The streamlines of Stokes
equations never intersect.  Therefore, it is physically impossible for
vesicles to collide with one another or with solid walls in finite
time.  Lubrication theory can be used to study the forces between
nearly touching vesicles, and therefore compute estimates of the
distance between the vesicles, but we choose not to discuss these
forces here.}  There has been a lot of work in collision
detection~\cite{jimenez-e13,redon-e04}, but they are not spectrally
accurate or require complex surface-surface intersections.  Most
methods construct piecewise approximations to the boundaries to
simplify the detection of collisions.  Since our vesicles and solid
walls are periodic, we can use some nice properties of layer potentials
to create a spectrally accurate collision detector.  We use the classic
potential theory result~\cite{kellogg}
\begin{align}
  I_{\Gamma}({\mathbf{x}}) = \frac{1}{2\pi}\int_{\Gamma}\pderiv{}{{\mathbf{n}}_{{\mathbf{y}}}}\log|{\mathbf{x}} - {\mathbf{y}}| ds_{{\mathbf{y}}} = 
  \begin{cases}
    1 & {\mathbf{x}} \in \Omega, \\ 
    \frac{1}{2} & {\mathbf{x}} \in \Gamma, \\ 
    0 & {\mathbf{x}} \in \mathbf{R}^{2} \backslash \overline{\Omega}, 
  \end{cases}
  \label{e:layerpotentials}
\end{align}
where $\Omega \subset \mathbf{R}^{2}$ is bounded with a smooth boundary
$\Gamma$, and ${\mathbf{n}}_{{\mathbf{y}}}$ is the outward pointing normal.  We
apply~\eqref{e:layerpotentials} to the configuration of the vesicles
and the solid walls by computing $I_{X}({\mathbf{x}}) = I_{\gamma_{1}}({\mathbf{x}}) +
\cdots + I_{\gamma_{M}} + I_{\Gamma}$ for all ${\mathbf{x}} \in \gamma$.  The
maximum value of $I_{X}({\mathbf{x}})$ is $3/2$ for confined flows and $1/2$ for
unconfined flows if vesicles have not crossed and have not left the
domain.  The maximum value of $I_{X}({\mathbf{x}})$ will be at least $3/2$ for
unconfined flows, and at least $5/2$ for confined flows if vesicles
have crossed or left the domain.  We note that $I_{X}({\mathbf{x}})$ can be
evaluated accurately using our near-singular integration scheme, and in
linear time using the FMM. 

We can use $I_{X}({\mathbf{x}})$ to create a simple adaptive time-stepping
method.  At each time step, we check for collisions.  If a collision has
occurred, we initialize the simulation at a previous time step with a
time step half the size.  We are currently implementing a more rigorous
adaptive time stepping method that uses ideas from spectral deferred
correction methods~\cite{dut:gre:rok,hua:jia:min,min}.  At each time
step, we estimate the local truncation error and then try to commit a
constant amount of error per time step by adjusting the time step size.
We leave the details and the implementation as future work.  Future work
also includes adaptivity in space.  If vesicles are sufficiently far
from one another, their hydrodynamic interaction can be represented with
a moderate number of points, but as they approach one another, the
hydrodynamic interaction strengthens and this requires a finer spatial
resolution.



\subsection{Shear Flow}
Here we consider two vesicles in the shear flow $(y,0)$.  We initially
place the left vesicle's center slightly above the x-axis and place the
right vesicle's center exactly on the x-axis.  Initially, we expect the
right vesicle to undergo tank-treading~\cite{kantsler,misbah2006} while
the other vesicle will travel from left to right.  Here we study the
effect of the reduced area and the viscosity contrast on the dynamics
and errors.  We expect vesicles with larger viscosity contrast and
larger reduced areas to behave more like rigid bodies.  In particular,
we expect them to come closer to one another.
Figure~\ref{f:shear:reducedAreas} indicate that a higher viscosity
contrast results in the vesicles coming closer together.  Moreover, for
larger reduced areas, the effect of the viscosity contrast is less
pronounced.  We also see that vesicles with a larger reduced area come
closer together.

\begin{figure}[htps]
\centering
\begin{tabular}{m{1cm}CCCCC}
$\nu=1$ &
\input{shearOrder1RA1viscCont1Time1.tikz} &
\input{shearOrder1RA1viscCont1Time2.tikz} &
\input{shearOrder1RA1viscCont1Time3.tikz} &
\input{shearOrder1RA1viscCont1Time4.tikz} &
\input{shearOrder1RA1viscCont1Time5.tikz} \\
$\nu=4$ &
\input{shearOrder1RA1viscCont4Time1.tikz} &
\input{shearOrder1RA1viscCont4Time2.tikz} &
\input{shearOrder1RA1viscCont4Time3.tikz} &
\input{shearOrder1RA1viscCont4Time4.tikz} &
\input{shearOrder1RA1viscCont4Time5.tikz} \\
$\nu=1$ &
\input{shearOrder1RA2viscCont1Time1.tikz} &
\input{shearOrder1RA2viscCont1Time2.tikz} &
\input{shearOrder1RA2viscCont1Time3.tikz} &
\input{shearOrder1RA2viscCont1Time4.tikz} &
\input{shearOrder1RA2viscCont1Time5.tikz} \\
$\nu=4$ &
\input{shearOrder1RA2viscCont4Time1.tikz} &
\input{shearOrder1RA2viscCont4Time2.tikz} &
\input{shearOrder1RA2viscCont4Time3.tikz} &
\input{shearOrder1RA2viscCont4Time4.tikz} &
\input{shearOrder1RA2viscCont4Time5.tikz}
\end{tabular}
\begin{tabular}{cc}
\input{shearRA1GapSizeOrder1.tikz} &
\input{shearRA2GapSizeOrder1.tikz}
\end{tabular}
\mcaption{Two vesicles are submerged in a shear flow.  In the first and
second rows, the vesicles have reduced area 0.65, and in the third and
fourth rows, the vesicles have reduced area 0.99.  Each simulation is
done with viscosity contrasts $\nu=1$ and $\nu=4$.  The bottom plots
show the distance between the vesicles for both reduced areas and both
viscosity contrasts.  Comparing semi-implicit inter-vesicle interactions
with both first- and second-order time stepping, the distance between
the vesicles never differ by more than 2.2\% for both reduced areas and
both viscosity contrasts.}{f:shear:reducedAreas} 
\end{figure}

We now do a convergence study for both values of the reduced area and
viscosity contrast.  We vary the number of points per vesicle, the
number of time steps, and the time integrator.  We report the errors in
area and length at the time horizon $t=12$ in Tables~\ref{t:shear:RA65}
and~\ref{t:shear:RA99}.  First, we see that first- and second-order
convergence is achieved for all combinations of reduced area and
viscosity contrast.  Second, semi-implicit inter-vesicle interactions
do not improve the accuracy of first-order time stepping, but, we will
see later that this time integrator is more stable.

\begin{table}[htps]
\centering
\begin{tabular}{ccc|cc|cc} 
Time & & & \multicolumn{2}{c|}{$\nu=1$} & \multicolumn{2}{c}{$\nu=4$} \\
Integrator   & $N$ & $\Delta t$ & $e_{A}$   & $e_{L}$   & 
$e_{A}$   & $e_{L}$ \\ 
\hline

Explicit (1) & 32  & 0.04       & $6.73e-2$ & $5.21e-2$ &
$4.07e-2$ & $2.41e-2$ \\ 

Explicit (1) & 64  & 0.02       & $3.36e-2$ & $2.62e-2$ &
$2.02e-2$ & $1.20e-2$ \\ 

Explicit (1) & 128 & 0.01       & $1.68e-2$ & $1.31e-2$ &
$1.01e-2$ & $5.99e-2$ \\ 
\hline

Implicit (1) & 32  & 0.04       & $6.73e-2$ & $5.21e-2$ &
$4.06e-2$ & $2.40e-2$ \\ 

Implicit (1) & 64  & 0.02       & $3.36e-2$ & $2.62e-2$ &
$2.02e-2$ & $1.20e-2$ \\ 

Implicit (1) & 128 & 0.01       & $1.68e-2$ & $1.31e-2$ &
$1.02-2$ & $5.98e-3$ \\ 
\hline

Implicit (2) & 32  & 0.04       & $1.59e-4$ & $3.21e-4$ &
$2.69e-5$ & $6.33e-5$ \\ 

Implicit (2) & 64  & 0.02       & $3.87e-6$ & $6.52e-5$ &
$3.78e-6$ & $1.20e-5$ \\ 

Implicit (2) & 128 & 0.01       & $2.79e-7$ & $1.47e-5$ &
$3.85e-7$ & $2.28e-6$ 
\end{tabular}
\mcaption{The errors in area and length at $t=12$ for different time
integrators applied to two vesicles of reduced area $0.65$ in a shear
flow.  The parenthetic values are the time stepping order.  Results are
presented for no viscosity contrast and a viscosity contrast of
4.}{t:shear:RA65}

\begin{tabular}{ccc|cc|cc} 
Time & & & \multicolumn{2}{c|}{$\nu=1$} & \multicolumn{2}{c}{$\nu=4$} \\
Integrator   & $N$ & $\Delta t$ & $e_{A}$   & $e_{L}$   & 
$e_{A}$   & $e_{L}$ \\ 
\hline

Explicit (1) & 32  & 0.04       & $1.05e-1$ & $5.48e-2$ &
$9.67e-2$ & $4.94e-2$ \\ 

Explicit (1) & 64  & 0.02       & $5.15e-2$ & $2.70e-2$ &
$4.78e-2$ & $2.45e-2$ \\ 

Explicit (1) & 128 & 0.01       & $2.55e-2$ & $1.34e-2$ &
$2.37e-2$ & $1.22e-2$ \\ 
\hline

Implicit (1) & 32  & 0.04       & $1.05e-1$ & $5.48e-2$ &
$9.69e-2$ & $4.93e-2$ \\ 

Implicit (1) & 64  & 0.02       & $5.16e-2$ & $2.70e-2$ &
$4.77e-2$ & $2.45e-2$ \\ 

Implicit (1) & 128 & 0.01       & $2.55e-2$ & $1.34e-2$ &
$2.36e-2$ & $1.22e-2$ \\ 
\hline

Implicit (2) & 32  & 0.04       & $2.40e-5$ & $1.13e-4$ &
$5.85e-5$ & $4.68e-5$ \\ 

Implicit (2) & 64  & 0.02       & $1.65e-6$ & $1.50e-5$ &
$7.05e-6$ & $5.86e-6$ \\ 

Implicit (2) & 128 & 0.01       & $1.56e-7$ & $1.92e-6$ &
$8.62e-7$ & $7.26e-7$ 

\end{tabular}

\mcaption{The errors in area and length at $t=12$ for different time
integrators applied to two vesicles of reduced area $0.99$ in a shear
flow.  The parenthetic values are the time stepping order.  Results are
presented for no viscosity contrast and a viscosity contrast of
4.}{t:shear:RA99}
\end{table}

\subsection{Taylor-Green Flow}
\label{s:taylorGreen}
We consider the Taylor-Green flow ${\mathbf{u}} = (\sin x \cos y,-\cos x \sin
y)$ with 9 vesicles each with $N=64$ points of reduced area 0.65 placed
inside the periodic cell $(0,2\pi)^{2}$ (left plot of
Figure~\ref{f:taylorGreen:collision}).  We consider first-order IMEX
Euler with both semi-implicit and explicit inter-vesicle interactions.
We find the largest time step of both methods so that the vesicles do
not cross before a time horizon of $T = 1$.  With explicit
interactions, the vesicles cross if $74$ time steps are taken but not
if $75$ time steps are taken.  However, with semi-implicit
interactions, the simulation successfully completes even with one time
step.  Of course the error is too large and a time step size this large
is unreasonable.

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{ccc}
    \input{taylorGreenSnapsExplicitTime1.tikz} &
    \input{taylorGreenSnapsExplicitTime2.tikz} &
    \input{taylorGreenSnapsImplicitTime2.tikz} \\
    &
    \input{taylorGreenSnapsExplicitTime2Zoomed.tikz} &
    \input{taylorGreenSnapsImplicitTime2Zoomed.tikz}
  \end{tabular}
\end{center}
\mcaption{We compare first-order time stepping with explicit and
semi-implicit inter-vesicle interactions for a time step size of 0.02
and 64 points per vesicle with the initial configuration in the
left-most plot.  With explicit interactions (middle plot), the vesicles
cross after only 13 time steps.  However, with semi-implicit
interactions (right plot), the vesicles do not cross.  The bottom plots
show magnifications of the region where the explicit solver first
crosses.}{f:taylorGreen:collision}
\end{figure}

For this example, we found that using explicit interactions with 75
time steps requires approximately the same amount of CPU time as using
semi-implicit interactions with 30 time steps.  Therefore, if one
permits the error to be approximately 2.5 times larger, using
semi-implicit inter-vesicle interactions is justified.

We do a convergence study for first- and second-order time stepping
with semi-implicit inter-vesicle interactions.  If these interactions
are treated explicitly, the coarsest simulation results in the vesicles
crossing.  We achieve the desired first-order convergence, but
second-order convergence is not achieved at these resolutions.
However, we ran the simulation with $\Delta t = 3.13e-4$ and $\Delta t
= 1.56e-4$ and this is where second-order convergence is first
observed; the error in area reduced by $4.67$ and the error in length
reduced by $3.96$.

We have also observed a qualitative difference in the
behaviour of the errors of first- and second-order time stepping.  The
errors of the first-order methods continue to grow throughout the
simulation, while for second-order time stepping, the errors
plateau.  Thus, given a fixed time step size, second-order time stepping
can simulate longer time horizons than first-order time stepping.  We
demonstrate this behaviour in the bottom plot of
Figure~\ref{f:taylorGreen:summary} where we plot the errors in area and
length for the first- and second-order time stepping.  Snapshots of the
simulation are in the top plots of Figure~\ref{f:taylorGreen:summary}.

\begin{table}[htps]
\begin{center}
\begin{tabular}{cc|cc|cc}
& & \multicolumn{2}{c|}{First-order} & \multicolumn{2}{c}{Second-order} \\
$N$ & $\Delta t$ & $e_{A}$ & $e_{L}$ & $e_{A}$ & $e_{L}$ \\
\hline
32  & $2e-2$ & $5.24e-2$ & $3.03e-2$ & $5.25e-4$ & $3.98e-3$ \\
64  & $1e-2$ & $2.77e-2$ & $1.64e-2$ & $2.08e-4$ & $1.96e-3$ \\
128 & $5e-3$ & $1.43e-2$ & $8.68e-3$ & $7.76e-5$ & $9.42e-4$ \\
256 &$2.5e-3$& $7.24e-3$ & $4.53e-3$ & $2.73e-5$ & $4.20e-4$ \\
\end{tabular}
\mcaption{The errors in area and length at $t=5$ for a Taylor-Green
flow.  All runs were done with semi-implicit inter-vesicle interactions
since explicit interactions with $N=32$ and $\Delta t = 2e-2$ results
in crossing vesicles.}{t:taylorGreen:convergence} 
\end{center}
\end{table}

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{cccccc}
    \input{taylorGreenSnapsTime1.tikz} &
    \input{taylorGreenSnapsTime2.tikz} &
    \input{taylorGreenSnapsTime3.tikz} &
    \input{taylorGreenSnapsTime4.tikz} &
    \input{taylorGreenSnapsTime5.tikz} &
    \input{taylorGreenSnapsTime6.tikz}
  \end{tabular}
  \begin{tabular}{cc}
    \input{taylorGreenErrorsOrder1.tikz} & 
    \input{taylorGreenErrorsOrder2.tikz} 
  \end{tabular}
\end{center}
\mcaption{The top plots show several snapshots of a Taylor-Green
simulation with a long time horizon using second-order time stepping.
In the bottom plots are the errors in area and length for the first- and
the second-order time integrators.  Notice that the errors for the
first-order method continue to grow and the simulation stops near $t=20$
(the simulation is stopped when a 10\% error has been incurred), whereas
the errors for the second-order method plateau and the simulation is
able to run much longer.}{f:taylorGreen:summary} 
\end{figure}



\subsection{Extensional Flow}
We consider two vesicles with reduced area $0.65$ symmetrically placed
around the origin (top plots of Figure~\ref{f:extensional:summary}).
The background velocity is given by ${\mathbf{u}} = (-x,y)$.  We first run the
simulation with and without near-singular integration.  We compute the
distance between the vesicles, and the errors in area and length
(bottom plots of Figure~\ref{f:extensional:summary}).  We see that when
$t \approx 5$, the simulation without near-singular integration begins
to introduce large errors in the area of the vesicles.  This indicates
that the approximation of the single-layer potential is not adequate.
Then, around $t \approx 12$, the vesicles are too close for the
$N$-point trapezoid rule and the simulations fails.  This agrees with
our error estimate from Appendix~\ref{A:AppendixA}  since for this
simulation, an arclength term is approximately 0.24.  However, with
near-singular integration, the simulation successfully runs to
completion.

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{ccccc}
  \input{extensionalSnapsTime1.tikz} &
  \input{extensionalSnapsTime2.tikz} &
  \input{extensionalSnapsTime3.tikz} &
  \input{extensionalSnapsTime4.tikz} &
  \input{extensionalSnapsTime5.tikz} 
  \end{tabular}
  \begin{tabular}{cc}
  \input{extensionalGapSize.tikz} &
  \input{extensionalErrors.tikz}
  \end{tabular}
\end{center}
\mcaption{Top: Two vesicles discretized with $N=32$ points at several
time steps using near-singular integration.  Bottom left: The distance
between the two vesicles both with (solid) and without (dashed)
near-singular integration.  Bottom right: The error in area (blue) and
length (red) with (solid) and without (dashed) near-singular
integration.}{f:extensional:summary} 
\end{figure}

We report the errors in area and length at $t=24$ in
Table~\ref{t:extensional:convergence}.  Also reported is the error in
the distance between the vesicles at $t=24$, for which we take the
``true'' distance from an overrefined simulation.  As before, we
achieve first-order convergence and there is little difference between
the accuracy when comparing explicit and semi-implicit inter-vesicle
interactions.  We also achieve second-order convergence, but it
requires a finer time step size when compared to the first-order
methods before the asymptotic rate is achieved.

In Table~\ref{t:extensional:convergence}, we also list the required CPU
time per time step relative to the smallest simulation.  When the
number of points is doubled from 32 to 64, the amount of work per time
step goes up by less than a factor of two.  The increase in CPU time
when $N$ is doubled from 64 to 128 is attributed to extra work required
by near-singular integration since, to simplify the implementation, we
only work with powers of 2.  Therefore, we upsample 64 to 256, but 128
is upsampled to 4096.  We also see that semi-implicit inter-vesicle
interactions are computationally not much more expensive.  This is
because a lot of work is spent precomputing the block-diagonal
preconditioner which is done once per time step for all the time
integrators.

\begin{table}[htps] 
\begin{centering} 
\begin{tabular}{ccccccc} 
Integrator & $N$ &   $\Delta t$ & $e_{A}$ & $e_{L}$ & $e_{\mathrm{gap}}$ & CPU Time \\ \hline
Explicit (1) & $32$  & $0.04$ & $1.99e-4$ & $6.53e-4$ & $1.97e-1$ & --- \\  
Explicit (1) & $64$  & $0.02$ & $6.34e-5$ & $3.42e-4$ & $3.73e-3$ & $1.15$ \\ 
Explicit (1) & $128$ & $0.01$ & $3.13e-5$ & $1.78e-4$ & $1.16e-3$ & $2.83$ \\ 
\hline
Implicit (1) & $32$  & $0.04$ & $2.12e-4$ & $6.64e-4$ & $1.92e-1$ & $2.65$ \\  
Implicit (1) & $64$  & $0.02$ & $6.43e-5$ & $3.46e-4$ & $2.23e-3$ & $2.72$ \\ 
Implicit (1) & $128$ & $0.01$ & $3.15e-5$ & $1.79e-4$ & $4.35e-4$ & $7.40$ \\ 
\hline
Implicit (2) & $32$  & $0.04$ & $3.26e-4$ & $1.38e-4$ & $3.96e-2$ & $2.43$ \\  
Implicit (2) & $64$  & $0.02$ & $4.14e-6$ & $4.45e-5$ & $8.32e-2$ & $2.89$ \\ 
Implicit (2) & $128$ & $0.01$ & $5.08e-7$ & $1.27e-5$ & $2.69e-2$ & $8.13$ \\
Implicit (2) & $128$ & $0.005$& $1.02e-7$ & $3.13e-6$ & $1.27e-4$ & $7.70$ \\
\end{tabular}
\mcaption{The errors in area, length, and the error in the gap size at
$t=24$, and the CPU time per time step for the extensional flow.  The
parenthetic values are the order of the time integrator.  The error in
the gap size is found by computing a ``true'' solution with $N=256$ points per
vesicle and a time step size so that the error in area is $9.0e-9$ and
in length is $7.3e-11$.  The CPU times are relative to the simulation
with $N=32$, $\Delta t = 0.04$ which took approximately
$2.62e-1$ seconds per time step.}{t:extensional:convergence}
\end{centering}
\end{table}



\subsection{Stenosis}

We consider a single vesicle of reduced area $0.65$ in a constricted
tube.  The solid wall is parameterized as ${\mathbf{x}}(\theta) =
10r(\theta)\cos(\theta)$ and $y(\theta) =
3r(\theta)\sin(\theta)\eta(\theta)$, where $r(\theta) =
(\cos(\theta)^{8} + \sin(\theta)^{8})^{-1/8}$, and
\begin{align*}
  \eta(\theta) = \left\{
    \begin{array}{cl}
      \frac{1-0.6\cos(x)}{1.6} & |x(\theta)| \leq \pi, \\
      1 & |x(\theta)| > \pi.
    \end{array}
  \right.
\end{align*}
For this geometry, the error of the trapezoid rule is $\mathcal{O}(10^{-4})$
with $N_{\mathrm{wall}}=256$ points and is $\mathcal{O}(10^{-6})$ with
$N_{\mathrm{wall}}=512$ points.  We impose a Gaussian-shaped boundary
condition at the intake and outake of the tube which is plotted along
with the geometry in Figure~\ref{f:stenosis:snaps}.  The vesicle's
initial height is $2.3$ times larger than the size of the constriction.
We report a convergence study for second-order time stepping in
Table~\ref{t:stenosis:errors} and we exceed the expected rate of
convergence.  We note that for this example, $N_{\mathrm{wall}}$ can not
be less than 256.  If it is, the Lagrange interpolation points required
by near-singular integration for wall-to-vesicle interactions are on
both sides of the solid wall.  Since the double-layer potential has a
jump across the solid wall, this creates a discontinuity in the velocity
field and the method becomes unstable.

This example is ideally situated for adaptive time stepping.  As the
vesicle passes through the constriction, large jumps in the errors are
committed.  In the simulation with $\Delta t = 5.0e-3$, the error in
length jumps by four orders of magnitude, and the error in area jumps
by three orders of magnitude near $t=3$.  These jumps are less severe
in the simulation with $\Delta t = 2.5e-3$, which has a jump of one
order of magnitude in length and two orders of magnitude in area.

\begin{figure}[htps] 
\centering 
\begin{tabular}{ccc}
\input{stenosisSnapsTime1.tikz} &
\input{stenosisSnapsTime2.tikz} &
\input{stenosisSnapsTime3.tikz} \\
\input{stenosisSnapsTime4.tikz} &
\input{stenosisSnapsTime5.tikz} &
\input{stenosisSnapsTime6.tikz} \\
\input{stenosisSnapsTime7.tikz} &
\input{stenosisSnapsTime8.tikz} &
\input{stenosisSnapsTime9.tikz}
\end{tabular} 
\mcaption{Snapshots of a single vesicle in a constricted tube.  The
boundary condition on the solid wall (plotted in blue) at the inlet and
outlet is Gaussian-shaped, and on the remainder of the solid wall is
zero.  The vesicle is initially $2.3$ times larger than the size of the
constriction.  Results from~\cite{rah:vee:bir} were only able to
simulate vesicles that were $1.12$ times larger than the
constriction.}{f:stenosis:snaps}
\end{figure}

\begin{table}[htps] 
\begin{centering} 
  \begin{tabular}{cccccc} 
    $N$ & $N_{\mathrm{wall}}$ & $\Delta t$ & 
      $e_{A}$ & $e_{L}$ \\ \hline 
    128 & 256 & $5.0e-3$ & $6.55e-4$ & $1.17e-2$ \\ 
    256 & 512 & $2.5e-3$ & $9.59e-7$ & $1.05e-6$  
  \end{tabular} 
  \mcaption{The errors in area and length at $t=8$ for the second-order
  stenosis flow.}{t:stenosis:errors}
  \end{centering}
\end{table}

We also report the minimum distance, and the time when the minimum
distance occurs, between the vesicle and the solid wall in
Table~\ref{t:stenosis:gapSize}.  In order to focus the error on the time
step size, we take a fine spatial resolution of $N=256$ and
$N_{\mathrm{wall}}=512$.  We report the minimum over time of the
distance 
\begin{align*}
  d(\gamma,\Gamma) = \inf_{{\mathbf{x}} \in \gamma, {\mathbf{y}} \in \Gamma} \|{\mathbf{x}} - {\mathbf{y}}\|
\end{align*}
rather than the error since an exact solution is not available.  The
time integrators predict different locations, $t_{\min}$, where the
minimum occurs, but they both are converging.  For the first-order
results, if we compare the minimum distances with the distance of the
most accurate second-order simulation at $t=5.26$, which is $9.6369e-2$,
then we see that first-order convergence has been achieved.  Without a
much more accurate numerical solution, we are unable to compute
convergence rates for the second-order results.  However, it appears
that we have resolved the distance up to at least three digits of
accuracy.

\begin{table}[htps] 
  \begin{centering} 
  \begin{tabular}{c|cc|cc} 
  & \multicolumn{2}{c|}{First-order} &
    \multicolumn{2}{c}{Second-order} \\
    $\Delta t$ & $t_{\min}$ & $d(\Gamma,\gamma)$ &  
      $t_{\min}$ & $d(\Gamma,\gamma)$  \\ \hline
    $1.0e-2$ & $4.37$ & $1.0318e-1$ & $4.34$ & $1.0196e-1$ \\
    $5.0e-3$ & $4.41$ & $1.0262e-1$ & $4.39$ & $9.8055e-2$ \\
    $2.5e-3$ & $4.42$ & $1.0143e-1$ & $5.77$ & $9.5809e-2$ \\
    $1.3e-3$ & $5.26$ & $9.9541e-2$ & $5.77$ & $9.5804e-2$ \\
    $6.3e-4$ & $5.26$ & $9.7938e-2$ & $5.77$ & $9.5803e-2$ \\ 
    $3.1e-4$ & $5.26$ & $9.7137e-2$ & $5.77$ & $9.5802e-2$ 
  \end{tabular}
 \mcaption{The minimum distance, $d(\Gamma,\gamma)$, between the vesicle
 and the solid wall, and the time, $t_{\min}$, when this distance is
 achieved.  Rather than forming an overrefined ''true" solution, we
 report the actual distances.}{t:stenosis:gapSize}
 \end{centering}
\end{table}



\subsection{Couette Apparatus}
\label{s:couette}
Here we consider two simulations of a Couette apparatus with different
volume fractions.  We first consider a Couette apparatus where the
inner boundary, which is slightly off-centered, is rotating with
constant angular velocity and the outer boundary is fixed.  We randomly
place 42 vesicles of reduced area 0.65 inside the apparatus which
corresponds to a volume fraction of 20\%.  Each vesicle is discretized
with 64 points, each solid wall with 64 points, the time step size is
$0.01$, and the time horizon is $T=30$.  Snapshots of the simulation
are in the top plots of Figure~\ref{f:couetteLow}.  Using first-order
time stepping with explicit inter-vesicle interactions, the error in
area is $1.70e-2$ and the error in length is $7.94e-2$ at $t=30$.
Semi-implicit inter-vesicle interactions are not required for this time
step size.  The simulation took a little less than 6 hours, or about
7.1 seconds per time step.

We repeat this experiment with second-order time stepping which
requires semi-implicit inter-vesicle interactions.  The simulation took
about 83 hours, or about 100 seconds per time step.  The errors in area
and length at $t=30$ are $2.95e-3$ and $7.13e-4$.  To achieve a similar
accuracy with first-order time stepping, a time step size 11 times
smaller would be required which would require more than 83 hours of CPU
time due to the double-layer potentials.  Moreover, the errors in area
and length for second-order time stepping have plateaued while for
first-order time stepping, they continue to grow (bottom plots of
Figure~\ref{f:couetteLow}).  Again, this means that with second-order
time stepping, we can take much longer time horizons without requiring
a smaller time step size.

\begin{figure}[htps]
\begin{center}
  \begin{tabular}{c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}}
    \input{couetteSnapsTime1.tikz} &
    \input{couetteSnapsTime2.tikz} &
    \input{couetteSnapsTime3.tikz} &
    \input{couetteSnapsTime4.tikz} &
    \input{couetteSnapsTime5.tikz} &
    \input{couetteSnapsTime6.tikz} 
  \end{tabular}
  \begin{tabular}{cc}
    \input{couetteErrorsOrder1.tikz} & 
    \input{couetteErrorsOrder2.tikz} 
  \end{tabular}
\end{center}
\mcaption{The top plots show snapshots of a Couette flow with 42
vesicles and a volume fraction of 20\%.  One vesicle is coloured in blue
to view its time history.  The outer boundary is fixed while the inner
boundary completes one revolution every 10 time units.  The bottom plots
show the errors in area and length for the first- and second-order time
integrators.  Notice that the errors for first-order method  continue to
grow whereas the errors for the second-order method
plateau.}{f:couetteLow}
\end{figure}

We now consider a Couette apparatus with the inner boundary exactly in
the middle of the outer boundary.  We use 150 vesicles of reduced area
0.65 resulting in a volume fraction of approximately $40.5\%$.  With
this higher concentration, a finer discretization of the solid walls is
required to resolve its density function ${\boldsymbol\eta}$.  Therefore, we have
discretized both boundaries with 256 points.  We also used our
collision detection algorithm to alert us if vesicles had crossed.
With $\Delta t = 0.01$ and explicit inter-vesicle interactions, the
first-order method results in crossing vesicles.  Therefore, we report
the results using first-order semi-implicit inter-vesicle interactions
which does not result in vesicles crossing.  Snapshots are illustrated
in Figure~\ref{f:couetteHigh}.  On average, each time step required 27
GMRES iterations. The total simulation took about 10 days, but this
number will drop dramatically once the double-layer potential is
implemented with a fast summation method.  At $t=10$, the error in area
is $9.62e-3$ and the error in length is $1.92e-2$.  We also ran this
simulation with second-order time stepping and achieved an error in
area of $9.39e-4$ and an error in length of $4.19e-4$.
\begin{figure}[htps]
\begin{center}
  \begin{tabular}{c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}}
    \input{couetteDenseSnapsTime1.tikz} &
    \input{couetteDenseSnapsTime2.tikz} &
    \input{couetteDenseSnapsTime3.tikz} &
    \input{couetteDenseSnapsTime4.tikz} &
    \input{couetteDenseSnapsTime5.tikz} &
    \input{couetteDenseSnapsTime6.tikz} 
  \end{tabular}
\end{center}
\mcaption{Snapshots of a Couette flow with 150 vesicles.  One vesicle
is coloured in blue to view its time history.  The inner boundary
completes one full revolution every 10 time units.}{f:couetteHigh}
\end{figure}


