# Segmentation of Roads in Satellite Images
using specially modified U-Net CNNs

Jonas Bokstaller, Yihang She, Zhehan Fu, Tommaso Macrì
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The image classification problem has been deeply investigated by the research
community, with computer vision algorithms and with the help of Neural
Networks. The aim of this paper is to build an image classifier for satellite
images of urban scenes that identifies the portions of the images in which a
road is located, separating these portions from the rest. Unlike conventional
computer vision algorithms, convolutional neural networks (CNNs) provide
accurate and reliable results on this task. Our novel approach uses a sliding
window to extract patches out of the whole image, data augmentation for
generating more training/testing data and lastly a series of specially
modified U-Net CNNs. This proposed technique outperforms all other baselines
tested in terms of mean F-score metric.

##  I Introduction

The task of image segmentation consists of choosing the right label for each
pixel such that the image is segmented into different parts like in this work
a satellite image into road segments and background (e.g buildings, trees,
grass). This task is more important then ever because of the recent
developments of autonomous driving vehicles, which requires precise guidance
for save driving.

Advances in computational infrastructure and breakthroughs in artificial
neural network architectures (ANNs) are emerging. These lead to the use of
machine learning to tackle this partitioning task. In order to process the
amount of information in a single image efficiently, we need to use
convolutional neural networks (CNNs) [1] because they split, downscale and
upscale the image while applying different filters in an efficient way. This
enables the use of such networks for real time image segmentation with high
accuracy.

As training data we have satellite images as well as the corresponding ground-
truth map which labels each pixel as road or background. The task of the paper
is to predict the label of new satellite images in chunks of 16x16 pixels. The
test images have dimensions 608x608 pixels, therefore we have to predict 1444
chunks per image. In order to compare the performance of different models and
approaches, mean $F_{1}$ score is used. The road label is assigned to a chunk
if the proportion of road pixels exceed the threshold of 25% in the chunk.

Figure 1: Example from the provided training data. On the left there is the
input satellite image and on the right its respective ground-truth.

This paper gives the reader a brief overview over the different methods we
implemented, the difficulties we had encountered during testing and our final
solution which consists of a set of specially modified U-Net CNNs.

##  II Models and Methods

###  II-A General Workflow

Given the training set of 100 images and corresponding masks with size
$400\times 400$, the task requires classifying the $16\times 16$ patches of
the test set into either road or background. To conduct such classification,
we have tested three baseline models:

  1. 1.

The sliding window approach, which conducts a binary classification for each
patch with the given window size using convolutional neural networks.

  2. 2.

The segmentation approach based on U-Net to generate pixel-level
classification.

  3. 3.

Transfer learning strategy for segmentation. The encoder of U-Net is adapted
from the base model of MobileNetV2. Pre-trained parameters on ImageNet were
taken for initialization.

Besides, several data preprocessing techniques and training strategies were
applied. Our final solution involves averaging outputs of three models based
on U-Net, which result from improvements on network architecture of original
U-net model. This combination has significantly improved F-score on both
validation and test set.

###  II-B Data Pre-processing and Augmentation

As the number of aerial images in training set is relatively small (with only
100 images and corresponding masks), we generated 900 patches of $256\times
256$ in size from them. These patches are used throughout the second and third
baseline solution and the final one. Furthermore, image data augmentation
techniques were used to overcome the impacts of limited training data size. We
augmented the images by rotation ($rotation\\_range=0.2$), shifting the width
and height of the image center ($width\\_shift\\_range=0.2$,
$height\\_shift\\_range=0.2$), shearing and zooming the image
($shear\\_range=0.2$, $zoom\\_range=0.2$), and random flipping in both
horizontal and vertical directions. The corresponding masks were also
transformed with the same operations.

###  II-C Network Architecture

####  II-C1 Baseline solution 1: sliding-window approach

This approach is based on the work of Gouk et al.[2] The main idea behind this
approach is to classify each patch with padding region around it to provide
the model with enough context information. Then we use this window and scan
over the image in order to classify each patch. The patch plus the padding is
then fed into a convolutional neural network with 5 layers in total. To avoid
overfitting, Leaky ReLu was used as activation function, and Dropout was added
with a probability of 0.5.[3]

Figure 2: Patch with padding around it called sliding window approach. The
padding is used to provide the model with enough contextual information. The
stride of the sliding window is its own size, namely 16x16.

To further avoid overfitting and to increase the number of windows as training
data, we used a data augmentation strategy different from the aforementioned
one. First, a patch is randomly selected as training instance and the window
inside that patch is chosen. Then it underwent random flipping and rotation so
that the model’s generalization ability is further enhanced. If the patch is
at the border, the padding is just a mirrored image around its edge.

####  II-C2 Baseline solution 2: U-Net

We implemented the original U-Net architecture[4] on this task. (see Figure 3)
This is a small and easily trainable fully convolutional neural network
designed originally for binary segmentation of biomedical images. The network
is consisted of two paths, the encoding and decoding paths. In the first half
or the encoding path, operations are 2D convolutions followed by a non-linear
activation function, ELU.[5] Feature map size is retained in each convolution
by padding. Then max pooling is used to reduce the feature map in the encoding
path. After these operations, number of feature channels is doubled. Then the
U-Net uses an expansion path (decoding path) to up-sample the feature maps and
eventually create a segmentation mask same in size as input image. This path
consists of sequences of up sampling in the form of transpose convolutions,
followed by concatenations with feature maps copied from the encoding path.
Finally, a convolution of kernel size (1,1) and a Softmax function is applied.
This maps each entry in input aerial image to one of two classes, road or non-
road.

Figure 3: The structure of the modified U-net architecture.

####  II-C3 Baseline solution 3: U-Net with MobileNetV2 trained on ImageNet

In order to extract robust features, We further applied transfer learning,
specifically using pre-trained parts of a lightweight model, MobileNetV2[6] as
encoder path of U-Net. This is especially applicable in this task because it
has a rather limited training set, and the method does not require much
computational resource. The base model to transfer, MobileNetV2 is trained on
a large dataset of another task, in this case the ImageNet datset[7], so that
the weights of its convolutional layers (especially the first few ones) are
more robust in extracting features that can be applied to the road extraction
task. The fundamental architecture of the baseline model is based on U-Net. We
adapted its encoder from the convolutional layers of pre-trained MobileNetV2.
The output layer of the encoder is the layer $block\\_13\\_expand\\_relu$ of
MobileNetV2. Four output layers of MobileNetV2 ($input\\_image$,
$block1\\_expand\\_relu$, $block\\_3\\_expand\\_relu$,
$block\\_6\\_expand\\_relu$) were concatenated to the transposed convolutional
layers of the encoder to build skip connections.

####  II-C4 Final solution: improvements on U-Net

We made two major improvements on U-Net architecture and implemented three
different models (Figure 4). The final submitted result was generated through
averaging the values of masks generated by three models (after Softmax
function) and classification was done with a threshold of 0.25.

Figure 4: The major improvements on U-Net architecture: (a) addition of
identity mapping before final nonlinear activation in each convolution block,
as in ResNet; (b) Substituting the original bottleneck with dilated
convolutions.

First, we added local skip connections (similar to that of ResNet[8]) to the
convolution blocks in both the encoding and decoding pathways of U-Net. This
is done through addition of the current layer input’s identity mapping to its
output values. Each modified convolution block takes the following general
form:

$\begin{array}[]{l}\mathbf{y}_{l}=h\left(\mathbf{x}_{l}\right)+\mathcal{F}\left(\mathbf{x}_{l},\mathcal{W}_{l}\right)\\\
\mathbf{x}_{l+1}=f\left(\mathbf{y}_{l}\right)\end{array}$

where $\mathbf{x}_{l}$ and $\mathbf{x}_{l+1}$ are $l$-level’s input and output
tensors, $\mathcal{F}(\cdot)$ is the residual function, which is the original
convolutions followed by batch normalization and nonlinear activation.
$f\left(\mathbf{y}_{l}\right)$ is the activation function, and
$h\left(\mathbf{x}_{l}\right)$ is the identity mapping function, which is a
simple 2D convolution with kernel size 1 in our solution. This added operation
addresses the gradient degradation problem by letting it backpropagate through
this additional identity mapping function. It also facilitates training, and
utilizes features of lower semantic information extracted in the previous
layers.

The second improvement is substitution of the central bottleneck in U-Net with
a series of dilated convolution operations. The summation of these
convolutions are fed into the decoder path. Previous work [9], [10] reports
that adding dilated convolution modules can solve the problem of degraded
picture resolution. This is achieved through expanded receptive fields, thus
maintaining per-pixel level classification accuracy while generating large-
scale feature maps with rich context information.

In all, we implemented three models based on the above improvements. They are
summarized in the table below.

TABLE I: Models used in final solution Model name | Description
---|---
unet-32 | 32 channels on the first layer, local skip connections
unet-64 | 64 channels on the first layer, local skip connections
unet-dilated | parallel dilated convolutions as bottleneck, local skip connections

###  II-D Loss function and evaluation metrics

The models were trained to optimize the smoothed dice loss, defined as
follows:

$L=1-\frac{2|y\cap\hat{y}|+\epsilon}{|y|+|\hat{y}|+\epsilon}$

where $y$ is the mask values predicted by model, $\hat{y}$ is the ground-
truth, and $\epsilon$ is the smoothing coefficient, selected as 1. Model
performance is evaluated based on the following metrics: IoU (intersection
over union) and F-1 score ($t_{p}$ is true positive, $f_{p}$ and $f_{n}$ are
false positive and false negative ratios):

$IoU=\frac{|y\cap\hat{y}|}{|y\cup\hat{y}|}$
$F_{1}=\frac{\mathrm{t_{p}}}{\mathrm{t_{p}}+\frac{1}{2}(\mathrm{f_{p}}+\mathrm{f_{n}})}$

###  II-E Training

We split the data into training ($80\%$) and validation set ($20\%$). The
initial learning rate is 0.0001 and it will be reduced with a factor of 0.5 if
the validation loss is less than 0.0002 after every 5 epochs. The model was
built in Keras[11] and trained with a batch size of 8 with 100 epochs in
maximum. When training the unet-32 and unet-64 models (described in Table I),
dropout was added after each concatenation with probability 0.2. Batch
normalization was performed to accelerate training (as illustrated in Figure
4), in the position after each 2D convolution and before nonlinear activation
functions. Baseline solutions were trained using Google Colaboratory[12] and
final solution models were trained on a cluster using one NVIDIA®V100 GPU and
28 Intel®Xeon®E5-2690 v4 CPUs.

##  III Results

###  III-A Key results from baselines

####  III-A1 Sliding window approach results

Hyperparameters were fine-tuned for this approach and mean F-score on test set
was around 0.85.

####  III-A2 Transfer learning results

The IoU on the validation set finally achieves 0.7991, which cannot meet
typical standards of road extraction. Given that the outcome of transfer
learning also depends on the similarity between the pretrained dataset and the
specific task, the trained model conditioned on ImageNet might only be able to
find a suboptimal parameters for road segmentation.

###  III-B Final solution results

We observed a certain degree of improvement on all three models proposed above
than the traditional U-Net architecture. Detailed performance metrics on
validation and test set are listed below:

TABLE II: Performance metrics of proposed models Model | val_loss | val_f1 | test_f1
---|---|---|---
unet-32 | 0.0826 | 0.9174 | 0.8780
unet-64 | 0.0368 | 0.9632 | 0.8959
unet-dilated | 0.0334 | 0.9666 | 0.9015

