Journal of Computational Physics 211 (2006) 531–550

www.elsevier.com/locate/jcp

A Cartesian grid embedded boundary method for the

heat equation and Poisson
s equation in three dimensions
b,*,1,2 a,3 b,1,2 b,2
Peter Schwartz , Michael Barad , Phillip Colella , Terry Ligocki
a
Department of Civil and Environmental Engineering, University of California, Davis, CA 95616, United States
b
Applied Numerical Algorithms Group/CRD, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Mailstop 50-A 1148,
Berkeley, CA 94720, United States

Received 2 November 2004; received in revised form 2 June 2005; accepted 7 June 2005
Available online 3 August 2005

Abstract

We present an algorithm for solving Poisson
s equation and the heat equation on irregular domains in three dimen-
sions. Our work uses the Cartesian grid embedded boundary algorithm for 2D problems of Johansen and Colella [A
Cartesian grid embedded boundary method for Poisson
s equation on irregular domains, J. Comput. Phys. 147(2)
(1998) 60–85] and extends work of McCorquodale, Colella and Johansen [A Cartesian grid embedded boundary
method for the heat equation on irregular domains, J. Comput. Phys. 173 (2001) 620–635]. Our method is based on
a finite-volume discretization of the operator, on the control volumes formed by intersecting the Cartesian grid cells
with the domain, combined with a second-order accurate discretization of the fluxes. The resulting method provides
uniformly second-order accurate solutions and gradients and is amenable to geometric multigrid solvers.

2005 Elsevier Inc. All rights reserved.

PACS: 02.60.Lj; 02.70.Bf; 41.05.+e; 41.20.Cv

Keywords: Poisson equation; Heat equation; Multigrid methods

*
Corresponding author. Tel.: +1 510 486 7507.
E-mail address: poschwartz@lbl.gov (P. Schwartz).
1
Supported by the DARPA BioComp program.
2
Supported at the Lawrence Berkeley National Laboratory by the U.S Department of Energy: Director, Office of Science, Office of
Advanced Scientific Computing, Mathematical, Information, and Computing Sciences Division under Contract DE-AC03-76SF00098.
3
Supported by the Computational Science Graduate Fellowship program of the Department of Energy, under Grant No. DE-FG02-
97ER25308.

0021-9991/$ - see front matter
2005 Elsevier Inc. All rights reserved.
doi:10.1016/j.jcp.2005.06.010
532 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

1. Introduction

In this paper, we present Cartesian grid embedded boundary methods for solving Poisson
s equation and
the heat equation on irregular domains in three dimensions. In this approach, the irregular domain is dis-
cretized as a collection of control volumes formed by the intersection of the domain with cubic Cartesian
grid cells. The primary unknowns are defined at the centers of the Cartesian cells, while the Laplacian is
approximated by a finite-volume discretization on each of the regular or irregular control volumes. This
approach was successfully used in [7] to solve elliptic PDE in two dimensions and extended to parabolic
PDE in two dimensions in [12]. There have been a number of other investigators using similar approaches
working in two dimensions [17] or in the context of specific applications problems in three dimensions
[9,14]. In this paper, we present a systematic extension of the approach in [7] to three dimensions. The prin-
cipal new feature of the algorithm is the use of bilinear interpolation to approximate the fluxes on irregular
faces to obtain a stable and formally consistent discretization of the Laplacian. In addition, we introduce an
alternate lower-order stencil for computing the flux on a Dirichlet boundary at locations where the bound-
ary geometry is underresolved on the grid. The latter feature permits the use of multigrid coarsening of the
grid to much coarser levels than in [7]. The resulting method is second-order accurate in space for Poisson
s
equation, and second-order accurate in space and time for the heat equation.
Cartesian grid methods for elliptic PDE have a long history beginning with the nodal-point schemes in
Shortley and Weller [15]. More recently, these methods have been rediscovered and extended by a number
of authors [6,5,10,11]. The approach we
re taking differs in that it is based on finite-volume discretizations.
For hyperbolic PDE, such methods also have a long history [13] (For a recent survey see [4].) For nonlinear
hyperbolic problems, the local conservation property of such methods is essential to obtain correct weak
solutions. For elliptic and parabolic problems, there are several reasons why local conservation is desirable.
For elliptic problems that arise as part of the calculation of hyperbolic fluxes, the need for local conserva-
tion is driven by the requirement for the hyperbolic problem. More generally, there are many examples in
which the use of locally conservative methods for elliptic and parabolic equations arising in fluid dynamics
provides more accuracy and greater robustness in marginally resolved calculations. Finally, the use of
finite-volume methods for elliptic PDE on locally refined meshes leads to easily computable solvability con-
ditions for Neumann and periodic problems. Corresponding solvability conditions are not known for stan-
dard nodal-point methods on locally refined grids, such as those described in [1,11]. A disadvantage of the
finite-volume approach is that formally consistent discretizations do not lead to symmetric linear systems.
This is in contrast to the finite-element based approach to Cartesian grid methods in [18], for which the
variational formulation guarantees the symmetry of the resulting linear system. In the present case, the sta-
bility of finite-volume methods typically must be studied empirically.

2. Embedded boundary discretization of the Laplacian in 3D

The underlying discretization of space is given by rectangular control volumes on a Cartesian grid:

i ¼ ½ih; ði þ uÞh
; i 2 ZD , where D is the dimensionality of the problem, h is the mesh spacing, and u is
the vector whose entries are all ones. In the case of a fixed, irregular domain X, the geometry is represented
by the intersection of X with the Cartesian grid. We obtain control volumes Vi =
i \ X and faces Ai12ed ,
which are the intersection of oVi with the coordinate planes fx : xd ¼ ðid  12Þhg. Here ed is the unit vector
in the d-direction. We also define ABi to be the intersection of the boundary of the irregular domain with the
Cartesian control volume: ABi ¼ oX \
i . We will assume here that there is a one-to-one correspondence
between the control volumes and faces and the corresponding geometric entities on the underlying
Cartesian grid. The description can be generalized to allow for boundaries whose width is less than the
mesh spacing or boundaries with sharp trailing edges.
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 533

In order to construct finite difference methods, we will need only a small number of real-valued quanti-
ties that are derived from these geometric objects.

 Areas and volumes are expressed in dimensionless terms: volume fractions ji = |Vi|hD, face apertures
ai12ed ¼ jAi12ed jhðD1Þ and boundary apertures aBi ¼ jABi jhðD1Þ . We assume that we can compute esti-
mates of the dimensionless quantities that are accurate to O(h2).
 The locations of centroids and the average outward normal to the boundary are given exactly by:
Z
1
Face centroid : xiþ12ed ¼ x dA;
jAiþ12ed j A 1
iþ ed
2
Z
B 1
Boundary face centroid : xi ¼ B x dA;
jAi j ABi
Z
B 1
Outward normal : ni ¼ B nB dA;
jAi j ABi
where nB is the outward normal to oX, defined for each point on oX. Again, we assume that we can com-
pute estimates of these quantities that are accurate to O(h2).

Using just these quantities, we can define conservative discretizations for the divergence operator. Let
~
F ¼ ðF 1 ; . . . ; F D Þ be a function of x. Then
Z Z " ! #
1 1 1 X X D
r~ F  r~
F dV ¼ ~
F  n dA  ai12ed F d ðxi12ed Þ þ aBi nBi  ~
F ðxBi Þ ;
jV i j V i jV i j oV i ji h ¼þ; d¼1

ð1Þ
where (1) is obtained by replacing the integrals of the normal components of the vector field ~
F with the
values at the face centroids.
We first consider Poisson
s equation on an irregular domain X:
Dw ¼ q on X;
ow
¼ gN on oX;
on ð2Þ
or
w ¼ gD on oX.
We define a discrete variable /, /i  w(ih). Using the discretization of the divergence defined in (1), we can
define a discretization of Poisson
s equation as follows:
ðDh /Þi ¼ qi ; ð3Þ
" ! #
1 X X
D
ðDh /Þi ¼ ai12ed F di1ed þ aBi F B ; ð4Þ
ji h ¼þ; d¼1
2

where qi = q(ih), and the fluxes Fd and FB are linear combinations of /i and the boundary values. In prac-
tice, we avoid problems arising from arbitrarily small values of ji in the denominator by solving

ji ðDh /Þi ¼ ji qi . ð5Þ

The fluxes are given by bilinear interpolation of centered differences. Explicitly, bilinear interpolation of
fluxes can be written as an iteration of linear interpolation of fluxes in the two directions that are not
534 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