Averaging the values of masks predicted by these models yield higher $F_{1}$
score on test set (0.9027), which was chosen for submission.

##  IV Discussion

In this work we first discovered that by using encoder-decoder architecture
such as U-net, semantic information in different levels are effectively
coupled to give clearer and more accurate results than a classification of
each patch with a simple 5-layer CNN. By further substituting bottleneck of
U-Net with parallel dilated convolution structure and adding residual-like
skip connections, we largely exceeded performance of baseline solutions.
However, there is still room of improvement for each of the models proposed.
The gap between score on test set and validation set suggests that the problem
of lacking training data is still not resolved. Possible candidate methods are
image sharpening with filters and color jittering. Also, different selections
of layer parameters transferred into U-Net as well as different tasks the base
model is trained on needs to be explored further. This can affect overall
model performance and the amount of time required to train until convergence.

## References

  * [1] Y. Lecun, P. Haffner, and Y. Bengio, “Object recognition with gradient-based learning,” 08 2000.
  * [2] H. Gouk and A. Blake, “Fast sliding window classification with convolutional neural networks,” 11 2014, pp. 114–118.
  * [3] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: a simple way to prevent neural networks from overfitting,” _The journal of machine learning research_ , vol. 15, no. 1, pp. 1929–1958, 2014.
  * [4] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [5] D.-A. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (elus),” _arXiv preprint arXiv:1511.07289_ , 2015.
  * [6] M. Sandler, A. Howard, M. Zhu, A. Zhmoginov, and L.-C. Chen, “Mobilenetv2: Inverted residuals and linear bottlenecks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018, pp. 4510–4520.
  * [7] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.
  * [8] K. He, X. Zhang, S. Ren, and J. Sun, “Identity mappings in deep residual networks,” in _European conference on computer vision_. Springer, 2016, pp. 630–645.
  * [9] S. Piao and J. Liu, “Accuracy improvement of unet based on dilated convolution,” in _Journal of Physics: Conference Series_ , vol. 1345, no. 5. IOP Publishing, 2019, p. 052066\.
  * [10] F. Yu and V. Koltun, “Multi-scale context aggregation by dilated convolutions,” _arXiv preprint arXiv:1511.07122_ , 2015.
  * [11] F. Chollet _et al._ (2015) Keras. [Online]. Available: https://github.com/fchollet/keras
  * [12] E. Bisong, “Google colaboratory,” in _Building Machine Learning and Deep Learning Models on Google Cloud Platform_. Springer, 2019, pp. 59–64.


this dataset comprises of a simple logistic regression as the main model, with
only the scalar logit outputs of the logistic regression being conditioned
using FiLM. The MLPs $M_{\mu},M_{\sigma}$ both have two hidden layers, each of
size 4, and ReLU activations. We train the model for a maximum of 1000 epochs,
with early stopping based on validation losses. Training these simple model
only requires one CPU and takes roughly 2 minutes.

Ground truth trade-off curve To obtain the ground truth trade-off curve
$\tau^{\ast}_{\textup{fair}}$, we consider the family of classifiers

$\displaystyle h_{c}(X)=\mathbbm{1}(X>c)$

for $c\in\mathbb{R}$. Next, we calculate the trade-offs achieved by this model
family for a fine grid of $c$ values between -3 and 3, using a dataset of size
500,000 obtained using the data-generating mechanism described above. The
large dataset size ensures that the finite sample errors in accuracy and
fairness violation values are negligible. This allows us to reliably plot the
trade-off curve $\tau^{\ast}_{\textup{fair}}$.

####  E.6.1 Results

Figure 27 shows the results for the synthetic data setup for three different
fairness violations, obtained using a calibration dataset
$\mathcal{D}_{\textup{cal}}$ of size 2000. It can be seen that for each
fairness violation considered, the YOTO trade-off curve aligns very well with
the ground-truth trade-off curve $\tau^{\ast}_{\textup{fair}}$. Additionally,
we also consider four different confidence intervals obtained using our
methodology, and Figure 27 shows that each of the four confidence intervals
considered contain the ground-truth trade-off curve. This empirically verifies
the validity of our confidence intervals in this synthetic setting.