normal to the face. For example, given the face with outward normal e1, with centroid x, define the linearly
interpolated flux in the d (d 6¼ 1) direction by:
ð/iþe1  /i Þ ð/iþe1 ed  /ied Þ
F diþ1e1 ¼ g þ ð1  gÞ ;
2 h h
jx  ed j
g¼1 ; ð6Þ
h


þ x  ed > 0;
¼
 x  ed 6 0.
The bilinear interpolation of the flux for the face with normal e1 can be written:
F iþ12e1 ¼ gF diþ1e1 þ ð1  gÞF die 0 þ1e1 ;
2 d 2

jx  ed 0 j
g¼1 ; ð7Þ

h
þ x  ed 0 > 0;
¼
 x  ed 0 6 0;
where d 0 6¼ d, d 0 6¼ 1. See Fig. 1. We note that the particular choice of bilinear interpolation for computing
the fluxes is a nontrivial one for obtaining a stable method. In particular, we also tried using simple linear
interpolation based on three of the faces in Fig. 1, omitting the face offset along the diagonal. We found
that such a method is unstable for some configurations of adjacent small control volumes, in the sense that
point Jacobi fails to converge for any value of the relaxation parameter. No such instability was observed
for the bilinear scheme. For that reason, we have chosen to reduce order to a piecewise constant interpolant
of the fluxes if all four faces required for a bilinear interpolant are unavailable.

2.1. Boundary conditions

For Neumann boundary conditions the flux on the boundary is specified. For Dirichlet boundary con-
ditions further calculations are necessary. Our primary method, for use when the geometry is well resolved,
is a generalization of the methods described in [7]. Fig. 2 shows how this generalizes to 3D
 
o/ 1 d2 B d1
 ð/  /I1 Þ  ð/B  /I2 Þ . ð8Þ
on d 2  d 1 d 1 d2

e3 e2

e1

Fig. 1. 3D bilinear flux.

 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 535

φI2

φI1

φB
P2

P1

Fig. 2. Diagram of the second-order stencil for the gradient normal to the interface.

Here, we have used /B for the value of / on the boundary, which is given by the Dirichlet boundary con-
dition. Interpolation from cell centered values determines /I1 and /I2 at distance d1 and d2 away from the
boundary, respectively.
If all 18 cells are available in Fig. 2, we make an order O(h2) estimate $/ as follows. Depending on the
orientation of the normal, two planes are chosen, P1 and P2. Using biquadratic interpolation, two values
/I1 and /I2 are calculated, each requiring 9 values.
The gradient is then calculated by fitting a quadratic to the interpolated values and the value at the inter-
face. We chose the planes P1, P2 to be perpendicular to ed, where d is given by
fd : nBd P nB‘ ; ‘ ¼ 1; 2; 3g. ð9Þ
In cases where the requisite eighteen cells are not available, we employ a lower order stencil to estimate
the flux to O(h). In 3D, this lower order stencil contains at most eight points including the centroid of the
embedded boundary. These eight points are chosen as follows. We associate each one of the eight
possible configurations of plus or minus signs of the components of the normal vector with one of
the octants in the coordinate system with origin at the centroid. The stencil consists of the cell-centers
of the seven nearest cells in this octant, excluding the cell-center of the control volume containing the
boundary centroid itself.
From these points we create an overdetermined linear system to estimate $/ as follows:
A  r/ ¼ d/; ð10Þ
where
A ¼ ðdx1 ; dx2 ; . . . ; dx7 ÞT ;
d/ ¼ ðd/1 ; d/2 ; . . . ; d/7 ÞT ;
dxm ¼ xm  xBi ;
d/m ¼ /m  /B .
We determine $/ using least squares by solving the normal equations
AT A  r/ ¼ AT d/.
536 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

Provided that A contains three linearly independent rows, AT A is invertible. This is always the case pro-
vided the set {xm} contains all points of the form (i + ed)h, plus at least one other point. If that is not the
case, we set $/ ” 0. These methods lead to a condition number for Dh that is bounded independent of j,
and comparable to that of the uniform grid algorithm.

2.2. Truncation and solution error

We define the truncation error in the usual fashion: si ¼ qi  Dh /exact

i , where /exact
i ¼ wðihÞ. We then have
the following asymptotic error estimates for the truncation error. At regular cells
si ¼ Oðh2 Þ. ð11Þ
If i is an irregular cell, and the flux on the boundary is second-order accurate, as in (8), then
 
h
si ¼ O . ð12Þ
ji
If we use the flux computation given by (10), the flux on the boundary is first-order accurate, and we have
 
1
si ¼ O . ð13Þ
ji
We refer to methods that satisfy truncation error estimates of the form (12) on the irregular control vol-
umes as being formally consistent. We also define the solution error i ¼ /i  /exact
i .
There is one apparent problem with this truncation error estimate: it is only first-order accurate at the
boundary. Nonetheless, we observe robust second-order convergence of the solution in max norm. These
two facts can be reconciled using a modified equation analysis [7]. Both methods of calculating the gradient
for Dirichlet boundary conditions, (8) and (10), lead to a second-order solution error. This is because, for
Dirichlet boundary conditions, solution error is two orders of accuracy more than the truncation error on
the boundary. On the other hand, in the following section we show that it is necessary to use (8) to obtain
second-order accurate values for $/ at the boundary.

3. Discretizing the heat equation

We now consider the heat equation on a moving domain. w : R3 · [0,1] ! R is the unknown and
f : R3 · [0,1] ! R is the source term. On a time-dependent domain, we solve:
ow
¼ Dw þ f ;
ot
ow
¼ gN ðx; tÞ; x 2 oXðtÞ; ð14Þ
on
or
w ¼ gD ðx; tÞ; x 2 oXðtÞ.

3.1. Time discretization for a fixed domain

First, we consider the case of a fixed domain, i.e., X(t) = X independent of time. We define discrete vari-
ables, /i(t) and fi(t),
/i ðtÞ  wðih; tÞ; f i ðtÞ  f ðih; tÞ. ð15Þ
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 537

This leads to a semi-discrete system of ODEs

d/i
¼ Dh /i þ fi . ð16Þ
dt
We discretize this system in time using the L0-stable method [16], which was also described in [12], as out-
lined below.
Denote by I the identity operator, and by DhI and DhH the discrete Laplacian with inhomogeneous and
homogeneous boundary conditions, respectively. We split the timestep Dt such that:

l1 þ l2 þ l3 ¼ Dt;
l1 þ l2 þ l4 ¼ Dt=2.

The update at step n uses the boundary values at the old and new times and also at an intermediate time tint
1 1
/nþ1 ¼ ðI  l1 DhI ðtnew ÞÞ ðI  l2 DhI ðtint ÞÞ  ½ðI þ l3 DhI ðtold ÞÞ/n þ ðI þ l4 DhH Þf ðtavg ÞDt
; ð17Þ

where

(x new,t new)

(x new,t int)

x
(x old ,t old ) (x new,t old)

Fig. 3. Boundary conditions for the equivalent problem are extrapolated for (xnew,tint) and (xnew,told).

Fig. 4. Centers of cells in X(told) are shown with solid circles and centers of cells in X(tnew)  X(told) are shown with unfilled circles. To
estimate the value of /n at one of these latter points, we extrapolate quadratically from values at the centers of three other cells in
X(told) forming a line with the new cell center. We pick whichever such line is closest in direction to the normal at time tnew.
538 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

told ¼ nDt;
tnew ¼ ðn þ 1ÞDt ¼ told þ l1 þ l2 þ l3 ;
tint ¼ tnew  l1 ¼ told þ l2 þ l3 ;
tavg ¼ ðtold þ tnew Þ=2 ¼ told þ l1 þ l2 þ l4 .

For a second-order L0-stable method, following [16], we pick a > 1/2 and:
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
a a2  4a þ 2
l1 ¼ Dt;
2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
a þ a2  4a þ 2
l2 ¼ Dt;
2
l3 ¼ ð1  aÞDt;
 
1
l4 ¼  a Dt.
2
pffiffiffi
For a method that uses real arithmetic only, the truncation error is minimized by taking a ¼ 2  2  ,
where  is machine precision.
It was shown in [12] for Dirichlet boundary conditions that Crank–Nicolson time discretization exhib-
ited oscillatory behavior and furthermore was unstable to some types of forcing at a moving boundary. In

2
10

1
10

0
10
Truncation Error

–1
10

–2
10

–3
10
16 32 64 128 256
Grid Size