Additionally, in Figure 28 we plot how the worst-case values for $\Delta(h)$
change (relative to the optimal trade-off $\tau^{\ast}_{\textup{fair}}$) as
the number of training data $\mathcal{D}_{\textup{tr}}$ increases. Here, we
use $h_{\lambda}$ as a short-hand notation for the YOTO model
$h_{\theta}(\cdot,\lambda)$ and the quantity on the $y$-axis,
$\max_{\lambda\in\Lambda}\frac{\Delta(h_{\lambda})}{\tau^{\ast}_{\textup{fair}}(\textup{acc}(h_{\lambda})},$
can be intuitively thought of as the worst-case error percentage between YOTO
trade-off and the optimal trade-off. The figure empirically verifies our
result in Theorem 3.4 by showing that as the number of training data
increases, the error term declines across all fairness metrics.

Figure 27: Results for the synthetic dataset with ground truth trade-off
curves $\tau^{\ast}_{\textup{fair}}$. Figure 28: Plot showing how $\Delta(h)$
decreases (relative to the ground truth trade-off value
$\tau^{\ast}_{\textup{fair}}(\textup{acc}(h))$) as the training data size
$|\mathcal{D}_{\textup{tr}}|$ increases. Here, we plot the worst (i.e.
largest) value of
$\Delta(h_{\lambda})/\tau^{\ast}_{\textup{fair}}(\textup{acc}(h_{\lambda}))$
achieved by our YOTO model over a grid of $\lambda$ values in $[0,5]$.


We present a new discretization approach to advection-diffusion problems with Robin boundary conditions on complex, time-dependent domains. The method is based on second order cut cell finite volume methods introduced by Bochkov  <cit.> to discretize the Laplace operator and Robin boundary condition. To overcome the small cell problem, we use a splitting scheme along with a semi-Lagrangian method to treat advection. We demonstrate second order accuracy in the $L^1$, $L^2$, and $L^\infty$ norms for both analytic test problems and numerical convergence studies. We also demonstrate the ability of the scheme to convert one chemical species to another across a moving boundary.

§ INTRODUCTION

The convective transport and diffusion of chemical species occurs in a broad range of systems. Many applications involve chemical concentrations that evolve within complex, time-dependent regions, including blood flow and clotting in the cardiovascular system <cit.>, particulate and chemical vapor transport in the lungs <cit.>, and drug absorption in the digestive tract <cit.>. In some cases, critical interactions occur between fluid-phase and structure-bound chemicals. These interactions can appear in the model equations as a Robin boundary condition for the fluid-bound chemical. For example, the interaction between circulating proteins and membrane-bound proteins play a pivotal role in thrombus formation <cit.>. Modeling these interactions becomes even more challenging as one considers the motion of the flow domain itself.

The numerical simulation of PDEs in complex domains has garnered significant attention for decades. Embedded boundary methods are popular approaches to such problems in which a fixed rectangular Cartesian grid is overlayed on the complex structure. Embedded boundary approaches typically alter the PDE to include an additional source term that is non-zero only near the boundary. For instance, the immersed boundary method <cit.> uses an integral transform with a regularized delta function kernel to enforce  boundary conditions along irregular interfaces immersed in a background Cartesian grid. These methods typically are designed for Dirichlet boundary conditions, and they have to be modified with special interpolation procedures to allow for other types of boundary conditions <cit.>. Volume penalization <cit.> methods and diffuse domain <cit.> methods both introduce phase field models to track the interface. A smoothed version of the interface is then used to modify the original PDE to account for the boundary conditions, ultimately yielding an approach similar to the immersed boundary method. These methods all have the effect of smoothing the boundary of the complex domain over several grid cells. Regularization lowers solution accuracy near the boundary, and this can limit the effectiveness to applications in which boundary interactions are important. One alternative approach is the immersed interface method <cit.>, which derives jump conditions across the interface and then builds these jump conditions into the discretized equations. Using this method to impose boundary conditions requires relating the boundary conditions to jump conditions along the interface. Many applications have succesfully used the immersed interface method, but to our knowledge, the jump conditions for arbitrary spatially dependent Robin boundary conditions have not yet been determined. A closely related approach is the ghost fluid method <cit.>, which uses an interpolation procedure to build stencils near the interface that incorporate the boundary conditions. This approach can lead to non-symmetric and, in some cases, ill-conditioned systems that require specialized solvers <cit.>. Another alternative to these methods is the immersed boundary smooth extension (IBSE) method <cit.>, which can solve PDEs on complex geometries by embedding the geometry in a simpler region and solving the PDE on the extended domain that now includes both “physical” and “non-physical” subdomains. In the IBSE approach, a body force is incorporated in the “non-physical" regime to extend the physical solution smoothly outside the physical domain. This allows for high order accuracy to the boundary, but at the cost of solving an additional multiharmonic problem, with the order depending on the number of interface conditions.

The approach that we use here to impose boundary conditions along a complex boundary while retaining a Cartesian grid discretization framework is based on a cut-cell finite volume formulation. In these flux-based methods, fluxes are carefully calculated to account for the portion of cells that are inside the physical domain. This approach allows for accurate solutions along the boundary, but comes at the expense of computing cell geometries at the interface. This expense can be alleviated through the use of a level set function to track the surface interface <cit.>. Recent work by Helgadottir  <cit.> demonstrated the ability to impose Dirichlet, Neumann, and Robin boundary conditions for Poisson problems. Cut-cell methods can suffer from the so-called “small cell problem," however, in which cells with small volumes necessitate the use of extremely small time step sizes for conditionally stable time stepping schemes. This becomes a serious problem for hyperbolic equations for which accurate, efficient, and unconditionally stable implicit time-stepping schemes are difficult to create. Approaches to alleviate the small cell problem for advective PDEs involve using an implicit method only for the cut-cells <cit.>, merging small cells with their neighbors to effectively create a larger cell <cit.>, or partitioning fluxes into “shielded" and “unshielded" zones based on the geometry of the cut-cell <cit.>. These methods have been succesfully deployed in two spatial dimensions, but extending them to three spatial dimensions remains a challenge. Recently, cut-cell methods that do not suffer from the small cell problem have been introduced to handle moving boundaries <cit.>; however, a formulation involving Robin boundary conditions that achieves second order accuracy has yet to be developed. The major contribution of this study is the construction of such a method with second order accuracy.

Herein, we develop a cut-cell method for advection-diffusion equations on moving domains that allows for the imposition of general Robin boundary conditions. To avoid the small cell problem, we introduce a split cut-cell semi-Lagrangian scheme, in which the diffusive operator is handled using well established finite volume methods, and the advective operator is treated using a semi-Lagrangian method. The benefit of splitting the two operators is two fold. First, in the diffusive solve, the flux from the boundary does not need to account for the change in location of the boundary. This allows us to leverage recent work by Papac  <cit.> and Bochkov  <cit.> to solve a Poisson-like problem with stationary boundaries. Second, the semi-Lagrangian scheme for advection has no CFL stability constraint, so the small cell problem is no longer an issue. We demonstrate that this method can accurately resolve concentrations near the boundary for both Robin and Neumann boundary conditions, including inhomogeneous Robin boundary conditions. In addition, we show that this method is able to handle conversion of one concentration into another across the boundary, effectively modelling surface reactions.

§ CONTINUOUS EQUATIONS

We consider the advection and diffusion of a chemical species with concentration $q\parens{\xx,t}$ in an arbitrary domain $\Omega_t$, embedded in a larger rectangular region $\mathcal{B}$. Both $q\parens{\xx,t}$ and $\Omega_t$ are transported by the same velocity field $\uu\parens{\xx,t}$. We assume that $q\parens{\xx,t}$ diffuses with diffusion coefficient $D$. We consider the particular case in which $q\parens{\xx,t}$ satisfies Robin boundary conditions on the boundary $\Gamma_t$ of $\Omega_t$, so that

\begin{align}
\frac{\partial q\parens{\xx,t}}{\partial t} + \grad\cdot\Big(\uu\parens{\xx,t} q\parens{\xx,t} - D\grad q\parens{\xx,t}\Big) &= f\parens{\xx,t}, \xx\in\Omega_t, \label{eq:main}\\
D\grad q\parens{\xx,t}\cdot\nn\parens{\xx,t} + a\parens{\xx,t} q\parens{\xx,t} &= g\parens{\xx,t}, \xx\in \Gamma_t, \label{eq:bc}
\end{align}

in which $\nn\parens{\xx,t}$ is the outward unit normal of $\Gamma_t$ and $f\parens{\xx,t}$ is a given volumetric source function. In the implementation, $q\parens{\xx,t}$ is defined only inside $\Omega_t$. The bounding region $\mathcal{B}$ is used only to define the cells contained within $\Omega_t$.

We describe the boundary $\Gamma_t$ using a signed distance function $\phi\parens{\xx,t}$ such that

\begin{equation}
\Gamma_t = \{\xx\in\BB |~ \phi\parens{\xx,t} = 0\}.
\end{equation}

The signed distance function is passively advected by the prescribed velocity,

\begin{equation} \label{eq:level_set_adv}
\frac{\partial \phi\parens{\xx,t}}{\partial t} + \grad\cdot\parens{\uu\parens{\xx,t} \phi\parens{\xx,t}} = 0, \xx\in\mathcal{B}.
\end{equation}

Because $\phi\parens{\xx,t}$ will not in general remain a signed distance function under advection by (<ref>), a reinitialization procedure is used to maintain the signed distance property. This can be achieved by computing the steady-state solution to the Hamilton-Jacobi equation

\begin{align}\label{eq:reinitialization}
\frac{\partial \hat{\phi}\parens{\xx,\tau}}{\partial \tau} + \text{sgn}\parens{\hat{\phi}\parens{\xx,\tau}}\parens{\norm*{\grad\hat{\phi}\parens{\xx,\tau}} - 1} &= 0, \\
\hat{\phi}\parens{\xx,0} &= \phi\parens{\xx,t}, \nonumber
\end{align}

after which $\phi\parens{\xx,t}$ is set to the steady state solution.

§ NUMERICAL METHODS

To simplify notation, we describe the method in two spatial dimensions. Extensions to a third spatial dimension are straightforward, with the exception of calculating cut-cell geometries <cit.>.

We use a splitting scheme to split equation (<ref>) into a diffusion step,

\begin{align}
\frac{\partial q\parens{\xx,t}}{\partial t} - \grad\cdot\parens{D\grad q\parens{\xx,t}} &= f\parens{\xx,t}, \xx\in\Omega_t \label{eq:main_diff}\\
D\grad q\parens{\xx,t}\cdot\nn\parens{\xx,t} + a\parens{\xx,t} q\parens{\xx,t} &= g\parens{\xx,t}, \xx\in \Gamma_t, \label{eq:bc_diff}
\end{align}

where the domain $\Omega_t$ remains fixed, followed by an advection step,

\begin{align}
\frac{\partial \phi\parens{\xx,t}}{\partial t} &+ \grad\cdot\parens{\uu\phi\parens{\xx,t}} = 0 \mbox{ for } \xx\in\BB \\
\frac{\partial q\parens{\xx,t}}{\partial t} &+ \grad\cdot\parens{\uu q\parens{\xx,t}} = 0 \mbox{ for } \xx\in\Omega_t,
\end{align}

in which we evolve both the level set and the function $q\parens{\xx,t}$. We note that because $\Omega_t$ and $q\parens{\xx,t}$ evolve with the same velocity, no boundary condition is needed with equation (<ref>). The use of this splitting procedure inccurs an additional error cost. This cost can be reduced to second order temporal accuracy using Strang splitting <cit.>, which will be described later.

We overlay a Cartesian grid on top of $\BB$ such that $\BB$ consists of rectangular grid cells $\cc_{i,j}$ and $\BB = \cup \cc_{i,j}$ with a grid cell spacing of $\Delta x = \Delta y = h$. The concentration field is approximated at the cell centroid of each full or partial cell contained within $\Omega_t$.

In the following sections, we describe the discretization of equations (<ref>) and (<ref>). In what follows, unless otherwise noted, $q_{i,j}$ refers to the concentration at the cell centroid of the cell $\cc_{i,j}\cap\Omega_t$.

§.§ Diffusion Step

To solve the diffusion step from equations (<ref>), we employ a cut-cell finite volume method based on the approaches of Papac  <cit.> and Arias  <cit.>. We summarize the derivation here, and refer the interested reader to a more detailed description in previous work. Integrating equation (<ref>) over a cell $\cc_{i,j}$ that is entirely or partially interior to $\Omega_t$ and dividing by the volume of the cell, we get

\begin{equation}\label{eq:integrate_diff}
\frac{1}{\size{\cc_{i,j}\cap\Omega_t}}\int_{\cc_{i,j}\cap\Omega_t} \frac{\partial q\parens{\xx,t}}{\partial t} \mathrm{d}\xx = \frac{1}{\size{\cc_{i,j}\cap\Omega_t}}\int_{\cc_{i,j}\cap\Omega_t} D\Lap q\parens{\xx,t} \mathrm{d}\xx.
\end{equation}

We define $Q_{i,j}$ as the cell average of $q\parens{\xx,t}$ in the cell $\cc_{i,j} \cap\Omega_t$. Replacing the cell average in the left side of equation (<ref>) and employing the divergence theorem on the right-hand side, we get

\begin{equation}\label{eq:diff_integral}
\frac{d Q_{i,j}}{d t} = \frac{1}{\size{\cc_{i,j}\cap\Omega_t}}\int_{\partial\parens{\cc_{i,j}\cap\Omega_t}} D\grad q\parens{\xx,t}\nn\cdot d\AA,
\end{equation}

in which $\nn$ is the outward unit normal of $\partial\parens{\cc_{i,j}\cap\Omega_t}$. We can further divide the integral in equation (<ref>) into an integral over the cell boundary $\partial \cc_{i,j}\cap\Omega_t$ and an integral over the physical boundary $\Gamma_t\cap\cc_{i,j}$

\begin{equation}
\int_{\partial\parens{\cc_{i,j}\cap\Omega_t}} D\grad q\parens{\xx,t}\cdot\nn dA = \parens{\int_{\partial \cc_{i,j}\cap\Omega_t} + \int_{\cc_{i,j}\cap\Gamma_t}} D\grad q\parens{\xx,t}\nn\cdot d\AA.
\end{equation}

We approximate the first integral by

\begin{align}\label{eq:first_integral}
\int_{\partial\cc_{i,j}\cap\Omega_t} D\grad q\parens{\xx,t}\cdot\nn dA \approx &L^\text{g}_{i+\frac{1}{2},j}\frac{\hat{q}_{i+1,j}-\hat{q}_{i,j}}{\Delta x} - L^\text{g}_{i-\frac{1}{2},j}\frac{\hat{q}_{i,j} - \hat{q}_{i-1,j}}{\Delta x} \\
+ &L^\text{g}_{i,j+\frac{1}{2}}\frac{\hat{q}_{i,j+1}-\hat{q}_{i,j}}{\Delta y} - L^\text{g}_{i,j-\frac{1}{2}}\frac{\hat{q}_{i,j}-\hat{q}_{i,j-1}}{\Delta y}, \nonumber
\end{align}

in which $\hat{q}_{i,j}$ is the point-wise concentration at the center of the cell $\cc_{i,j}$ and $L^\text{g}_{i+\frac{1}{2},j}$ is the length fraction of the face $\parens{i+\frac{1}{2}}\times \left[j-\frac{1}{2},j+\frac{1}{2}\right]$ covered by the irregular domain, see Figure <ref>. It is challenging to compute $L^\text{g}_{i+\frac{1}{2},j}$ exactly. The smoothness of the level set is determined by the zero contour. For shapes that have sharp features, the level set will inherit these features. However, to achieve second order accuracy for smooth level sets, it suffices to use a linear approximation

\begin{equation}\label{eq:L}
L^\text{g}_{i+\frac{1}{2},j} = \left\{ \begin{array}{cc}
\Delta y \size{\frac{\phi_{i+\frac{1}{2},j-\frac{1}{2}}}{\phi_{i+\frac{1}{2},j-\frac{1}{2}}-\phi_{i+\frac{1}{2},j+\frac{1}{2}}}} &\mbox{ if } \phi_{i+\frac{1}{2},j-\frac{1}{2}}\cdot\phi_{i+\frac{1}{2},j+\frac{1}{2}} < 0 \\
\Delta y &\mbox{ if } \phi_{i+\frac{1}{2},j-\frac{1}{2}} < 0 \mbox{ and } \phi_{i+\frac{1}{2},j+\frac{1}{2}} < 0 \\
0 &\mbox{ otherwise.}

\end{array}\right.
\end{equation}

While in the evolution of the level set $\phi\parens{\xx,t}$, the degrees of freedom live at cell centers, in equation (<ref>), we require the value of $\phi\parens{\xx,t}$ at cell nodes. In our computations, we use a simple bi-linear interpolant to find these nodal values. In three spatial dimensions, evaluating the corresponding quantity $L^\text{g}_{i+\frac{1}{2},j,k}$ would involve computing the surface area.

dy$\Delta y$
[width=0.35,hsep=-1.25em,vsep=1em,pos=ul](a)Graphs/drawing.eps    [width=0.35,hsep=-1.25em,vsep=1em,pos=ul](b)Graphs/close_up.eps

Depiction of the nomenclature used in the diffusion discretization for the cell length fractions (fig:discretization:diff) and for the boundary fluxes (fig:discretization:boundary).

We emphasize that $\hat{q}_{i,j}$ refers to the cell center of the (possibly cut) grid cell $\cc_{i,j}$ regardless of the location of the physical boundary $\Gamma_t$. In cases where the cell center does not correspond to the location of degrees of freedom, e.g. near cut-cells, we must reconstruct these values. Here, we perform this reconstruction using either a moving least squares (MLS) approximation, or a radial basis function (RBF) interpolant. Both procedures are described in more detail in Section <ref>. We note that this reconstruction is not required in the methods of Arias  <cit.> and Bochkov  <cit.>, as those methods always define the degrees of freedom at cell centers.

As done in Bochkov  <cit.>, the second integral is approximated using a linear approximation to $q\parens{\xx,t}$ on the boundary in the direction normal to $\Gamma_t$

\begin{equation}\label{eq:second_integral}
\int_{\cc_{i,j}\cap\Gamma_t} D\grad q\parens{\xx,t}\cdot\nn dA = \int_{\cc_{i,j}\cap\Gamma_t} \parens{g\parens{\xx,t} - a q\parens{\xx,t}}dA \approx L^\Gamma_{i,j}\parens{g\parens{\rr_{i,j}, t} - a q\parens{\rr_{i,j}, t}},
\end{equation}

in which $\rr_{i,j}$ is the closest point on $\Gamma_t$ to $\xx_{i,j}$ and $L^\Gamma_{i,j} = \size{\cc_{i,j}\cap\Gamma_t}$ (see Figure <ref>). The value $q\parens{\rr_{i,j}, t}$ is found using a Taylor series expansion

\begin{equation}\label{eq:bdry_1}
q\parens{\xx_{i,j}, t} = q\parens{\rr_{i,j}, t} + d_{i,j}\frac{\partial q\parens{\rr_{i,j}, t}}{\partial \nn} + \bigOh{h^2},
\end{equation}

in which $d_{i,j} = \frac{\phi\parens{\xx_{i,j}, t}}{\size{\grad\phi\parens{\xx_{i,j}, t}}}$ is the distance between $\rr_{i,j}$ and cell center $\xx_{i,j}$. On the boundary, we have that

\begin{equation}\label{eq:bdry_2}
D\frac{\partial q\parens{\rr_{i,j},t}}{\partial \nn} + a q\parens{\rr_{i,j}, t} = g\parens{\rr_{i,j},t}.
\end{equation}

We solve equations (<ref>) and (<ref>), dropping the $\bigOh{h^2}$ terms, for $q\parens{\rr_{i,j},t}$ and use this value in the approximation to the integral in equation (<ref>). As shown previously <cit.>, this system is well posed for a sufficiently refined grid.

The cut-cell volume $\size{\cc_{i,j}\cap\Omega_t}$ and physical boundary length $\size{\cc_{i,j}\cap\Gamma_t}$ are found by decomposing the cut-cell and boundary into simplices, for which analytic formulas for the volume exist <cit.>.

We approximate equation (<ref>) using the implicit trapezoidal rule:

\begin{align}\label{eq:diff_timestepping}
\frac{1}{\Delta t} Q_{i,j}^{n+1} &- \frac{D}{2\size{\cc_{i,j}\cap\Omega^-}}\left(L^\text{g}_{i+\frac{1}{2},j}\frac{q^{n+1}_{i+1,j}-q^{n+1}_{i,j}}{\Delta x} - L^\text{g}_{i-\frac{1}{2},j}\frac{q^{n+1}_{i,j} - q^{n+1}_{i-1,j}}{\Delta x}\right.& \nonumber \\
&\left.+ L^\text{g}_{i,j+\frac{1}{2}}\frac{q^{n+1}_{i,j+1}-q^{n+1}_{i,j}}{\Delta y} - L^\text{g}_{i,j-\frac{1}{2}}\frac{q^{n+1}_{i,j}-q^{n+1}_{i,j-1}}{\Delta y}\right)  \nonumber \\
&= \frac{1}{\Delta t} Q_{i,j}^n+ \frac{D}{2\size{\cc_{i,j}\cap\Omega^-}}\left(L^\text{g}_{i+\frac{1}{2},j}\frac{q^n_{i+1,j}-q^n_{i,j}}{\Delta x} - L^\text{g}_{i-\frac{1}{2},j}\frac{q^n_{i,j} - q^n_{i-1,j}}{\Delta x}\right. \\
&\left.+ L^\text{g}_{i,j+\frac{1}{2}}\frac{q^n_{i,j+1}-q^n_{i,j}}{\Delta y} - L^\text{g}_{i,j-\frac{1}{2}}\frac{q^n_{i,j}-q^n_{i,j-1}}{\Delta y}\right)+ \frac{1}{\size{\cc_{i,j}\cap\Omega^-}}f_{i,j}^{n+\frac{1}{2}}. \nonumber
\end{align}

In our computational examples, we use GMRES without preconditioning to solve this system of equations. The solver typically converges to a solution with tolerance $10^{-12}$ after approximatley 20 to 50 iterations.

§.§ Advection Step

For the advection step, we use a semi-Lagrangian method <cit.> to advance both the level set $\phi\parens{\xx,t}$ and the concentration field $q\parens{\xx,t}$. We solve equation (<ref>) by first finding the preimage $\XX_{i,j} = \bchi^{-1}\parens{\xx_{i,j},t^{n+1}}$ of the fluid parcel located at $\xx_{i,j}$ at time $t^{n+1}$. The mapping $\bchi\parens{\xx,t}$ satisfies the differential equation

\begin{equation}\label{eq:material_eq}
\frac{\partial\bchi\parens{\xx,t}}{\partial t} = \uu\parens{\xx,t} \mbox{ for } \xx\in\BB.
\end{equation}

The preimage $\XX_{i,j}$ can be found by integrating equation (<ref>) backwards in time. We use an explicit two step Runge Kutta method:

\begin{align}
\XX^\star_{i,j} = \xx_{i,j} - \frac{\Delta t}{2}\uu\parens{\xx_{i,j}, t^{n+1}}, \\
\XX_{i,j}^n = \xx_{i,j} - \Delta t \uu\parens{\XX^\star_{i,j},t^{n+\frac{1}{2}}}.
\end{align}

Having found the preimage, we can interpolate the solution at time $t^n$ at the preimage locations $\XX_{i,j}^n$.

Our interpolation procedure consists of one of two possible methods depending on whether $\XX_{i,j}^n$ is near cut-cells. Away from cells pierced by the zero level set, we use a tensor product of special piecewise Hermite polynomials called Z-splines <cit.>. The Z-spline $Z_m\parens{x}$ interpolating the data $\parens{x_i,f_i = f\parens{x_i}}_{i = 1}^n$ is a piecewise polynomial function that satisfies

\begin{align}
Z_m\parens{x}&\in C^m\parens{\left[x_1,x_n\right]}, \\
\left.\frac{d^p}{d x^p} Z_m\parens{x}\right\rvert_{x_j} &= f_{m,j}^p \mbox{ for } p=0,\ldots,m \mbox { and } j=1,\ldots,n, \\
Z_m\parens{x}&\in \pi_{2m+1}\parens{\left[x_i,x_{i+1}\right]} \mbox{ for } i=1,\ldots,n-1,
\end{align}

in which $f_{m,j}^p$ is the approximation of the $p^\text{th}$ order derivative of $f(x)$ computed from high-order finite differences of $f_j$ using $2m + 1$ points, and $\pi_{2m+1}$ is the space of polynomials of degree less than or equal to $2m+1$.

It is possible to define $Z_m$ for a general set of data points (or abscissae) via the cardinal Z-splines <cit.>

\begin{equation}\label{eq:interpolant}
Z_m\parens{x} = \sum_i f_i\tilde{Z}_m\parens{x - x_i}.
\end{equation}

The cardinal Z-splines have two key properties. First, the cardianl Z-splines have compact support, $\tilde{Z}_m\parens{x} = 0$ for $\size{x} > m+1$, which allows for efficient local evaluation of interpolants. Second, cardinal Z-splines are interpolatory, which makes the interpolant trivial to form. Because data values are defined on a regular Cartesian grid, the computation of the full interpolant uses a tensor product of cardinal Z-splines

\begin{equation}\label{eq:interpolant_full}
Z_m\parens{\xx} = \sum_i\sum_j q_{i,j}\tilde{Z}_m\parens{x - x_i}\tilde{Z}_m\parens{y - y_j}.
\end{equation}

In our computations, we use the quintic Z-splines, which are defined in terms of

\begin{equation}
Z_2\parens{x} = \left\{ \begin{array}{cc}
1 - \frac{15}{12} x^2 - \frac{35}{12} x^3 + \frac{63}{12} x^4 - \frac{25}{12} x^5 &\mbox{ if } \size{x} \leq 1 \\
-4 + \frac{75}{4} x - \frac{245}{8} x^2 + \frac{545}{24} x^3 - \frac{63}{8} x^4 + \frac{25}{24} x^5 &\mbox{ if } 1 < \size{x} \leq 2 \\
18 - \frac{153}{4} x + \frac{255}{8} x^2 - \frac{313}{24} x^3 + \frac{21}{8} x^4 - \frac{5}{24} x^5 &\mbox{ if } 2 < \size{x} \leq 3 \\
0 &\mbox{ otherwise.}
\end{array}\right.
\end{equation}

In cases where we do not have enough points to define the Z-spline, e.g., near cut-cells, we use the procedure described in Section <ref>. We note that although the interpolation procedure described here has a discontinuous switch between operators, this does not appear to affect the overall convergence rates in our numerical tests of the methodology.

We note that this form of the semi-Lagrangian method is not conservative, but, in our experiments, the change in the amount of material was very small. Conservative verions of semi-Lagrangian methods have been developed <cit.>. The core change is to advect grid cells instead of individual points, and then integrate the resulting interpolating polynomial over the grid cell. Because conservation is not critical to our ultimate problem of interest, we use the simpler non-conservative approach in this work.

The evolution of the level set and the concentration use the same procedure, except for the choice of reference grid. The reference grid for the level set consists of the entire domain $\BB$, whereas the reference grid for the concentration consists of the domain $\Omega_{t^{n+1}} = \left\{ \xx\in\BB | \phi\parens{\xx,t^{n+1}} < 0\right\}$. Therefore, we update the level set prior to updating the concentration.

§.§ Moving Least Squares and Radial Basis Function Interpolation

Near cut-cells, it is necessary to form interpolants on unstructured data. In this study, we compare the accuracy of a moving least squares (MLS) approximation and a local radial basis function (RBF) interpolant. This procedure is used both in the diffusion step, to extrapolate data from cut-cell centroids to full-cell centers, and also in the advection step, for cells where the Z-spline interpolant can not be formed. We describe the procedure in the context of reconstructing a function $f\parens{\xx}$ from the data points $\{\xx_i,f\parens{\xx_i}\}$ with $i = 1,2,\ldots,n$ in which $n$ is an arbitrary number of points.

The MLS approximation is formed by finding the best approximation $q\parens{\xx}$ at a point $\xx_\text{c}$ of the form

\begin{equation}
q\parens{\xx} = \sum_{j = 1}^{N} c_j p_j\parens{\xx}
\end{equation}

in which the $p_j\parens{\xx}$ form a basis for the space of polynomials up to a certain degree and $N$ is the number of polynomials in the basis. The approximation $q\parens{\xx}$ is chosen to be optimal with respect to the standard weighted $L^2$ inner product with weight function $w\parens{\xx}$

\begin{equation}
\langle f\parens{\xx},g\parens{\xx}\rangle = \int_\Omega f\parens{\xx}g\parens{\xx}w\parens{\xx}d\xx.
\end{equation}

We use a stencil width of approximately two grid cells. It has been shown that the weight function must be singular at the data locations $\xx_i$ for the MLS approximation to interpolate the data <cit.>. Here, we use

\begin{equation}
w\parens{\xx} = e^{-\norm*{\xx - \xx_\text{c}}^2}.
\end{equation}

Because this weighting function is not singular at $\xx_\text{c}$, the reconstructed polynomial will not be interpolatory. The MLS calculation reduces to solving a linear system of the form

\begin{equation}
\Lambda A \cc = \Lambda \ff
\end{equation}

in which $A$ is a matrix whose entries consist of $a_{i,j} = p_j\parens{\xx_i}$ and $\Lambda$ is a diagonal matrix with $\lambda_{i,i} = w\parens{\xx_i}$. In preliminary computational tests, we observed that using a lower order reconstruction for the diffusion step is more robust without affecting the order of accuracy while quadratic reconstructions can be used in the advection step. Consequently, we use a quadratic polynomial in the reconstructions for the advection step and a linear reconstruction with the diffusion step for the tests reported in section <ref>.

The RBF approximation is constructed via a polyharmonic spline of the form

\begin{equation}
q\parens{\xx} = \sum_{j = 1}^n \lambda_j\phi\parens{\norm*{\xx-\xx_j}} + \sum_{j = 1}^s \beta_j p_j\parens{\xx}
\end{equation}

in which $\phi\parens{\xx}$ is a polyharmonic radial basis function of degree $m$ and $p_j\parens{\xx}$ are a set of $s$ polynomial basis functions <cit.>. The degree of the polynomials $k$ is chosen such that $m = 2k + 1$. The coefficients are chosen so that

\begin{align}
q\parens{\xx_j} &= f_j \mbox{ for } j = 1,\ldots,n,\\
\sum_{k=1}^n\lambda_j p_k\parens{\xx_j} &= 0  \mbox{ for } j = 1,\ldots,n.
\end{align}

This leads to the linear system for the coefficients

\begin{equation}
\left(\begin{array}{cc}
A & P \\
P^T & 0 \\
\end{array}\right)
\left(\begin{array}{c}
\lambda \\
\beta
\end{array}\right) =
\left(\begin{array}{c}
\ff \\
\end{array}\right)
\end{equation}

in which $A$ is a square matrix with elements $A_{i,j} = \phi\parens{\norm*{\xx_i - \xx_j}}$, $P$ is a rectangular matrix with elements $P_{i,j} = p_j\parens{\xx_i}$. In our computations, we use the polyharmonic spline $\phi\parens{\xx} = \norm*{\xx}^3$ combined with linear polynomials.

§.§ Full procedure and Implementation

We summarize the full procedure to advance the solution from time $t^{n}$ to time $t^{n+1}$:

* If needed, reinitialize the level set $\phi\parens{\xx,t}$ into the signed distance function by iterating equation (<ref>) to steady state using the procedure as described by Nangia  <cit.>.

* Update the diffusion equation to half time $t^{n+\frac{1}{2}}$ using the methods described in Section <ref>.

* Update the level set using the prescribed velocity $\uu$ to time $t^{n+1}$ by the methods described in Section <ref>.

* Advect the concentration using the prescribed velocity $\uu$ by the methods described in Section <ref>.

* Update the diffusion equation to full time $t^{n+1}$ using the concentration and level sets from the previous two steps.

The above procedure is implemented using the SAMRAI <cit.> infrastructure, which provides an efficient, parallelized environment for structured adaptive mesh refinement. The diffusion solve is computed using matrix free solvers with operators provided by IBAMR <cit.> and Krylov methods provided by PETSc <cit.>.

§ RESULTS

Here we demonstrate the capabilities of the method both using a prescribed level set and an level set advected with the fluid. We start with diffusion dominated examples before exploring inclusion of advection and spatially varying boundary conditions.

§.§ Diffusion from a point source

We consider the advection-diffusion of a point source within a disk of radius $R$. The disk is advected with velocity $\uu = \parens{\cos\parens{\frac{\pi}{4}},\sin\parens{\frac{\pi}{4}}}$. The exact concentration is

\begin{equation}\label{eq:exact}
q\parens{\xx,t} = \frac{10}{4 D \parens{t+\frac{1}{2}}}e^{-\frac{\norm*{\xx - \xx_\text{c}\parens{t}}^2}{4D\parens{t+\frac{1}{2}}}} \mbox{ for } \xx \in \BB,
\end{equation}

in which $\xx_\text{c}\parens{t}$ is the center of the disk. We apply Robin boundary conditions of the form

\begin{equation}\label{eq:bdry_conds}
D\frac{\partial q}{\partial\nn} = g\parens{\xx,t} - a\parens{\xx,t} q\parens{\xx,t}.
\end{equation}

We set $a\parens{\xx,t} = 1$ and use the method of manufactured solutions to determine the value of $g\parens{\xx,t}$. For this example, we specify the location of the level set $\phi$ at each timestep. This procedure allows us test the convergence of the advection-diffusion step without having to account for numerical error in the evolution of the level set function. We note that whereas this example includes a trivial semi-Lagrangian backwards integration step, the interpolation procedure and diffusion step are non-trivial.

The extended domain $\BB$ is the box $[0,12]\times[0,12]$ and is discretized using $N$ points in each direction. We use a disk radius of $R = 1$ and initial center $\xx_\text{c}(0) = \parens{1.521,1.503}$. The choice of center slightly offsets the disk from the background grid so that different cut-cells are generated on opposite sides of the disk. We use a diffusion coefficient of $D = 0.01$. The simulations are run at an advective CFL number of $C_{\text{CFL}} = \frac{\Delta x}{\size{\uu}\Delta t} = 0.5$.

Figure <ref> shows convergence plots, which demonstrate second order accuracy in the $L^1$, $L^2$, and $L^\infty$ norms. The coarsest simulation uses $N = 128$ grid points in each direction of the background Cartesian mesh, which corresponds to approximately 26 grid cells covering the diameter of the disk. Despite the relatively coarse description of the disk, we still see errors that are on the order of one percent.

We also consider the case where the level set $\phi$ is advected with the same velocity. In this case, numerical issues make it difficult to assess the pointwise error on the Cartesian grid because the piecewise linear reconstruction used to determine the cut-cells gives different cut-cell geometries than that of a prescribed level set. Instead, we compare the solution at the boundary to the imposed boundary condition. We compute the solution on the boundary using a moving least squares linear extrapolation from cut-cells to the cell boundary. Further, to assess the convergence rate, we compute the error of the surface integral

\begin{equation}\label{eq:Error_moving}
E = \size{\int_0^{2\pi} q\parens{\xx,t} d\theta - \int_0^{2\pi} \tilde{q}\parens{\xx,t} d\theta},
\end{equation}

in which $\theta$ is the angle along the boundary from the center of the disk, $q\parens{\xx,t}$ is the exact solution from (<ref>), and $\tilde{q}\parens{\xx,t}$ is the approximate solution computed using the moving least squares extrapolation as described in Section <ref>. A convergence plot is shown in Figure <ref> and indicates that the method achieves between first and second order convergence rates. In both the evolved and prescribed level set cases, the RBF reconstruction yields errors that are more than an order of magnitude smaller than the MLS reconstructions.

Convergence rates for an advecting and diffusion concentration field initialized from a point source in a constant flow using either a RBF (fig:adv_diff_pt_source:rbf) or MLS (fig:adv_diff_pt_source:mls) reconstruction. Also shown is the convergence rate for the integral in equation (<ref>) using either an RBF or MLS reconstruction (fig:adv_diff_pt_source:bdry). In all cases, the error using RBF reconstructions are more than an order of magnitude smaller than MLS reconstructions. Simulations are run to a final time of $T = 10$ and at a CFL number of $C_\text{CFL} = 0.5$. We use a diffusion coefficient of $D = 0.01$.

§.§ Solid body rotation

In this section, we again consider diffusion from a point source, but instead apply a solid body rotation to the disk. In this case, the semi-Lagrangian sub-step contains a non-trivial integration in time, and therefore tests all components of the solver. We again prescribe the initial condition as in equation (<ref>), and apply boundary conditions as in equation (<ref>).

The computational domain $\BB$ is $[-4,4]\times[-4,4]$ and is again discretized using $N$ points in each direction. We use a disk radius of $R = 1$ and initial center $\xx_c(0) = \parens{1.521, 1.503}$. We use a diffusion coefficient of $D = 0.1$. The prescribed velocity is $\uu = 2\pi\parens{-y,x}$. Simulations are run at an advective CFL number of $0.5$ to a final time of $T = 1$. Figure <ref> shows the initial and final solution. We can see that the solution maintains symmetry throughout the simulation. Figure <ref> shows convergence rates for both RBF and MLS reconstructions. Figures <ref> and <ref> show convergence rates in which the level set is prescribed. We see second order convergence rates in all norms for the prescribed level set case. Figure <ref> shows the convergence rate for the solution extrapolated to the boundary for the case in which the level set is evolved with the fluid velocity. In this case, we see approximately second order convergence rates. In both of these cases, the error for the RBF reconstruction is almost an order of magnitude lower than that of the MLS reconstruction.

Initial (fig:rotational_init_final:initial) and final (fig:rotational_init_final:final) configuration for a diffusing concentration field initialized from a point source in a rigid body rotational flow field. The zero contour of $\phi$ is shown in red. Simulations are run at an advective CFL number of $C_\text{CFL} = 0.5$ to a final time of $T = 1$.

Convergence rates for advecting and diffusion concentration field initialized from a point source in a solid body rotational flow using either a RBF (fig:adv_diff_rotational:rbf) or MLS (fig:adv_diff_rotational:mls) reconstruction. Also shown is the convergence rate for the integral in equation (<ref>) using both an RBF or MLS reconstruction (fig:adv_diff_rotational:bdry). Simulations are run at an advective CFL number of $C_\text{CFL} = 0.5$ to a final time of $T = 1$.

§.§ Advection diffusion in oscillatory Couette flow

We now consider the advection and diffusion of a concentration inside a vesicle under oscillating Couette flow with no flux boundary conditions. Specifically, we specify the rotational component of the velocity field as

\begin{align}
u_\theta\parens{r} &= \parens{a r + \frac{b}{r}}\sin\parens{\pi t}, \\
a &= \frac{\Omega_2 R_2^2 - \Omega_1 R_1^2}{R_2^2 - R_1^2}, \\
b &= \frac{\parens{\Omega_1 - \Omega_2}R_1^2R_2^2}{R_2^2 - R_1^2},
\end{align}

where $R_1$ and $R_2$ are the radii of the coaxial cylinders, and $\Omega_1$ and $\Omega_2$ are the angular velocities of the cylinders. Here, we set the radii to be $R_1 = 0.5$ and $R_2 = 3.75$ with respective angular velocities $\Omega_1 = 10.0$ and $\Omega = 1.0$.

The vesicle is described by the level set $\phi_\text{v}\parens{\xx,t}$ so that the initial condition is

\begin{equation}
\phi_\text{v}\parens{\xx,0} = \norm*{\xx - \xx_{\text{v},\text{c}}} - R
\end{equation}

where $R = 1$ is the radius of the disk, and $\xx_{\text{v},\text{c}} = \parens{1.521,1.503}$ is the center of the disk. We initialize the concentration field $q_\text{v}$ via

\begin{equation}\label{eq:no_flux_inside}
q_\text{v}\parens{\xx,0} = \parens{\cos\parens{\pi \norm*{\xx-\xx_{\text{v},\text{c}}}} + 1}^2
\end{equation}

For comparison, we define another concentration field $q_\text{o}\parens{\xx,t}$ outside the vesicle, but inside the cylinders. The second level set is defined by

\begin{equation}
\phi_\text{o}\parens{\xx,t} = \max\parens{-\phi_\text{v}\parens{\xx,t},\norm*{\xx} - R_1, R_2 - \norm*{\xx}}.
\end{equation}

The initial condition is given by

\begin{equation}
q_\text{o}\parens{\xx,0} = \left\{\begin{array}{cc}
\parens{\cos\parens{\pi \norm*{\xx-\xx_{\text{o},\text{c}}}} + 1}^2 &\mbox{ if } \norm*{\xx-\xx_{\text{o},\text{c}}} \leq 1 \\
0 &\mbox{ otherwise,}
\end{array}\right.
\end{equation}

in which $\xx_{\text{o},\text{c}} = - \xx_{\text{v},\text{c}}$.

The computational domain $\BB$ is the box $[-4, 4]\times [-4,4]$ and is discretized using $N$ points in each direction, and the diffusion coefficient for both the interior and exterior concentration fields is $D = 0.05$. Figure <ref> shows the solution at five different time values during the simulation. Over the course of the simulations, the concentration field inside the disk stays completely contained within the disk, while the outside concentration field is free to diffuse between the coaxial cylinders.

Figures <ref> and <ref> show numerical convergence results for the concentrations inside the vesicle and outside the vesicle but inside the two cylinders respectively. We estimate the convergence rate $r$ using Richardson extrapolation via

\begin{equation}
r = \log_2\frac{\norm*{q_{4h} - q_{2h}}}{\norm*{q_{2h} - q_{h}}},
\end{equation}

in which $q_h$ is the solution with a grid spacing of $h$. To compute the difference $\norm*{q_{2h} - q_{h}}$, we first interpolate the fine solution $q_h$ onto the coarser grid $q_{2h}$. Because the interpolation from a fine grid to a coarse grid is nontrivial near cut-cells, we exclude all cut-cells and any cells that neighbor cut-cells in the computation of the norm. To assess the error on the cut-cells, we perform a reconstruction of the the solution on the surface of the disk by extrapolating the solution to the disk in each cut-cell, then computing a cubic spline of these data points. Figure <ref> shows convergence of this reconstruction for the solution inside the vessicle. We obtain second order convergence rates using RBFs at all points in time for the finest grids. In contrast, the MLS reconstruction does not show optimal convergence rates for the range of grid spacings considered. We expect that under additional grid refinement the MLS approximation will settle to second order accuracy, but this is not achieved over the range of grid spacings considered herein. The dips in convergence rates in Figure <ref> occur when significant values of $q_\text{o}$ begin to appear in cut cells. Additionally, Figures <ref> and <ref> show the norms of the differences. To quantify the difference in accuracy between MLS and RBF reconstructions, we estimate the error coefficient $C$ in the expression

\begin{equation}
\norm*{q_h - q} \approx C h^p,
\end{equation}

in which $q$ is the exact solution and $p$ is the order of accuracy. We approximate $C$ by comparing differences between two grids of refinement $h$ and $2h$ through the following equation

\begin{equation}\label{eq:C}
C \approx \frac{\norm*{q_h - q_{2h}}}{h^p\parens{1 - 2^p}}.
\end{equation}

For both MLS and RBF reconstructions, we have $p = 2$. Table <ref> shows our estimates of the $C$ coefficient for both reconstructions. In all cases, the RBF reconstruction is more accurate.

[width=0.45,hsep=0em,vsep=2em,pos=ul](a)Graphs/couette/couette_0.eps  [width=0.45,hsep=0em,vsep=2em,pos=ul](b)Graphs/couette/couette_1.eps

[width=0.45,hsep=0em,vsep=2em,pos=ul](c)Graphs/couette/couette_2.eps  [width=0.45,hsep=0em,vsep=2em,pos=ul](d)Graphs/couette/couette_3.eps

Concentrations $q_\text{o}$ (orange) and $q_\text{v}$ (blue) at time points $0.0$ (fig:two_couette:0), $0.5$ (fig:two_couette:1), $1.0$ (fig:two_couette:2), $1.5$ (fig:two_couette:3), and $2.0$ (fig:two_couette:4). The zero contour of the level set $\phi_\text{o}$ is shown in green. The diffusion coefficient is fixed at $D = 0.05$ with $1024$ points in each direction. We use a fixed time-step corresponding to a maximum CFL number of $C_\text{CFL} = 0.5$ to a final time of $T = 2$.

[width=0.45,hsep=-1em,vsep=2em,pos=ul](a)Graphs/couette/inside_l1_rates.eps  [width=0.45,hsep=-1em,vsep=2em,pos=ul](b)Graphs/couette/inside_l2_rates.eps

[width=0.45,hsep=-1em,vsep=2em,pos=ul](c)Graphs/couette/inside_max_rates.eps  [width=0.45,hsep=-1em,vsep=2em,pos=ul](d)Graphs/couette/boundary_l1_rates.eps

Numerical convergence study for the interior of the disk in the $L^1$ (fig:couette_in_converge:l1_rate), $L^2$ (fig:couette_in_converge:l2_rate), and $L^\infty$ (fig:couette_in_converge:max_rate) norms. The convergence rates are computed from simulations with $N = 128,256,512$ points (blue) and $N = 256,512,1024$ points (orange). Also shown is the convergence rate for the extrapolation of the solution to the boundary (fig:couette_in_converge:bdry_l1_rate).

Numerical convergence study for the exterior of the disk but inside the two cylinders in the $L^1$ (fig:couette_out_converge:l1_rate), $L^2$ (fig:couette_out_converge:l2_rate), and $L^\infty$ (fig:couette_out_converge:max_rate) norms. The convergence study is computed from simulations with $N = 128,256,512$ points and $N = 256,512,1024$ points. The dip in convergence rates near the end of the simulation occur when there is a significant value of $q_\text{o}$ near the boundaries.

Norms of solution differences for the concentration field inside the domain for the $L^1$ (fig:couette_in_norms:l1_norms), $L^2$ (fig:couette_in_norms:l2_norms), and $L^\infty$ (fig:couette_in_norms:max_norms) norms. The dashed lines correspond to MLS reconstruction while the solid lines correspond to a RBF reconstruction. Also shown are the norms of solution differences for the extrapolation to the boundary (fig:couette_in_norms:bdry_norms).

Norms of solution differences for the concentration field outside the domain for the $L^1$ (fig:couette_out_norms:l1_norms), $L^2$ (fig:couette_out_norms:l2_norms), and $L^\infty$ (fig:couette_out_norms:max_norms) norms. The dashed lines correspond to MLS reconstruction while the solid lines correspond to a RBF reconstruction.

$q_\text{o}$             $q_\text{v}$

$L^1$ norm     $L^2$ norm     $L^\infty$ norm     $L^1$ norm     $L^2$ norm     $L^\infty$ norm

$C_\text{MLS}$     41.01     26.71     30.37     5.85     3.49     6.68

$C_\text{RBF}$      16.20      11.12     12.88     3.44     1.60     2.50

The leading order error coefficient $C$ for both MLS and RBF reconstructions is estimated using equation (<ref>). In all cases, the RBF reconstruction yields superior accuracy.

§.§ Choice of time step size

Previous sections use time step sizes based on a CFL number of $0.5$. We note that each substep consists of an unconditionally stable method. While the affect on stability of the splitting error is unclear, we briefly explore the choice of time step size in this section. We again use the oscillatory Couette flow example with radial basis function reconstructions as in the previous section; however, we now vary the CFL number. Recall we defined the CFL number as $C_{\text{CFL}} = \frac{\Delta x}{\size{\uu}\Delta t}$ in which $\size{\uu}$ is the maximum velocity over both space and time so that we used a fixed time step size throughout the entire simulation. As before, we estimate the $C$ coefficient (see equation (<ref>)), and the values are reported in table <ref>. We observe a decrease in the error as the CFL number increases. This is not surprising as the error term for the semi-Lagrangian method contains a term that scales like $\bigOh{\frac{1}{\Delta t}}$, so that the error can decrease for larger timesteps <cit.>, although this reduction in error will eventually break down.

$C_\text{CFL}$     0.5     1.0     1.5     2.0     2.5     3.0

$L^1$ norm          16.20     12.44     10.08     8.18     7.04     6.98

$L^2$ norm          11.12     8.95      7.47      6.38     5.69     5.25

$L^\infty$ norm     12.88     11.81     10.55     9.43     8.46     7.61

$C_\text{CFL}$     0.5     1.0     1.5     2.0     2.5     3.0

$L^1$ norm          3.44     3.01     2.86     2.64     2.40     2.18

$L^2$ norm          1.60     1.47     1.38     1.24     1.08     0.95

$L^\infty$ norm     2.50     2.65     2.49     2.17     1.81     1.60

The leading order error coefficient $C$ as estimated by equation (<ref>) as a function of the CFL number $C_\text{CFL}$. We see a decrease in error as the CFL number is increased.

§.§ Interaction of two concentrations in oscillatory Couette flow

We now consider the case of two fluid phase chemicals advecting and diffusing on two different domains, but interacting through a common boundary. Specifically, we consider oscillatory Couette flow of two chemicals, one of which is contained within a vesicle that passively advects with the flow and is initially in the shape of a disk. The other chemical exists outside the vesicle, but between the disks defining the domain of interest. Specifically, we solve the equations

∂q_v/∂t + ·q_v    = D q_v,        ∈Ω_v

-D∂q_v/∂    = κq_v - q_o,        ∈Γ_v

∂q_o/∂t + ·q_o    = D q_o,        ∈Ω_o

-D∂q_o/∂    = κq_o - q_v,        ∈Γ_v

-D∂q_o/∂    = 0,        ∈Γ_o∖Γ_v,

in which $\uu$ is defined in the previous section and $\Omega_\text{v}$ and $\Omega_\text{o}$ are the domains for the concentration inside and outside the vesicle, respectively. The initial concentration for $q_\text{v}$ is given by equation (<ref>), and for $q_\text{o}$ is initially uniformly zero.

We modify the time discretization of the diffusion step to use a modified trapezoidal rule for the boundary conditions. We first solve equation (<ref>) for initial approximations $\tilde{q}_\text{o}$ and $\tilde{q}_\text{v}$ using explicit approximations for the boundary conditions

\begin{align}
-D\frac{\partial \tilde{q}_\text{o}}{\partial \nn} = \kappa\parens{\tilde{q}_\text{o} - q^n_\text{v}}, \text{ and } \\
-D\frac{\partial \tilde{q}_\text{v}}{\partial \nn} = \kappa\parens{\tilde{q}_\text{v} - q^n_\text{o}}.
\end{align}

We then solve equation (<ref>) for $q_\text{o}^{n+1}$ and $q_\text{v}^{n+1}$ again using the intermediate result when evaluating the boundary conditions

\begin{align}
-D\frac{\partial q^{n+1}_\text{o}}{\partial \nn} = \kappa\parens{q_\text{o}^{n+1} - \tilde{q}_\text{v}}, \text{ and } \\
-D\frac{\partial q^{n+1}_\text{v}}{\partial \nn} = \kappa\parens{q_\text{v}^{n+1} - \tilde{q}_\text{o}},
\end{align}

in which $\tilde{q}_\text{v}$ and $\tilde{q}_\text{o}$ are the results from the intermediate result.

For this example, we perform only RBF reconstructions at the boundaries. We set $\kappa = 1$. The remaining parameters are the same as in the previous section. Concentration values for various points in time are shown in Figure <ref>. As $q_\text{v}$ reaches the boundary $\Gamma_\text{v}$, it is converted to $q_\text{o}$ at a rate proportional to the difference in concentrations.

As in the previous section, we perform a numerical convergence study for both the interior and exterior fields. Figures <ref> and <ref> show the convergence rates as a function of time. The convergence rates as we refine the grid appear to be converging toward second order, although the method has not yet settled down into it's asymptotic regime. We also perform a convergence study for the solution extrapolated to the boundary, which is shown in Figure <ref>. We again see approximate second order rates for this reconstruction. A plot of the norms of the solution differences is shown in Figures <ref> and <ref>. The differences in the finest grid show a difference of less than one percent at the final time.

To assess conservation, we perform a convergence study on the total amount of the two concentration fields. Figure <ref> shows the total amounts of both fields as a function of time, while the sum of both amounts remains relatively constant. Despite this method not being conservative, we see a loss of less than one percent over the course of the finest simulation. The total amount of concentration converges at a rate that is between first and second order accuracy.

Concentrations $q_\text{o}$ (orange) and $q_\text{v}$ (blue) at times $0.0$ (fig:reacting_couette:0), $0.5$ (fig:reacting_couette:1), $1.0$ (fig:reacting_couette:2), $1.5$ (fig:reacting_couette:3), and $2.0$ (fig:reacting_couette:4). The zero contour of the level set $\phi_\text{o}$ is shown in green. The diffusion coefficient is fixed at $D = 0.05$ with $1024$ points in each direction. We use a fixed time-step corresponding to a maximum CFL number of $C_\text{CFL} = 0.5$ to a final time of $T = 2$.

Numerical convergence study for the interior of the disk in the $L^1$ (fig:reaction_in_converge:l1_rates), $L^2$ (fig:reaction_in_converge:l2_rates), and $L^\infty$ (fig:reaction_in_converge:max_rates) norms. The convergence rates are computed from simulations with $N = 128,256, \text{ and } 512$ points (blue), $N = 256,512,\text{ and } 1024$ points (orange), and $N = 512,1024, \text{ and } 2048$ points (green). Also shown are the convergence rates for the extrapolation of the solution to the boundary (fig:reaction_in_converge:bdry_rates).

Numerical convergence study for the exterior of the disk but inside the two cylinders in the $L^1$ (fig:reaction_out_converge:l1_rates), $L^2$ (fig:reaction_out_converge:l2_rates), and $L^\infty$ (fig:reaction_out_converge:max_rates) norms. The convergence study is computed from simulations with $N = 128,256,\text{ and }512$ points (blue), $N = 256,512, \text{ and } 1024$ points (orange), and $N = 512,1024,\text{ and } 2048$ points (green).

Norms of solution differences for the concentration field inside the domain for the $L^1$ (fig:reaction_in_norms:l1_norms), $L^2$ (fig:reaction_in_norms:l2_norms), and $L^\infty$ (fig:reaction_in_norms:max_norms) norms. Also shown are the norms of the differences between solutions extrapolated to the boundary (fig:reaction_in_norms:bdry_norms).

Norms of solution differences for the concentration field outside the domain for the $L^1$ (fig:reaction_out_norms:l1_norms), $L^2$ (fig:reaction_out_norms:l2_norms), and $L^\infty$ (fig:reaction_out_norms:max_norms) norms.

Total concentration of both $q_\text{v}$ and $q_\text{o}$ for interactions with a oscillatory Couette flow (fig:total_amount:0). The total amount converges at a rate that is between first and second order (fig:total_amount:1).

§ CONCLUSION

We have presented a numerical method to simulate advection-diffusion problems with Robin boundary conditions on irregular, evolving domains. The method shows second order convergence in the $L^1$, $L^2$, and $L^\infty$ norms. We have also demonstrated the ability to accurately reconstruct solution values on the boundary, achieving second order accurate results. We have further demonstrated the method to be able to capture transformation of one concentration into another concentration with all interaction mediated through boundary conditions. Although only two dimensional tests are presented, we expect this method to be easily extended to three spatial dimensions. The use of radial basis functions when compared to moving least squares gives optimal convergence rates and more accurate solutions for the grid spacings considered. Z-splines were used in this study for computational efficiency, although any suitable reconstruction procedure should work. This method has many possible applications, including models with chemical transport in evolving bodies, such as esophageal transport, oxygen flow in the lungs, and blood flow.

§ ACKNOWLEDGEMENTS

The authors thank Dr. Varun Shankar for helpful discussions on both radial basis functions and semi-Lagrangian schemes. This research was supported through NHLBI award 5U01HL143336 and NSF awards DMS 1664645, OAC 1450327, and OAC 1931516.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

journaltitleJournal of Computational Physics
titlePoisson equations in irregular domains with Robin boundary conditions - Solver with second-order accurate gradients

Argonne National Laboratory

numberANL-95/11 - Revision 3.13
titlePETSc Users Manual

titlePETSc Web page

Birkhäuser Press

booktitleModern Software Tools in Scientific Computing
titleEfficient Management of Parallelism in Object Oriented Numerical Software Libraries

journaltitleJournal of Scientific Computing
titleComparison of Moving Least Squares and RBF+poly for Interpolation and Derivative Approximation

abstractWe present a forward semi-Lagrangian numerical method for systems of transport equations able to advect smooth and discontinuous fields with high-order accuracy. The numerical scheme is composed of an integration of the transport equations along the trajectory of material elements in a moving grid and a reconstruction of the fields in a reference regular mesh using a nonlinear mapping and adaptive moment-preserving interpolations. The nonlinear mapping allows for the arbitrary deformation of material elements. Additionally, interpolations can represent discontinuous fields using adaptive-order interpolation near jumps detected with a slope-limiter function. Due to the large number of operations during the interpolations, a serial implementation of this scheme is computationally expensive. The scheme has been accelerated in many-core parallel architectures using a thread per grid node and parallel data gathers. We present a series of tests that show the scheme to be an attractive option for simulating advection equations in multidimensions with high accuracy.
journaltitleSIAM Journal on Scientific Computing
titleMoments preserving and high-resolution semi-lagrangian advection scheme

journaltitleDrug Discovery Today
titleModels for drug absorption from the small intestine: where are we and where are we going?

journaltitleJournal of Computational Physics
titleSolving Poisson-type equations with Robin boundary conditions on piecewise smooth interfaces

Berlin, Heidelberg

Springer Berlin Heidelberg

booktitleNumerical Mathematics and Advanced Applications 2009
titleAn Immersed Interface Technique for the Numerical Solution of the Heat Equation on a Moving Domain

Elsevier Inc.

abstractThis paper proposes a finite difference discretization method for simulations of heat and mass transfer with Robin boundary conditions on irregular domains. The level set method is utilized to implicitly capture the irregular evolving interface, and the ghost fluid method to address variable discontinuities on the interface. Special care has been devoted to providing ghost values that are restricted by the Robin boundary conditions. Specifically, it is done in two steps: 1) calculate the normal derivative in cells adjacent to the interface by reconstructing a linear polynomial system; 2) successively extrapolate the normal derivative and the ghost value in the normal direction using a linear partial differential equation approach. This method produces second-order accurate solutions for both the Poisson and heat equations with Robin boundary conditions, and first-order accurate solutions for the Stefan problems. The solution gradients are of first-order accuracy, as expected. It is easy to implement in three-dimensional configurations, and can be straightforwardly generalized into higher-order variants. The method thus represents a promising tool for practical heat and mass transfer problems involving Robin boundary conditions.
journaltitleJournal of Computational Physics
titleA finite difference discretization method for heat and mass transfer with Robin boundary conditions on irregular domains
Finite difference scheme,Ghost value,Irregular domain,Level set method,Robin boundary condition,Stefan problem

titleSemi-Lagrangian Methods

titleSemi-Lagrangian Advection Methods and Their Applications in Geoscience
Earth sciences – Mathematical models

journaltitleJournal of Computational Physics
titleOn the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy

abstractCopyright $\copyright$ 2015 by Annual Reviews. All rights reserved.Intravascular blood clots form in an environment in which hydrodynamic forces dominate and in which fluid-mediated transport is the primary means of moving material. The clotting system has evolved to exploit fluid dynamic mechanisms and to overcome fluid dynamic challenges to ensure that clots that preserve vascular integrity can form over the wide range of flow conditions found in the circulation. Fluid-mediated interactions between the many large deformable red blood cells and the few small rigid platelets lead to high platelet concentrations near vessel walls where platelets contribute to clotting. Receptor-ligand pairs with diverse kinetic and mechanical characteristics work synergistically to arrest rapidly flowing cells on an injured vessel. Variations in hydrodynamic stresses switch on and off the function of key clotting polymers. Protein transport to, from, and within a developing clot determines whether and how fast it grows. We review ongoing experimental and modeling research to understand these and related phenomena.
journaltitleAnnual Review of Fluid Mechanics
titleFluid Mechanics of Blood Clot Formation

journaltitleToxicological Sciences
titleA Validated Hybrid Computational Fluid Dynamics-Physiologically Based Pharmacokinetic Model for Respiratory Tract Vapor Absorption in the Human and Rat and Its Application to Inhalation Dosimetry of Diacetyl

abstractLike many problems in biofluid mechanics, cardiac mechanics can be modeled as the dynamic interaction of a viscous incompressible fluid (the blood) and a (visco-)elastic structure (the muscular walls and the valves of the heart). The immersed boundary method is a mathematical formulation and numerical approach to such problems that was originally introduced to study blood flow through heart valves, and extensions of this work have yielded a three-dimensional model of the heart and great vessels. In the present work, we introduce a new adaptive version of the immersed boundary method. This adaptive scheme employs the same hierarchical structured grid approach (but a different numerical scheme) as the two-dimensional adaptive immersed boundary method of Roma et al. [A multilevel self adaptive version of the immersed boundary method, Ph.D. Thesis, Courant Institute of Mathematical Sciences, New York University, 1996; An adaptive version of the immersed boundary method, J. Comput. Phys. 153 (2) (1999) 509-534] and is based on a formally second order accurate (i.e., second order accurate for problems with sufficiently smooth solutions) version of the immersed boundary method that we have recently described [B.E. Griffith, C.S. Peskin, On the order of accuracy of the immersed boundary method: higher order convergence rates for sufficiently smooth problems, J. Comput. Phys. 208 (1) (2005) 75-105]. Actual second order convergence rates are obtained for both the uniform and adaptive methods by considering the interaction of a viscous incompressible flow and an anisotropic incompressible viscoelastic shell. We also present initial results from the application of this methodology to the three-dimensional simulation of blood flow in the heart and great vessels. The results obtained by the adaptive method show good qualitative agreement with simulation results obtained by earlier non-adaptive versions of the method, but the flow in the vicinity of the model heart valves indicates that the new methodology provides enhanced boundary layer resolution. Differences are also observed in the flow about the mitral valve leaflets. ?? 2006 Elsevier Inc. All rights reserved.
journaltitleJournal of Computational Physics
titleAn adaptive, formally second order accurate version of the immersed boundary method

Adaptive mesh refinement,Blood flow,Cardiac mechanics,Convergence,Fluid-structure interaction,Hemodynamics,Immersed boundary method,Projection method

abstractFluid–structure interaction is ubiquitous in nature and occurs at all biological scales. Immersed methods provide mathematical and computational frameworks for modeling fluid–structure systems. These methods, which typically use an Eulerian description of the fluid and a Lagrangian description of the structure, can treat thin immersed boundaries and volumetric bodies, and they can model structures that are flexible or rigid or that move with prescribed deformational kinematics. Immersed formulations do not require body-fitted discretizations and thereby avoid the frequent grid regeneration that can otherwise be required for models involving large deformations and displacements. This article reviews immersed methods for both elastic structures and structures with prescribed kinematics. It considers formulations using integral operators to connect the Eulerian and Lagrangian frames and methods that directly apply jump conditions along fluid–structure interfaces. Benchmark problems demonstrate the effectiveness of these methods, and selected applications at Reynolds numbers up to approximately 20,000 highlight their impact in biological and biomedical modeling and simulation.
journaltitleAnnual Review of Fluid Mechanics
titleImmersed Methods for Fluid-Structure Interaction

abstractWe consider the Poisson equation with mixed Dirichlet, Neumann and Robin boundary conditions on irregular domains. We describe a straightforward and efficient approach for imposing the mixed boundary conditions using a hybrid finite-volume/finite-difference approach, leveraging on the work of Gibou et al. (2002) [14], Ng et al. (2009) [30] and Papac et al. (2010) [33]. We utilize three different level set functions to represent the irregular boundary at which each of the three different boundary conditions must be imposed; as a consequence, this approach can be applied to moving boundaries. The method is straightforward to implement, produces a symmetric positive definite linear system and second-order accurate solutions in the L221e-norm in two and three spatial dimensions. Numerical examples illustrate the second-order accuracy and the robustness of the method.
journaltitleComputers and Fluids
titleImposing mixed Dirichlet-Neumann-Robin boundary conditions in a level-set framework
Finite difference,Level set,Mixed boundary conditions

abstractA major challenge facing software libraries for scientific computing is the ability to provide adequate flexibility to meet sophisticated, diverse, and evolving application requirements. Object-oriented design techniques are valuable tools for capturing characteristics of complex applications in a software architecture. In this paper, we describe certain prominent object-oriented features of the SAMRAI software library that have proven to be useful in application development. SAMRAI is used in a variety of applications and has demonstrated a substantial amount of code and design re-use in those applications. This flexibility and extensibility is illustrated with three different application codes. We emphasize two important features of our design. First, we describe the composition of complex numerical algorithms from smaller components which are usable in different applications. Second, we discuss the extension of existing framework components to satisfy new application needs.
journaltitleConcurrency and Computation: Practice and Experience
titleManaging application complexity in the SAMRAI object-oriented framework
Adaptive mesh refinement,Design patterns,Object-oriented programming
titleIBAMR: An adaptive and distributed-memory parallel implementation of the immersed boundary method

journaltitleNumerical Methods for Partial Differential Equations
titleConservative semi-Lagrangian advection on adaptive unstructured meshes

journaltitleJournal of Computational Physics
titleA volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles

abstractCut-cell meshes present an attractive alternative to terrain-following coordinates for the representation of topography within atmospheric flow simulations, particularly in regions of steep topographic gradients. in this paper, we present an explicit two-dimensional method for the numerical solution on such meshes of atmospheric flow equations including gravitational sources. This method is fully conservative and allows for time steps determined by the regular grid spacing, avoiding potential stability issues due to arbitrarily small boundary cells. We believe that the scheme is unique in that it is developed within a dimensionally split framework, in which each coordinate direction in the flow is solved independently at each time step. Other notable features of the scheme are: (i) its conceptual and practical simplicity, (ii) its flexibility with regard to the one-dimensional flux approximation scheme employed, and (iii) the well-balancing of the gravitational sources allowing for stable simulation of near-hydrostatic flows. The presented method is applied to a selection of test problems including buoyant bubble rise interacting with geometry and lee-wave generation due to topography. This journal is $\copyright$ 2009 The Royal Society.
journaltitlePhilosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences
titleWell-balanced compressible cut-cell simulation of atmospheric flow
Atmospheric,Conservative,Cut cell,Dimensional splitting,Numerical methods,Well balancing

journaltitleJournal of Computational Physics
titleAn immersed interface method for discrete surfaces

abstractThe body's response to vascular injury involves two intertwined processes: platelet aggregation and coagulation. Platelet aggregation is a predominantly physical process, whereby platelets clump together, and coagulation is a cascade of biochemical enzyme reactions. Thrombin, the major product of coagulation, directly couples the biochemical system to platelet aggregation by activating platelets and by cleaving fibrinogen into fibrin monomers that polymerize to form a mesh that stabilizes platelet aggregates. Together, the fibrin mesh and the platelet aggregates comprise a thrombus that can grow to occlusive diameters. Transport of coagulation proteins and platelets to and from an injury is controlled largely by the dynamics of the blood flow. To explore how blood flow affects the growth of thrombi and how the growing masses, in turn, feed back and affect the flow, we have developed the first spatial-temporal mathematical model of platelet aggregation and blood coagulation under flow that includes detailed descriptions of coagulation biochemistry, chemical activation and deposition of blood platelets, as well as the two-way interaction between the fluid dynamics and the growing platelet mass. We present this model and use it to explain what underlies the threshold behaviour of the coagulation system's production of thrombin and to show how wall shear rate and near-wall enhanced platelet concentrations affect the development of growing thrombi. By accounting for the porous nature of the thrombus, we also demonstrate how advective and diffusive transport to and within the thrombus affects its growth at different stages and spatial locations. $\copyright$ The author 2010. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
journaltitleMathematical Medicine and Biology
titleGrow with the flow: A spatial-temporal model of platelet deposition and blood coagulation under flow
Mathematical model,Thrombosis

journaltitleJournal of Computational Physics
titleAn unconditionally stable fully conservative semi-Lagrangian method

journaltitleCommunications in Mathematical Sciences
titleSolving pdes in complex geometries

journaltitleJournal of Computational Physics
titleA robust incompressible Navier-Stokes solver for high density ratio multiphase flows

abstractWe describe the implementation of an interpolation technique, which allows the accurate imposition of the Dirichlet, Neumann, and mixed (Robin) boundary conditions on complex geometries using the immersed-boundary technique on Cartesian grids, where the interface effects are transmitted through forcing functions. The scheme is general in that it does not involve any special treatment to handle either one of the three types of boundary conditions. The accuracy of the interpolation algorithm on the boundary is assessed using several two- and three-dimensional heat transfer problems: (1) forced convection over cylinders placed in an unbounded flow, (2) natural convection on a cylinder placed inside a cavity, (3) heat diffusion inside an annulus, and (4) forced convection around a stationary sphere. The results show that the scheme preserves the second-order accuracy of the equations solver and are in agreement with analytical and/or numerical data.
journaltitleJournal of Heat Transfer
titleA General Scheme for the Boundary Conditions in Convective and Diffusive Heat Transfer With Immersed Boundary Methods

journaltitleDrug Metabolism and Disposition
titleModeling of Intestinal Drug Absorption: Roles of Transporters and Metabolic Enzymes (for the Gillette Review Series)

journaltitleJournal of Computational Physics
titleEfficient symmetric discretization for the Poisson, heat and Stefan-type problems with Robin boundary conditions

Cambridge University Press

abstractThis paper is concerned with the mathematical structure of the immersed boundary (IB) method, which is intended for the computer simulation of fluid-structure interaction, especially in biological fluid dynamics. The IB formulation of such problems, derived here from the principle of least action, involves both Eulerian and Lagrangian variables, linked by the Dirac delta function. Spatial discretization of the IB equations is based on a fixed Cartesian mesh for the Eulerian variables, and a moving curvilinear mesh for the Lagrangian variables. The two types of variables are linked by interaction equations that involve a smoothed approximation to the Dirac delta function. Eulerian/Lagrangian identities govern the transfer of data from one mesh to the other. Temporal discretization is by a second-order Runge–Kutta method. Current and future research directions are pointed out, and applications of the IB method are briefly discussed.
journaltitleActa Numerica
titleThe immersed boundary method

journaltitleJournal of Computational Physics
titleThe smooth forcing extension method: A high-order technique for solving elliptic equations on complex domains

journaltitleKinetic & Related Models
titleSemilagrangian schemes applied to moving boundary problems for the BGK model of rarefied gas dynamics

abstractA cut-cell method for Cartesian meshes to simulate viscous compressible flows with moving boundaries is presented. We focus on eliminating unphysical oscillations occurring in Cartesian grid methods extended to moving-boundary problems. In these methods, cells either lie completely in the fluid or solid region or are intersected by the boundary. For the latter cells, the time dependent volume fraction lying in the fluid region can be so small that explicit time-integration schemes become unstable and a special treatment of these cells is necessary. When the boundary moves, a fluid cell may become a cut cell or a solid cell may become a small cell at the next time level. This causes an abrupt change in the discretization operator and a suddenly modified truncation error of the numerical scheme. This temporally discontinuous alteration is shown to act like an unphysical source term, which deteriorates the numerical solution, i.e., it generates unphysical oscillations in the hydrodynamic forces exerted on the moving boundary. We develop an accurate moving boundary formulation based on the varying discretization operators yielding a cut-cell method which avoids these discontinuities. Results for canonical two- and three-dimensional test cases evidence the accuracy and robustness of the newly developed scheme. $\copyright$ 2012 Elsevier Inc.
journaltitleJournal of Computational Physics
titleAn accurate moving boundary formulation in cut-cell methods
Compressible Navier-Stokes,Cut-cell method,Immersed boundary method,Moving boundary,Oscillating cylinder,Spurious force oscillations

journaltitleJournal of Aerosol Medicine
titleAnalysis of particle deposition in the turbinate and olfactory regions using a human nasal computational fluid dynamics model

journaltitleJournal of Computational Physics
titleImmersed Boundary Smooth Extension (IBSE): A high-order method for solving incompressible flows in arbitrary smooth domains

journaltitleSIAM Journal on Numerical Analysis
titleOn the Construction and Comparison of Difference Schemes

journaltitleCommunications in Applied Mathematics and Computational Science
titleA cut-cell method for simulating spatial models of biochemical reaction networks in arbitrary geometries

journaltitleJournal of Biomechanical Engineering
titleDevelopment of a CFD boundary condition to model transient vapor absorption in the respiratory airways

journaltitleJournal of Computational Physics
titleA source term method for Poisson problems on irregular domains

journaltitleJournal of Computational Physics
titleAn immersed interface method for simulating the interaction of a fluid with moving boundaries

journaltitleInternational Journal for Numerical Methods in Fluids
titleSimulations of chemical transport and reaction in a suspension of cells I: an augmented forcing point method for the stationary case

abstractIn this paper, we describe an approach for solving partial differential equations with general boundary conditions imposed on arbitrarily shaped boundaries. A continuous function, the domain parameter, is used to modify the original differential equations such that the equations are solved in the region where a domain parameter takes a specified value while boundary conditions are imposed on the region where the value of the domain parameter varies smoothly across a short distance. The mathematical derivations are straightforward and applicable to a wide variety of partial differential equations. To demonstrate the general applicability of the approach, we provide four examples herein: (1) the diffusion equation with both Neumann and Dirichlet boundary conditions; (2) the diffusion equation with both surface diffusion and reaction; (3) the mechanical equilibrium equation; and (4) the equation for phase transformation with the presence of additional boundaries. The solutions for several of these cases are validated against numerical solutions of the corresponding sharp-interface equations. The potential of the approach is demonstrated with five applications: surface-reaction-diffusion kinetics with a complex geometry, Kirkendall-effect-induced deformation, thermal stress in a complex geometry, phase transformations affected by substrate surfaces and relaxation of a droplet on irregular surfaces. $\copyright$ 2012 IOP Publishing Ltd.
journaltitleModelling and Simulation in Materials Science and Engineering
titleExtended smoothed boundary method for solving partial differential equations with general boundary conditions on complex boundaries

journaltitleAdvanced Drug Delivery Reviews
titleTransport approaches to the biopharmaceutical design of oral drug delivery systems: prediction of intestinal absorption

Elsevier Inc.

abstractAn immersed-boundary method, previously developed for flow-field simulation, is extended to study heat-transfer problems of flow over a circular cylinder. The Dirichlet-type isothermal temperature boundary conditions and the Neumann-type iso-heat-flux boundary conditions are both implemented. To verify the accuracy of the simulation method, L2-norm are computed to test the order of accuracy of the scheme. Numerical solutions are then further validated by comparing simulated temperature distributions and local convective-heat-transfer coefficients in flow over a stationary circular cylinder with data in the literature. Comparisons are made at different Reynolds-number flows and under different types of temperature boundary conditions. Finally, the effect of heat-transfer enhancement for flow over a transversely oscillating cylinder is investigated.
journaltitleInternational Journal of Heat and Fluid Flow
titleStudy of heat-transfer on the surface of a circular cylinder in flow using an immersed-boundary method
Flow over a cylinder,Heat convection,Immersed-boundary method