Fig. 5. Truncation error for Neumann boundary conditions. The symbol h denotes the L1 norm. The denotes L2, and the e denotes
*
the max norm. The reference line for first-order is given by – –, and the reference line for second-order is shown with Æ – Æ.
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 539

[8] this behavior was attributed to the combination of the neutral stability of Crank–Nicolson at high wave
numbers and the presence of eigenvalues of Dh with nontrivial imaginary parts, corresponding to eigen-
values with oscillatory behavior near the boundary.

3.2. Time discretization for a moving domain

In the moving case, the domain X is now a function of time, X = X(t), and the various geometric quan-
tities can also be computed in a time-dependent way: ji(t), aiþ12ed ðtÞ, and so on.
The timestep is assumed to satisfy a CFL condition with respect to the normal velocity v ¼ dx
dt
n

Dt
jvj < 1.
h
In the present approach, we solve the moving-boundary problem by defining an equivalent fixed-boundary
problem for each timestep. Specifically, we solve at each time step the discretization of the following fixed-
boundary problem:

owfixed
ðx; tÞ ¼ DDwfixed ðx; tÞ þ f ðx; tÞ; ð18Þ
ot

2
10

1
10

0
10
Truncation Error

–1
10

–2
10

–3
10
16 32 64 128 256
Grid Size

Fig. 6. Truncation error for Dirichlet boundary conditions. The symbol h denotes the L1 norm. The denotes L2, and the e denotes
*
the max norm. The reference line for first-order is given by – –, and the reference line for second-order is shown with Æ – Æ.
540 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

where
wfixed ¼ wfixed ðx; tÞ; x 2 Xðtnew Þ; told 6 t 6 tnew ;
fixed
ow
¼ gextrap
N ðx; tÞ; x 2 oXðtnew Þ;
on ð19Þ
or
wfixed ¼ gextrap
D ðx; tÞ; x 2 oXðtnew Þ.

The boundary conditions gextrap

N ;D on the fixed boundary are computed by extrapolating values from the mov-
ing boundary to the points on the fixed boundary oX(tnew) at times told and tint, as in (20). Fig. 3 shows a
one-dimensional schematic with time and Fig. 4 shows the boundary at two different times.
The steps required in setting up the fixed-boundary problem (19) are:

1. Extend the domain of /n to X(tnew) and define the newly uncovered values by extrapolation.
2. Compute boundary values for /n + 1 at ðxBi ðtnew Þ; told Þ and ðxBi ðtnew Þ; tint Þ.

In Step 1, to estimate the value of /n at the center xi(tnew) of a newly uncovered cell in X(tnew)  X(told),
we use a quadratic extrapolant from three other cells in X(told), such that the centers of these cells form a
line with xi(tnew). We choose whichever line passing through the centers of the new cell and one of its imme-
diate neighbors has a direction closest to that of the normal nBi ðtnew Þ. See Fig. 4. This differs from the
approach used in [5] in that we do not extrapolate in the direction of the normal.

–2
10

–3
10

–4
10
Gradient Error

–5
10

–6
10

–7
10
64 128 256
Grid Size

Fig. 7. Convergence of the first component of the gradient of the solution to Poisson
s equation with Dirichlet boundary conditions,
using (8) to compute the flux. Notation as in Fig. 6.
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 541

–1
10

–2
10

–3
10
Gradient Error

–4
10

–5
10

–6
10
64 128 256
Grid Size

Fig. 8. Convergence of the first component of the gradient of the solution to Poisson
s equation conditions, using (10) to compute the
flux on the boundary. Notation as in Fig. 6.

Fig. 9. Solution error for Poisson
s problem on a 643 grid.

542 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

Fig. 10. Solution error for Poisson
s problem on a 643 grid.

2
10

0
10

–2
10

–4
10
Residual

–6
10

–8
10

–10
10

–12
10
0 5 10 15 20 25
V Cycle Iteration

Fig. 11. Plot of the max norm of the partial volume-weighted residual versus V-cycle iteration. The graphs are for meshes of size 64
(h), 128 ( ), and 256 (Æ – Æ).
*
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 543

In Step 2, for each Vi 2 X(tnew)  X(told), we choose the j in the intersection of the 3 · 3 · 3 block of cells
with center i and X(told) that has the largest boundary area. For our extrapolation, we define

dði; jÞ ¼ ðxBi ðtnew Þ; told Þ  ðxBj ðtold Þ; told Þ.

We also need approximations to $w and $$w, which we estimate as follows.

~i ¼ rwðxi ðtnew Þ; told Þ þ OðhÞ. Each component of Gd is computed separately by differentiating the
Let G i
quadratic interpolant through three points chosen from /i, /ied , and /i2ed , where the sign of ± is chosen
so that all points are in X(tnew) and therefore where /n has been computed.
Explicitly, for the d component,
 
d 1 3 n n 1 n
Gi ¼   /i þ 2/ied  /i2ed
h 2 2

or
1
Gdi ¼ ð/niþed  /nied Þ
h
depending on whether i is on the end or in the middle of 3 cells.
We estimate second derivatives as follows. Let GGi  $$w(xi(tnew),told). The order of the error in this
approximation will vary between one and two, depending on the local geometry. For derivatives of the
form /d,d, we differentiate the quadratic extrapolant through /i, /ied , and /i2ed according to the stencil

2
10

0
10

–2
10

–4
10
Residual

–6
10

–8
10

–10
10

–12
10
0 2 4 6 8 10 12 14
W Cycle Iteration

Fig. 12. Plot of the max norm of the partial volume-weighted residual versus W-cycle iteration. The graphs are for meshes of size 64
(h), 128 ( ), and 256 (Æ – Æ).
*
544 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

1 n
Gd Gdi ¼ ð/iþed  2/ni þ /nied Þ.
h2
This term is second-order if i is in the center of the three point stencil and first-order otherwise. For mixed
derivatives, (/d;d 0 with d 6¼ d 0 ), we average one, two, three, or four estimates of the derivative each of which
is second-order accurate at points of the form i ± ed ± ed 0 . For example, assuming that i + ed, i + ed 0 , and
i + ed + ed 0 are all in X(tnew), we include the following in our average:
0 1 n
Gd Gdi ¼ ð/iþed þed 0  /niþed  /niþed 0 þ /niþed Þ.
h2
Given these estimates of $w and the matrix $$w at cell centers, we extrapolate the boundary conditions
using:
gextrap
N
~ i þ nB  ½ G
¼ gN ðxBi ðtold Þ; told Þ þ ðnBj  nBi Þ  G i
~G~i dði; jÞ

or ð20Þ
gextrap
D ¼ gD ðxBi ðtold Þ; told Þ ~i  dði; jÞ;
þG
where as remarked above, j is chosen from among the neighbors of i as the neighbor with the largest bound-
ary area.
Finally, we linearly interpolate the boundary conditions at ðxBi ðtnew Þ; tint Þ from the boundary conditions
at ðxBi ðtnew Þ; told Þ and ðxBi ðtnew Þ; tnew Þ. The former is calculated and the latter is given as part of the data of
the problem.

0
10

–1
10
Solution Error

–2
10

–3
10

64 128 256
Grid Size

Fig. 13. Solution error for Poisson
s equation with Dirichlet boundary conditions. Notation as in Fig. 6.
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 545

4. Numerical results

For our test problems, we compute the max norm of the solution error and truncation error. For a dis-
crete variable, n, the max norm is given by
knk1 ¼ max jni j
i

and the p-norm in 3D is given by

!1p
X p D
knkp ¼ jðni Þ h ji j .
i

4.1. Truncation error for the Laplacian and solution error for Poisson
s equation

We estimate the truncation error of the discretized Laplacian using the test function
f ðx; y; zÞ ¼ sinðxÞ sinð2yÞ sinð3zÞ. ð21Þ

Here, X is a sphere of radius r = 0.392 centered at the origin. Our discretized Laplacian has inhomogeneous
boundary conditions of either Neumann or Dirichlet type. Figs. 5 and 6 show that the discretization of the
operator has the accuracy anticipated by (11)–(13). If we set w = f, w satisfies:

0
10

–1
10
Solution Error

–2
10

–3
10

64 128 256
Grid Size

Fig. 14. Solution error for Poisson
s equation with Neumann boundary conditions. Notation as in Fig. 6.
546 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

Dw ¼ 14f on X;
ow of
¼ on oX;
on on ð22Þ
or
w¼f on oX.
Figs. 13 and 14 show the solution error for Poisson
s equation with Dirichlet boundary conditions (using the
higher order stencil) and Neumann boundary conditions. Figs. 7 and 8 show one of the differences between
the higher and lower order Dirichlet stencils. If the lower order stencil is used, $/  $/exact = O(h), as ex-
pected. However, for the higher order stencil $/  $/exact = O(h2). In other words, our numerical experi-
ments imply that for these cases, at least, our estimated solution has the asymptotic form / = /exact +
Ch2, where C is a smooth function for the higher order stencil, but not for the lower order stencil. In Figs.
9 and 10, we display the smoothness of the solution error.
In Figs. 11 and 12, we show the effectiveness of the multigrid V cycles and W cycles. We considered
mesh sizes of 16, 32, 64, 128, and 256 (only the last three are shown). Using the test problem (22), we
reduced the residual eleven orders of magnitude by using V cycles and by using W cycles. Fig. 11 shows
that for each doubling of the mesh size, four or five more V cycles are required to achieve the same
reduction of the residual. By contrast, increasing the mesh size does not increase the number of W
cycles required. On the other hand, W cycles require more work than V cycles and at these resolutions
solving this problem with V cycles takes less computational time than solving it with W cycles (see
Figs. 13 and 14).

–3
10

–4
10
Solution Error

–5
10

–6
10
64 128 256
Grid Size

Fig. 15. Solution error for the heat equation on the Neutrophil geometry with Neumann boundary conditions. Notation as in Fig. 6.
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 547

4.2. Solution error for the heat equation

Our test problems for the heat equation (14) have as their solution
 2 2 2
þy þz
4 exp  x 5ðtþ1Þ
wðx; y; z; tÞ ¼ ð23Þ
5pðt þ 1Þ
satisfying
wt ¼ Dw þ f ; ð24Þ
where
 
4ðx2 þ y 2 þ z2 þ 5ðt þ 1ÞÞ x2 þ y 2 þ z2
f ðx; y; z; tÞ ¼ 3
exp  .
125pðt þ 1Þ 5ðt þ 1Þ
In the case where the boundary is not moving we solve the Neumann problem on a computational domain
based on a neutrophil. Beginning with slice-data generated by confocal microscopy, we constructed an im-
plicit function with the following property: the surface of the neutrophil is implicitly defined by the set of
points at which the implicit function takes the value zero. From this implicit definition, we construct the
volume fractions, apertures, normals, and centroids necessary for Cartesian grid embedded boundary
methods.
We display the solution error in Fig. 18, the solution in Fig. 19, and a convergence study in Fig. 15.

–4
10

–5
10

–6
10
Solution Error

–7
10

–8
10

–9
10
16 32 64 128 256
Grid Size

Fig. 16. Solution error for the heat equation on a moving domain with Neumann boundary conditions. Notation as in Fig. 6.
548 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

–4
10

–5
10

–6
10
Solution Error

–7
10

–8
10

–9
10
16 32 64 128 256
Grid Size

Fig. 17. Solution error for the heat equation on a moving domain with Dirichlet boundary conditions. Notation as in Fig. 6.

Fig. 18. Solution error for the heat equation for the 2563 Neutrophil Geometry.
 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550 549

Fig. 19. Solution for the heat equation for the 2563 Neutrophil Geometry.

In the case of a moving boundary, we numerically solve the Neumann and Dirichlet problems on a
spherical domain, with boundary conditions of Neumann or Dirichlet type computed using the exact
solution.
The radius of the sphere changes with time, increasing at a prescribed speed. We advance the solution in
time from t = 0 to t = 0.1875 using a mesh spacing h and corresponding timestep Dt such that Dt/h = 0.5.
The number of timesteps equals 0.1875/Dt and h = 2n, where n = 4, . . . ,8. The solution error is shown in
Figs. 16–19.

5. Conclusion

We have described a Cartesian grid embedded boundary algorithm for solving Poisson
s equation and
the heat equation on irregular domains in three dimensions. The resulting method provides uniformly sec-
ond-order accurate solutions and gradients and is amenable to geometric multigrid solvers. In the future,
we plan on extending this approach in a variety of directions: to generalize this approach to free boundary
value problems, in which the discretizations used here are used on both sides of the boundary, combined
with jump relations at the boundary, to define a finite-volume discretization; and to combine this approach
for elliptic and parabolic problems with the approach described in [3] to solve a variety of problems in low-
Mach number flows with fixed and free boundaries. To carry this program out, it will be useful to extend
the approach described here to the case of adaptively refined meshes, for which the only outstanding issue is
the question of discretizing the operator when the embedded boundary crosses a boundary between refine-
ment levels. Finally, it would be interesting to attempt to extend the method described here to higher order,
e.g. fourth-order accuracy. A starting point for this could be the finite-volume formulation of fourth-order
Mehrstellen methods in [2].
550 P. Schwartz et al. / Journal of Computational Physics 211 (2006) 531–550

Acknowledgments

We thank Adam Arkin, Matt Onsum, and David Adalsteinsson for help with generating the neutrophil
geometry. Dan Graves, Dan Martin, Ted Sternberg, and Brian Van Straalen contributed helpful ideas as
well as information about the Chombo software library.

References

[1] Ann S. Almgren, Thomas Buttke, Phillip Colella, A fast adaptive vortex method in three dimensions, J. Comput. Phys. 113 (1994)
177–200.
[2] M. Barad, P. Colella, A fourth order accurate local refinement method for Poisson
s equation, J. Comput. Phys. 209 (1) (2005) 1–
18.
[3] P. Colella, D.T. Graves, B.J. Keen, D. Modiano, A Cartesian grid embedded boundary method for hyperbolic conservation laws,
J. Comput. Phys., in press.
[4] P. Colella, Volume-of-fluid methods for partial differential equations, in: E.F. Toro (Ed.), Godunov Methods: Theory and
Applications, Kluwer Academic/Plenum Publishers, New York, 2001, p. 3.
[5] Frederic Gibou, Ronald Fedkiw, A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains
with applications to the Stefan problem, J. Comput. Phys. 202 (2005) 577–601.
[6] Frederic Gibou, Ronald Fedkiw, Li-Tien Cheng, Myungjoo Kang, A second-order-accurate symmetric discretization of the
Poisson equation on irregular domains, J. Comput. Phys. 176 (2002) 205–227.
[7] H. Johansen, P. Colella, A Cartesian grid embedded boundary method for Poisson
s equation on irregular domains, J. Comput.
Phys. 147 (2) (1998) 60–85, December.
[8] Hans Svend Johansen, Cartesian Grid Embedded Boundary Finite Difference Methods for Elliptic and Parabolic Partial
Differential Equations on Irregular Domains. Ph.D. Thesis, University of California, Berkeley, 1997.
[9] M.P. Kirkpatrick, S.W. Armfield, J.H. Kent, A representation of curved boundaries for the solution of the Navier–Stokes
equations on a staggered three-dimensional Cartesian grid, J. Comput. Phys. 184 (2003) 1–36.
[10] Z. Li, A fast iterative method for elliptic interface problems, SIAM J. Numer. Anal. 35 (1998) 230.
[11] Peter McCorquodale, Phillip Colella, David P. Groteb, Jean-Luc Vay, A node-centered local refinement algorithm for Poisson
s
equation in complex geometries, J. Comput. Phys. 201 (2004) 34–60.
[12] P. McCorquodale, P. Colella, H. Johansen, A Cartesian grid embedded boundary method for the heat equation on irregular
domains, J. Comput. Phys. 173 (2001) 620–635.
[13] W.F. Noh, Cel: a time-dependent, two-space-dimension, coupled Eulerian–Lagrange code, SIAM J. Numer. Anal. 3 (1964).
[14] Stephane Popinet, Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries, J. Comput.
Phys. 190 (2003) 572–600.
[15] G.H. Shortley, R. Weller, The numerical solution of Laplace
s equation, J. Appl. Phys. 9 (1938) 334–344.
[16] E.H. Twizell, A.B. Gumel, M.A. Arigu, Second-order, l0-stable methods for the heat equation with time-dependent boundary
conditions, Adv. Comput. Math. 6 (1996) 333–352.
[17] H.S. Udaykumar, R. Mittal, P. Rampunggoon, A. Khanna, A sharp interface Cartesian grid method for simulating flows with
complex moving boundaries, J. Comput. Phys. 174 (2001) 345–380.
[18] David P. Young, Robin G. Melvin, Michael B. Bieterman, Forrester T. Johnson, Satish S. Samant, John E. Bussoletti, A locally
refined rectangular grid finite element method: application to computational fluid dynamics and computational physics, J.
Comput. Phys. 92 (1991) 1–66.