# Deriving RIP sensing matrices for sparsifying dictionaries

Jinn Ho and Wen-Liang Hwang

###### Abstract

Compressive sensing involves the inversion of a mapping
$SD\in\mathbb{R}^{m\times n}$, where $m<n$, $S$ is a sensing matrix, and $D$
is a sparisfying dictionary. The restricted isometry property is a powerful
sufficient condition for the inversion that guarantees the recovery of high-
dimensional sparse vectors from their low-dimensional embedding into a
Euclidean space via convex optimization. However, determining whether $SD$ has
the restricted isometry property for a given sparisfying dictionary is an NP-
hard problem, hampering the application of compressive sensing. This paper
provides a novel approach to resolving this problem. We demonstrate that it is
possible to derive a sensing matrix for any sparsifying dictionary with a high
probability of retaining the restricted isometry property. In numerical
experiments with sensing matrices for K-SVD, Parseval K-SVD, and wavelets, our
recovery performance was comparable to that of benchmarks obtained using
Gaussian and Bernoulli random sensing matrices for sparse vectors.

##  1 Introduction

The Johnson-Lindenstrauss lemma [1, 2] states that a set of points in a high-
dimensional space can be embedded into a space of far lower dimensionality
using an invertible Lipschitz function, such that distances between points are
nearly preserved. Compressive sensing (CS) links this result within the
framework of sparse representation. The problem of CS inversion involves the
recovery of sparse vectors in a high-dimensional Euclidean space from their
low-dimensional embedding via $SD$, where $S$ is a sensing matrix and $D$ a
sparsifying dictionary. In accordance with the amount of a priori knowledge
that is available, the issue of recovery has been formulated as a variety of
problems in which either both $S$ and $D$ are given (the CS recovery problem)
or neither is given (simultaneously learning of sensing matrix and sparsifying
dictionary [3]).

The restricted isometry property (RIP)[4] of order $k$ assumes a restricted
isometry constant $\delta_{k}\in[0,1)$, such that for any $k$-sparse vector
$x\in\mathbb{R}^{n}$,

$\displaystyle(1-\delta_{k})\|x\|_{2}^{2}\leq\|SDx\|_{2}^{2}\leq(1+\delta_{k})\|x\|_{2}^{2}.$

If $\delta_{k}$ is small, then $SD$ is nearly linearly isometric, which means
that it preserves the distance between any pair of $\frac{k}{2}$-sparse
vectors. Candes [5] showed that k-sparse vectors can be correctly recovered
via CS using convex optimization under the assumption that $SD$ has RIP with
restricted isometry constant $\delta_{2k}$, where $k$ is the sparsity level
when $\delta_{2k}<\sqrt{2}-1$. The latter bound can be improved to $0.307$
[6]. It has been demonstrated in [7] that a stable inversion can be achieved
for both $k$-sparse and compressible signals when $D$ is a basis (invertible
square matrix) and that $SD$ has the RIP for $3k$ sparse vectors. The main
issue to utilize the implications of the RIP is thus in constructing a sensing
matrix $S$ for a given sparsifying dictionary $D$ such that $SD$ has the RIP.

There are sensing matrices that satisfy with high probability the RIP for a
particular class of bases. In [7, 8, 9], the authors reported that a random
matrix $S\in\mathbb{R}^{m\times n}$ with independently identically distributed
(i.i.d.) entries (e.g., Gaussian and Bernoulli distributed) that satisfy the
concentration inequality is universal for any orthonormal matrix
$D\in\mathbb{R}^{n\times n}$, in the sense that there is a high probability
that $SD$ retains the RIP when $m\geq ck\log(\frac{n}{k})$, where $k$ is the
sparsity level of vectors in $\mathbb{R}^{n}$ to be recovered. Thus, a signal
that is sparse with respect to some orthonormal basis can be coupled with
random matrices to achieve the RIP.

In practice however, one more commonly encounters signals that are sparse with
respect to a redundant dictionary that may be far from orthonormal (e.g.,
Gabor-, curvelet-, wavelet-, or other data-driven learned dictionaries). One
solution to the above problem with respect to tight frames employs the notion
of D-RIP [10] 111Let $\Sigma_{k}$ be the the union of all subspaces spanned by
all subsets of $k$ columns of $D$. The sensing matrix $A$ obeys the RIP
adapted to $D$ with constant $\delta_{k}$ if
$(1-\delta_{k})\|x\|_{2}^{2}\leq\|Ax\|_{2}^{2}\leq(1+\delta_{k})\|x\|_{2}^{2}$
holds for all $x\in\Sigma_{k}$.. Assuming that the tight-frame coefficients
$D^{\top}x$ are sparse, CS-recovery of $x$ from $y=Ax+\eta$, where
$A\in\mathbb{R}^{m\times n}$ is a random Gaussian sensing matrix and the
energy of the additive noise $\eta$ is bounded by $\|\eta\|_{2}<\epsilon$, can
be achieved by solving

$\displaystyle\begin{cases}\min_{x}\|D^{\top}x\|_{1}\\\
\|y-Ax\|_{2}\leq\epsilon.\end{cases}$ (1)

It has been demonstrated in [10] that under the D-RIP framework, with high
probability, it is possible to obtain stable inversion of $k$-sparse or
compressible signals by solving (1) in cases where the number of measurements
$m$ is on the order of $k\log\frac{n}{k}$.

The nullspace property gives necessary and sufficient conditions on the
reconstruction of sparse signals using the convex optimization [11]. The
nullspace property is often difficult to check in practice. Mutual incoherence
between the sensing matrix $S$ and the dictionary $D$ is another frequently
encountered condition that facilitates CS-recovery, albeit with greater
restrictions on the sparsity level of the vectors to be recovered than for the
RIP [12, 8]. In a given application, the prevalence of one sufficient
condition over others is generally determined by a trade-off between the cost,
memory capacity, and computational capability of the relevant devices.

This paper follows the line of the RIP in deriving a novel approach to solving
the following problem: Derive a sensing matrix $S$ for a given spasifying
dictionary $D$, such that with high probability $SD$ retains the RIP. It is
widely believed that in seeking to satisfy the RIP, deriving random sensing
matrices uses less number of sensing measurements than deriving deterministic
sensing matrices [13]. Throughout the remainder of this paper, the symbol
$\mathcal{E}$ is used to indicate the random-row-selection operation without
repetition. We base our approach on the assumption that the dictionary
$D\in\mathbb{R}^{l\times n}$, with $l\leq n$, takes the form of $GAH$, where
$G\in\mathbb{R}^{l\times l}$ is invertible, $H\in\mathbb{R}^{n\times n}$ is
orthonormal, and the probability distribution associated with the generation
of $\mathcal{E}A\in\mathbb{R}^{m\times n}$ satisfies the concentration
inequality of [2, 14, 15], see Lemma 1. The embedding function takes the usual
form of $SD$, where $S$ refers to the sensing matrix. Clearly, if
$\mathcal{E}G^{-1}$ is the sensing matrix of $D$, then $SD=\mathcal{E}AH$ has
the RIP in the sense described in [16].

Clearly, the success of the proposed approach depends on whether $D$ can be
expressed as a $GAH$. Thus, the approach is referred to as the compressive
sensing factorization approach to highlight the fact that a sparsifying
dictionary is factorized into components that are used for compressing sensing
purpose. We demonstrate that if the matrix $A$ (of the same rank of $D$) is
also given, the answer to the problem is affirmative. In other words, given
$D\in\mathbb{R}^{l\times n}$ and $A\in\mathbb{R}^{l\times n}$, there exists an
invertible matrix $G$ and orthonormal matrix $H$, such that $D=GAH$, provided
that $D$ and $A$ have equal rank. If $l<n$, then the solutions $G$ and $H$ are
not unique, due to the existence of a non-zero null space spanned by the rows
of $A$. Despite their non-uniqueness, we demonstrate that both $G$ and $H$ can
be determined analytically. In numerical experiments for the sensing matrices
for K-SVD [17], Parseval K-SVD[18], and CDF 9-7 wavelets [19] we verify that
the recovery performance of our approach is comparable to that of benchmarks
achieved using Gaussian or Bernoulli random matrices as sensing matrices for
sparse vectors.

The remainder of the paper is organized as follows. Section 2 presents the
background of CS. Section 3 presents the technical underpinnings of our
approach and outlines the method by which we derive the sensing matrix for a
given dictionary. Section 4 presents experimental results on acquisition and
representation using real-world images. Concluding remarks are presented in
Section 5.

##  2 Background

The Johnson-Lindenstrauss lemma is concerned with embedding a set of $N$
points in $\mathbb{R}^{n}$ within a lower dimensional $\mathbb{R}^{m}$ (with
$m$ is as small as possible), aiming to approximately preserve the distances
between any two of the $N$ points.222Let $\epsilon\in(0,1)$. For every set of
$N$ points in $\mathbb{R}^{n}$ and $m>{\cal O}(\frac{\ln N}{\epsilon^{2}})$,
there exists a Lipschitz mapping $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$,
such that for any two points $u$ and $v$ in the set,
$(1-\epsilon)\|u-v\|^{2}_{2}\leq\|f(u)-f(v)\|^{2}_{2}\leq(1+\epsilon)\|u-v\|^{2}_{2}.$
Various statements have been made based on this lemma, and a concise
description of its evolution can be found in [20, 21]. CS is linked to the
Johnson-Lindenstrauss lemma in the framework of sparse representation, where
the Lipschitz embedding is $SD$.

Sparse representations [22, 8, 7] seek to determine the sparsest coefficient
vector (in terms of $\ell_{0}$-norm) from which signals can be synthesized
using an over-complete dictionary. The $\ell_{0}$-norm objective is usually
relaxed to the $\ell_{1}$-norm, in order to convert the non-convex into a
convex optimization problem. The problem has gradually evolved from the
recovery of sparse vectors based on observations, to deriving a sparsifying
dictionary based on a set of observation [17, 23, 24]. In [25], learning an
over-complete sparsifying dictionary after pre-conditioning empirical
observations was reduced to learning an orthogonal transform, based on which
the optimality of the solution and the complexity of the learning problem can
theoretically be guaranteed. In [26], the parameters (i.e., the size of the
dictionary and sparsity level) and the sparsifying dictionary are learned
adaptively. There are numerous situations in which the adopted dictionary
lacks the precision required for sparse representation. Pioneering work on the
sensitivity of an imprecise dictionary to CS recovery was performed in [27].
The sensitivity of the sparse recovery against dictionary perturbation in
terms of the restricted isometry constant under the CS paradigm is outlined in
[28].

When using CS for signal acquisition, signals that have sparse coefficients
with respect to $D$ can be mapped into a much lower-dimensional space via a
sensing matrix $S$. Despite the results in [16] \- indicating that
orthonormal, Gaussian and Bernoulli sensing matrices retain the RIP - deriving
the composition of a dictionary with a sensing matrix that retains the RIP is
an NP-hard problem and the main issue hindering the imposition of the RIP as a
constraint in CS problems. Nonetheless, researchers have devised ingenious
dictionary and sensing matrix pairs as well as algorithms to facilitate CS
recovery. It has been demonstrated that if the sparsity level $k$ of a vector
$x$ is smaller than a constant multiple of the reciprocal of the coherence of
$SD$, then it is possible to recover $x$ from $y=SDx$ [29]. Moreover,
numerical experiments have demonstrated that CS works well even when the
sparsity level exceeds the theoretically guaranteed bound. For example, Elad
[30] demonstrated that the optimal choice of sensing matrix can improve
incoherence with a pre-defined sparsifying dictionary, thereby substantially
improving CS reconstruction performance. A number of algorithms have been
developed with the aim of recovering both the sensing matrix and the
sparsifying dictionary from observations. For example, the method in [3] can
simultaneously learn the sensing matrix and sparsifying dictionary from given
samples $\\{(y_{i},z_{i})\\}_{i=1}^{N}$, for which $y_{i}=Dx_{i}$ and
$z_{i}=Sy_{i}+\eta$, where $\eta$ is additive noise contained in the sensing
system.

From the perspectives of computation and implementation, it is desirable to
use structured sensing matrices [31]. In [13], random-row-sampling
$\mathcal{E}$ is used as the sensing matrix for any bounded orthonormal
system. A well-known example of this approach is the random partial Fourier
matrix, in which rows of a Fourier matrix are randomly sampled without
repetition. Random partial circular and Toeplitz matrices, obtained by
applying $\mathcal{E}$ as a sensing matrix to circular and Toeplitz matrices,
are used in wireless communication and radar applications [32]. Romberg [33]
developed an approach in which $\mathcal{E}$ is used as the sensing matrix for
a random square matrix of the form $D=(F^{*}\Sigma
F)\psi\in\mathbb{R}^{n\times n}$, where $\psi$ is an orthonormal system, the
entries of the diagonal matrix $\Sigma$ are randomly selected numbers within
specified ranges, and $F$ is the Fourier matrix. The RIPless approach [34]
employs a sensing mechanism capable of approximately recovering sparse signals
from noisy measurements, as long as the random sensing measurements are
selected independently from a probability distribution with specific
properties.

The fact that Gaussian and Bernoulli sensing matrices do not provide rapid
matrix multiplications or economical storage precludes their adoption as
sensing matrices for large scale orthonormal systems. To the best of our
knowledge, matrix factorization was first used for a sensing matrix in the
structurally random method (SRM) [35], where a matrix product involving
$\mathcal{E}$ is used as sensing matrix. The SRM makes it possible to
construct a sensing matrix that can be calculated rapidly, thereby eliminating
the need to store the sensing matrix explicitly. This approach is applicable
to any sensing matrix of the form $\mathcal{E}WR$, where $W$ is an orthonormal
matrix and $R$ is either a random permutation or a diagonal matrix. The low
mutual coherence of SRM sensing matrices with any orthonormal matrix makes
them suitable sensing matrices for large scale orthonormal systems.

The conventional approach to solving an inverse problem is via algorithms that
iteratively approximate a solution to within a desired error bound, in
conjunction with a theoretical analysis to determine the rate of convergence.
This approach can hamper applications that require rapid decisions, or
expensive computations for each iteration. Deep neural networks (DNNs) have
been used to solve inverse problems [36, 37] as well as the CS recovery
problem, wherein a sparse vector, or a signal that allows for sparse transform
coefficient, is derived via forward inference from a low-dimensional input.
Numerical experiments have demonstrated that the DNN approach is well-suited
to the CS recovery problem [38, 39, 40], allowing the recovery of a sparse
vector using far fewer measurements than the conventional approach.
Unfortunately, the DNN approach, so far, does not provide theoretical
justification for the solution.

##  3 Sensing matrix construction for prescribed dictionaries

In this section we provide the details of our proposed approach. Recall that
$\mathcal{E}$ denotes the random-row-selection operation without repetition.
In the following, let $A\in\mathbb{R}^{l\times n}$ be a random matrix, and
suppose $\mathcal{E}A\in\mathbb{R}^{m\times n}$ satisfies the concentration
inequality of the following result of [16], which guarantees that with high
probability $\mathcal{E}A$ has the RIP.

######  Lemma 1 (Concentration inequality: [16], Theorem 5.2).

Suppose $\mathcal{E}\in\mathbb{R}^{m\times l}$ and $A\in\mathbb{R}^{l\times
n}$. For any $x$, if the probability distribution generating
$\mathcal{E}A\in\mathbb{R}^{m\times n}$ satisfies
$\mathbb{E}\|\mathcal{E}Ax\|_{2}^{2}=\|x\|_{2}^{2}$ and for some
$\epsilon\in(0,1)$ the concentration inequality
$\mathbb{P}(|\|\mathcal{E}Ax\|_{2}^{2}-\|x\|_{2}^{2}|\geq\epsilon\|x\|_{2}^{2})\leq
2e^{-lc(\epsilon)}$ holds for some $c(\epsilon)>0$ depending only on
$\epsilon$, where the probability is taken over all $A\in\mathbb{R}^{l\times
n}$, then there exist $c_{1},c_{2}>0$ (depending only on $\delta_{k}$) such
that the RIP, with the restricted isometry constant $\delta_{k}$ and
$k\leq\frac{c_{1}m}{\log(n/k)}$, holds for $\mathcal{E}A$ with probability at
least $1-2e^{-c_{2}m}$.

The random matrices whose entries are independent realizations of Gaussian
random variables or Bernoulli random variables satisfy the above lemma [16].
It is clear that, if $\mathcal{E}A$ satisfies the concentration inequality,
then so does $\mathcal{E}AH$ for any orthonormal $H$. The following corollary
is therefore immediate.

###### Corollary 2 .

Suppose that $\mathcal{E}A$ satisfies the assumptions of Lemma 1 and that the
sparsifying dictionary $D\in\mathbb{R}^{l\times n}$ can be decomposed as $GAH$
for some invertible $G$ and orthonormal $H$. Then for the sensing matrix
$S:=\mathcal{E}G^{-1}\in\mathbb{R}^{m\times l}$, the embedding
$SD=\mathcal{E}AH$ has the RIP in the sense of Lemma 1.

By the following lemma, the existence of a factorization $D=GAH$ with
invertible $G$ and orthonormal $H$, as assumed in Corollary 2, is equivalent
to $D$ and $A$ having equal rank.

###### Lemma 3 .

Let $A,D\in\mathbb{R}^{l\times n}$ with $l\leq n$. Then, the following
statements are equivalent: 333(o) $\\{WD\mid W\in\mathbb{R}^{l\times l}\text{
is invertible}\\}\cap\\{AH\mid H\in\mathbb{R}^{n\times n}\text{ is
orthonormal}\\}\neq\emptyset$.

  * (i)

There exists an invertible $W\in\mathbb{R}^{l\times l}$ and orthonormal
$H\in\mathbb{R}^{n\times n}$ such that $WD=AH$.

  * (ii)

There exists an invertible $W\in\mathbb{R}^{l\times l}$ such that
$AA^{\top}=WDD^{\top}W^{\top}$.

  * (iii)

$A$ and $D$ have equal rank.

###### Proof.

The implications (i)$\Rightarrow$(ii)$\Rightarrow$(iii) are straightforward.
To show (iii)$\Rightarrow$(ii), assume $A$ and $D$ have equal rank $k$. Since
row-rank equals column-rank, also the symmetric matrices $AA^{\top}$ and
$DD^{\top}$ have equal rank $k$. Thus the desired $W$ can be derived from the
spectral decompositions $Q_{A}\Sigma_{A}Q_{A}^{\top}$ of $AA^{\top}$ and
$Q_{D}\Sigma_{D}Q_{D}^{\top}$ of $DD^{\top}$, where $\Sigma_{A}$ and
$\Sigma_{D}$ are diagonal with $k$ nonzero entries, assumed to be in the same
positions. Since $\text{ rank}(AA^{\top})=\text{ rank}(DD^{\top})$, there
exists an invertible diagonal matrix $\Sigma_{S}$ such that

$\Sigma_{A}=\Sigma_{S}\Sigma_{D}\Sigma_{S}.$ (2)

Then $W:=Q_{A}\Sigma_{S}Q_{D}^{\top}$ is invertible, and

$\displaystyle WDD^{\top}W^{\top}$
$\displaystyle=WQ_{D}\Sigma_{D}Q_{D}^{\top}W^{\top}=Q_{A}\Sigma_{S}\Sigma_{D}\Sigma_{S}Q_{A}^{\top}=AA^{\top}.$
(3)

To show (ii)$\Rightarrow$(i), by the equivalence of (ii) and (iii) we can
assume $A$ and $D$ have equal rank $k\leq l$, and considering the spectral
decompositions and $W$ as above. We let $A^{+}$ to denote the Moore–Penrose
pseudo-inverse of matrix $A$. Let $N_{A},N_{D}\in\mathbb{R}^{n\times(n-k)}$
consist of pairwise orthonormal columns spanning the null spaces of $A$, resp.
$D$, i.e., $AN_{A}=DN_{D}=0_{l\times(n-k)}$ and
$N_{A}^{\top}N_{A}=N_{D}^{\top}N_{D}=I_{n-k}$. Let
$H:=A^{+}WD+N_{A}N_{D}^{\top}$. Then

$AH=AA^{+}WD+AN_{A}N_{D}^{\top}=AA^{+}WD=WD$

since $AA^{+}$ is the identity on the range of $A$, of which the range of $W$
is a subspace, by definition. Moreover, by (3),

$HH^{T}=A^{+}WDD^{\top}W^{\top}(A^{+})^{\top}+N_{A}N_{A}^{\top}=A^{+}(AA^{\top}(A^{+})^{\top})+N_{A}N_{A}^{\top}=A^{+}A+N_{A}N_{A}^{\top}=I_{n}.\qed$

Note that the factorizations in Lemma 3 are not unique, since the spectral
decompositions involved in their construction depend on the order in which the
spectra are considered. Letting $G:=W^{-1}$, Lemma 3 implies the following
result.

###### Theorem 4 .

Suppose $A,D\in\mathbb{R}^{l\times n}$ ($l\leq n$) have equal rank $k$. Then
$D=GAH$ for some invertible $G$ and orthonormal $H$ (and this decomposition is
not unique).

A sensing matrix $S$ for a given sparsifying dictionary $D$ can be constructed
from any $A$ of equal rank that satisfies the concentration inequality, and
any factorization $D=GAH$ as in Theorem 4. Indeed, letting $S:=DG^{-1}$,
Corollary 2 then implies that the operator $SD=DAH$ has the RIP.

Tight frames are dictionaries $D$ with the property that $DD^{\top}$ is a
multiple of the identity. They are ubiquitous in signal and image processing
[41], since then the condition number of $DD^{\top}$ is equal to $1$ and the
canonical dual (used to derive decomposition coefficients of a signal) is
given by $D^{\top}$ (up to a constant rescaling). If $D$ is a tight frame, we
may assume $DD^{\top}=I_{l}$ (by rescaling the columns of $D$ if necessary).
If $A$ and $D$ have full rank, a factorization as in Theorem 4 can then be
constructed as follows.444If $A,D$ have equal – but not full – rank, the
following result holds replacing the inverses by the pseudo-inverses.

###### Theorem 5 .

Suppose $A,D\in\mathbb{R}^{l\times n}$ have full rank $l\leq n$, and $D$ is a
tight frame. Then $D=GAH$ with invertible, respectively orthonormal, factors
defined by

$\displaystyle G:=O(AA^{\top})^{-1/2}\quad\text{ and }\quad
H:=A^{\top}G^{\top}D+N_{A}N_{D}^{\top},$

where $O\in\mathbb{R}^{l\times l}$ is an orthonormal matrix and the columns of
$N_{A},N_{D}\in\mathbb{R}^{n\times(n-l)}$ are any orthonormal bases of the
null spaces of $A$, resp. $D$.

###### Proof.

Note first that $GAA^{\top}G^{\top}=I_{l}$ and thus $D=GAH$. Further,
$AN_{A}=DN_{D}=0_{l\times(n-l)}$ and
$N_{A}^{\top}N_{A}=N_{D}^{\top}N_{D}=I_{n-l}$. Thus the tight-frame property
of $D$ implies that

$\displaystyle HH^{\top}$
$\displaystyle=A^{\top}G^{\top}DD^{\top}GA+N_{A}N_{A}^{\top}=A^{\top}G^{\top}GA+N_{A}N_{A}^{\top}=A^{\top}(AA^{\top})^{-1}A+N_{A}N_{A}^{\top}=I_{n}.\qed$

We conclude this section by detailing one particular factorization according
to Theorem 4; based only on orthonormal bases for the ranges of $A$ and $D$:
Given $A,D\in\mathbb{R}^{l\times n}$ of rank $k\leq l$, choose
$U,V\in\mathbb{R}^{n\times k}$, each with pairwise orthonormal columns, such
that

$\displaystyle AUU^{\top}=A\quad\text{ and }\quad DVV^{\top}=D.$ (4)

(Such $U$ and $V$ may be derived via Gram-Schmidt through constructing
orthonormal column vectors spanning the ranges of $A$ and $D$.) Next consider
$U_{\perp},V_{\perp}\in\mathbb{R}^{n\times(n-k)}$ such that the extended
matrices $[U\,\,U_{\perp}],[V\,\,V_{\perp}]\in\mathbb{R}^{n\times n}$ are
orthonormal, and define

$\displaystyle H:=UV^{\top}+U_{\perp}V_{\perp}^{\top}.$ (5)

Then, clearly, $H$ is orthonormal and $AH=AUV^{\top}$. The latter implies, by
(4),

$\displaystyle D-GAH=(DV-GAU)V^{\top}$

for any $G\in\mathbb{R}^{l\times l}$. Since $V^{T}$ has linearly independent
rows, it follows that,

$\displaystyle D-GAH=0\quad\text{ if and only if }\quad DV-GAU=0.$ (6)

We extend $DV,AU\in\mathbb{R}^{l\times k}$ (where $k\leq l$) to invertible
square matrices $\widehat{DV},\widehat{AU}$, appending columns. Then

$\displaystyle G:=\widehat{DV}(\widehat{AU})^{-1}\in\mathbb{R}^{l\times l}$
(7)

is invertible with $\widehat{DV}=G\widehat{AU}$, and thus $DV=GAU$. The latter
implies $D=GAH$ by (6). Note that, again, the definition of $H$ via (5) is not
unique, but depends on the Gram-Schmidt procedure to derive $U,V$. Moreover,
if $k<l$ the definition of $G$ via (7) is not unique; even for fixed $U,V$.

##  4 Experimental results

Experiments were conducted to compare CS recovery performance using the
proposed compressive sensing factorization approach with benchmarks derived
using Gaussian or Bernoulli random matrices as sensing matrices for sparse
vectors. The entries of $A\in\mathbb{R}^{l\times n}$ in our experiments are
i.i.d. Bernoulli and Gaussian random numbers where the Gaussian random
variables are $\mathcal{N}(0,n^{-1})$ and the Bernoulli random variables are
$\pm\frac{1}{\sqrt{n}}$ with equal probability. In all experiments, the sparse
vectors were derived using the compressive sampling matched pursuit (CoSaMP)
algorithm [42]. The formulation of CoSaMP for our approach, with sparsifying
dictionary $D=GAH$, is

$\displaystyle\begin{cases}\displaystyle\min_{x_{i}}\|z_{i}-\mathcal{E}_{i}G^{-1}Dx_{i}\|^{2}_{2}\\\
\|x_{i}\|_{0}\leq k,\end{cases}$ (8)

while for the benchmark methods used for comparison it is

$\displaystyle\begin{cases}\displaystyle\min_{x_{i}}\|z_{i}-\mathcal{E}_{i}Ax_{i}\|^{2}_{2}\\\
\|x_{i}\|_{0}\leq k.\end{cases}$ (9)

Note that $z_{i}$ in (8) and (9) are obtained from $\mathcal{E}_{i}Ax_{i}$ and
$\mathcal{E}_{i}$ in (8) and (9) are the same. Following the fact that if
$D=GAH$ is plugged into (8), the only difference is $AH$ in (8) and $A$ in (9)
and $AH$ and $A$ have the same distribution, we can derive that applying
sensing matrix $\mathcal{E}_{i}G^{-1}$ to dictionary $D$ for sparse recovery
(8) yields a performance similar to (in terms of probability) the sparse
recovery using (9). To verify this, we performed experiments on three
sparsifying dictionaries in $\mathbb{R}^{128\times 1024}$: $D_{wavelet}$,
$D_{KSVD}$, and $D_{PKSVD}$ (a variant of K-SVD optimized to be a Parseval
tight frame). The latter two dictionaries were respectively derived via the
K-SVD and PK-SVD algorithms, using the same set of training vectors obtained
from the $256\times 256$ gray-scale test image “Boat”555Parameters of the
K-SVD process were set at $K=1024$, $L=64$, $numIteration=50$,
$InitializationMethod=^{\prime}GivenMatrix^{\prime}$, and $errorGoal=10^{-8}$.
Parameters for PK-SVD were set at $\Phi_{0}=\Psi_{0}=D_{KSVD}$, $L=64$,
$maxIter=50$, $\rho=[0.1,10^{8},10^{8}]$, and $t=10^{-10}$. . The image was
divided into overlapping patches of $16\times 16$ pixels with stride $2$ (in
each dimension) with overlaps of $4$ pixels (in each dimension), which
resulted in $3,721$ training patches. Since the stride is $2$ in each
dimension, each patch forms a vector of size $128$. The mean of each patch was
normalized to zero, whereupon the resulting patches were transformed into
vectors (via the vec operation) to form the set of training vectors. The
dictionary $D_{wavelet}$ is derived from the CDF 9-7 wavelet. The first $640$
(i.e., each level has $128$ columns and there are $5$ levels) columns in
$D_{wavelet}$ were obtained from the first $5$-level of the wavelet, whereas
the remainder were generated randomly. The column norms of the sparsifying
dictionaries were normalized. Figures 1 and 2 present the sparsifying
dictionaries and corresponding matrices $G$ for various $A$.

Figures 3, 4, and 5 illustrate CS recovery performance via plots indicating
the probability of successfully recovering a sparse vector versus the CS
ratio. The figures compare the performance of our approach (8) against the
benchmark (9) using Gaussian and Bernoulli sensing matrices at various levels
of sparsity. The horizontal and vertical axes, respectively, indicate the CS
ratio (i.e., $\frac{m}{n}\times 100\%$, where $m$ is the number of rows of
$\mathcal{E}_{i}$ and $n=1024$) and the probability of successfully recovering
a sparse vector. We claim that the true sparse vector $x$ can be recovered as
long as the estimate $\hat{x}$ satisfies $\|\hat{x}-x\|_{1}<1024\times
10^{-2}$.

##  5 Conclusions

This paper outlines a novel approach to CS involving the construction of
sensing matrix for a sparsifying dictionary, such that the composition of
sensing matrix and dictionary has the RIP. The proposed approach is based on
factorizing a sparsifying dictionary $D$ with the help of a matrix $A$ of
equal rank, and satisfying the concentration inequality to achieve $D=GAH$.
The factorization solutions for $G$ and $H$ are not unique, and $G$ is related
to the sensing matrix of $D$. The non-uniqueness raises the question which
factorization solution would provide the most stable result, as indicated by
the ability to recover the sparse vector by solving the inverse problem of
$AH$ via convex optimization. As shown in Figure 2, a lack of structure in $G$
raises concerns pertaining to the applicability of this approach to large
scale systems with restricted computation and storage resources.

Acknowledgements:  The original proof of Lemma 3 and Theorem 5 were shortened
by Assist Prof Andreas Heinecke, National University of Singapore. We truly
appreciate his comments and suggestions on the presentation of this paper.

## References

  * [1] W. B. Johnson and J. Lindenstrauss, “Extensions of lipschitz mappings into a hilbert space 26,” Contemporary mathematics, vol. 26, 1984.
  * [2] S. Dasgupta and A. Gupta, “An elementary proof of the johnson-lindenstrauss lemma,” International Computer Science Institute, Technical Report, vol. 22, no. 1, pp. 1–5, 1999.
  * [3] J. M. Duarte-Carvajalino and G. Sapiro, “Learning to sense sparse signals: Simultaneous sensing matrix and sparsifying dictionary optimization,” IEEE Transactions on Image Processing, vol. 18, no. 7, pp. 1395–1408, 2009.
  * [4] E. J. Candes and T. Tao, “Decoding by linear programming,” IEEE transactions on information theory, vol. 51, no. 12, pp. 4203–4215, 2005.
  * [5] E. J. Candes, “The restricted isometry property and its implications for compressed sensing,” Comptes rendus mathematique, vol. 346, no. 9-10, pp. 589–592, 2008.
  * [6] T. T. Cai, L. Wang, and G. Xu, “New bounds for restricted isometry constants,” IEEE Transactions on Information Theory, vol. 56, no. 9, pp. 4388–4394, 2010.
  * [7] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on information theory, vol. 52, no. 2, pp. 489–509, 2006\.
  * [8] D. L. Donoho, “Compressed sensing,” IEEE Transactions on information theory, vol. 52, no. 4, pp. 1289–1306, 2006.
  * [9] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, “The johnson-lindenstrauss lemma meets compressed sensing,” preprint, vol. 100, no. 1, pp. 1–9, 2006\.
  * [10] E. J. Candes, Y. C. Eldar, D. Needell, and P. Randall, “Compressed sensing with coherent and redundant dictionaries,” Applied and Computational Harmonic Analysis, vol. 31, no. 1, pp. 59–73, 2011.
  * [11] A. Cohen, W. Dahmen, and R. DeVore, “Compressed sensing and best k-term approximation,” Journal of the American mathematical society, vol. 22, no. 1, pp. 211–231, 2009.
  * [12] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Transactions on information theory, vol. 53, no. 12, pp. 4655–4666, 2007.
  * [13] H. Rauhut, “Compressive sensing and structured random matrices,” Theoretical foundations and numerical methods for sparse recovery, vol. 9, pp. 1–92, 2010.
  * [14] D. Achlioptas, “Database-friendly random projections,” in Proceedings of the twentieth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database systems, pp. 274–281, 2001.
  * [15] F. Krahmer and R. Ward, “New and improved johnson–lindenstrauss embeddings via the restricted isometry property,” SIAM Journal on Mathematical Analysis, vol. 43, no. 3, pp. 1269–1281, 2011.
  * [16] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, “A simple proof of the restricted isometry property for random matrices,” Constructive Approximation, vol. 28, no. 3, pp. 253–263, 2008.
  * [17] M. Aharon, M. Elad, and A. Bruckstein, “K-svd: An algorithm for designing overcomplete dictionaries for sparse representation,” IEEE Transactions on signal processing, vol. 54, no. 11, pp. 4311–4322, 2006.
  * [18] W.-L. Hwang, P.-T. Huang, B.-C. Kung, J. Ho, and T.-L. Jong, “Frame-based sparse analysis and synthesis signal representations and parseval k-svd,” IEEE Transactions on Signal Processing, vol. 67, no. 12, pp. 3330–3343, 2019.
  * [19] A. Cohen, I. Daubechies, and J.-C. Feauveau, “Biorthogonal bases of compactly supported wavelets,” Communications on pure and applied mathematics, vol. 45, no. 5, pp. 485–560, 1992.
  * [20] P. Frankl and H. Maehara, “The johnson-lindenstrauss lemma and the sphericity of some graphs,” Journal of Combinatorial Theory, Series B, vol. 44, no. 3, pp. 355–362, 1988.
  * [21] P. Indyk and R. Motwani, “Approximate nearest neighbors: towards removing the curse of dimensionality,” in Proceedings of the thirtieth annual ACM symposium on Theory of computing, pp. 604–613, 1998.
  * [22] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM review, vol. 43, no. 1, pp. 129–159, 2001.
  * [23] S. Arora, A. Bhaskara, R. Ge, and T. Ma, “More algorithms for provable dictionary learning,” arXiv preprint arXiv:1401.0579, 2014.
  * [24] B. Barak, J. A. Kelner, and D. Steurer, “Dictionary learning and tensor decomposition via the sum-of-squares method,” in Proceedings of the forty-seventh annual ACM symposium on Theory of computing, pp. 143–151, 2015\.
  * [25] Y. Zhai, Z. Yang, Z. Liao, J. Wright, and Y. Ma, “Complete dictionary learning via l4-norm maximization over the orthogonal group,” Journal of Machine Learning Research, vol. 21, no. 165, pp. 1–68, 2020.
  * [26] K. Schnass, “Dictionary learning-from local towards global and adaptive,” arXiv preprint arXiv:1804.07101, 2018.
  * [27] M. A. Herman and D. Needell, “Mixed operators in compressed sensing,” in 2010 44th Annual Conference on Information Sciences and Systems (CISS), pp. 1–6, IEEE, 2010.
  * [28] J. Ho and W.-L. Hwang, “Ambiguity and sensitivity in imprecise dictionaries for compressed sensing,” Signal Processing, vol. 178.
  * [29] D. L. Donoho and M. Elad, “Optimally sparse representation in general (nonorthogonal) dictionaries via l1 minimization,” Proceedings of the National Academy of Sciences, vol. 100, no. 5, pp. 2197–2202, 2003.
  * [30] M. Elad, “Optimized projections for compressed sensing,” IEEE Transactions on Signal Processing, vol. 55, no. 12, pp. 5695–5702, 2007.
  * [31] R. A. DeVore, “Deterministic constructions of compressed sensing matrices,” Journal of complexity, vol. 23, no. 4-6, pp. 918–925, 2007.
  * [32] W. U. Bajwa, J. D. Haupt, G. M. Raz, S. J. Wright, and R. D. Nowak, “Toeplitz-structured compressed sensing matrices,” in 2007 IEEE/SP 14th Workshop on Statistical Signal Processing, pp. 294–298, IEEE, 2007.
  * [33] J. Romberg, “Compressive sensing by random convolution,” SIAM Journal on Imaging Sciences, vol. 2, no. 4, pp. 1098–1128, 2009.
  * [34] E. J. Candes and Y. Plan, “A probabilistic and ripless theory of compressed sensing,” arXiv preprint arXiv:1011.3854, 2010.
  * [35] T. T. Do, L. Gan, N. H. Nguyen, and T. D. Tran, “Fast and efficient compressive sensing using structurally random matrices,” IEEE Transactions on signal processing, vol. 60, no. 1, pp. 139–154, 2011.
  * [36] A. Mousavi and R. G. Baraniuk, “Learning to invert: Signal recovery via deep convolutional networks,” in 2017 IEEE international conference on acoustics, speech and signal processing (ICASSP), pp. 2272–2276, IEEE, 2017\.
  * [37] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” IEEE Signal Processing Magazine, vol. 38, no. 2, pp. 18–44, 2021.
  * [38] J. Zhang and B. Ghanem, “Ista-net: Interpretable optimization-inspired deep network for image compressive sensing,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1828–1837, 2018.
  * [39] A. Bora, A. Jalal, E. Price, and A. G. Dimakis, “Compressed sensing using generative models,” in International Conference on Machine Learning, pp. 537–546, PMLR, 2017.
  * [40] W.-L. Hwang and S.-S. Tung, “Learning dnn networks using un-rectifying relu with compressed sensing application,” arXiv preprint arXiv:2101.06940, 2021\.
  * [41] C. Bao, J.-F. Cai, and H. Ji, “Fast sparsity-based orthogonal dictionary learning for image restoration,” in Proceedings of the IEEE International Conference on Computer Vision, pp. 3384–3391, 2013.
  * [42] D. Needell and J. A. Tropp, “Cosamp: Iterative signal recovery from incomplete and inaccurate samples,” Applied and computational harmonic analysis, vol. 26, no. 3, pp. 301–321, 2009.

(a) $D_{KSVD}\in\mathbb{R}^{128\times 1024}$

(b) $D_{PKSVD}\in\mathbb{R}^{128\times 1024}$

(c) $D_{wavelet}\in\mathbb{R}^{128\times 1024}$

Figure 1:  Visualization of sparsifying dictionaries with integer values
ranging from $0$ to $255$.

(a1) $G_{KSVD}$ (Gaussian) (a2) $G_{KSVD}$ (Bernoulli)

(a1) $G_{PKSVD}$ (Gaussian) (a2) $G_{PKSVD}$ (Bernoulli)

(a1) $G_{wavelet}$ (Gaussian) (a2) $G_{wavelet}$ (Bernoulli)

Figure 2:  Matrix $G\in\mathbb{R}^{128\times 128}$ of several sparsifying
dictionaries $D$ for various $A$. Recall that $D=GAH$ and the sensing matrix
of $D$ is $\mathcal{E}G^{-1}$. In (a1), (b1), and (c1), $A$ is a Gaussian
random matrix. In (a2), (b2), and (c2), $A$ is a Bernoulli random matrix. In
(a1) and (a2), $D=D_{KSVD}$. In (b1) and (b2), $D=D_{PKSVD}$. In (c1) and
(c2), $D=D_{wavelet}$. Note that the structure in the sub-figures is barely
discernable.

(a1) k = 10 (a2) k = 10

(b1) k = 12 (b2) k = 12

(c1) k = 14 (c2) k = 14

Figure 3:  Comparisons of CS recovery performance (i.e., the probability of
sparse vector recovery versus CS ratio) using sparsifying dictionary
$D_{KSVD}$. Red and blue curves were respectively obtained using the benchmark
(9) and our approach (8). Sparse vectors $x_{i}$ were randomly generated and
each point on the curve is the average of $2,000$ probability measurements.
The positions of non-zero coefficients of $x_{i}$ are uniformly distributed
and the values of the non-zero coefficients of $x_{i}$ are uniformly
distributed in $[-1,1]$. In (a1), (b1), and (c1), $A$ is a Gaussian random
matrix. In (a2), (b2), and (c2), $A$ is a Bernoulli random matrix.

(a1) k = 10 (a2) k = 10

(b1) k = 12 (b2) k = 12

(c1) k= 14 (c2) k= 14

Figure 4:  Comparisons of CS recovery performance (i.e., the probability of
sparse vector recovery versus CS ratio) using sparsifying dictionary
$D_{PKSVD}$. Red and blue curves were respectively obtained using the
benchmark (9) and our approach (8). Sparse vectors $x_{i}$ were randomly
generated and each point on the curve is the average of $2,000$ probability
measurements. The positions of non-zero coefficients of $x_{i}$ are uniformly
distributed and the values of the non-zero coefficients of $x_{i}$ are
uniformly distributed in $[-1,1]$. In (a1), (b1), and (c1), $A$ is a Gaussian
random matrix. In (a2), (b2), and (c2), $A$ is a Bernoulli random matrix.

(a1) k = 10 (a2) k = 10

(b1) k = 12 (b2) k = 12

(c1) k = 14 (c2) k = 14

Figure 5:  Comparisons of CS recovery performance (i.e., the probability of
sparse vector recovery versus CS ratio) using sparsifying dictionary
$D_{wavelet}$. Red and blue curves were respectively obtained using the
benchmark (9) and our approach (8). Sparse vectors $x_{i}$ were randomly
generated and each point on the curve is the average of $2,000$ probability
measurements. The positions of non-zero coefficients of $x_{i}$ are uniformly
distributed and the values of the non-zero coefficients of $x_{i}$ are
uniformly distributed in $[-1,1]$. In (a1), (b1), and (c1), $A$ is a Gaussian
random matrix. In (a2), (b2), and (c2), $A$ is a Bernoulli random matrix.


depending on the magnitude of the imposed transverse magnetic field, the
evolution of the interface under the Rayleigh–Taylor instability varies. In
the case of $\hat{H}_{\mathrm{0}}<\hat{H}_{\mathrm{c}}$, since the magnitude
of the magnetic field is small compared to the critical magnetic field, the
growth of the Rayleigh–Taylor instability is expected to be similar to the
non-magnetic case, but with a slightly reduced growth rate over time. This
behaviour is captured by the implemented numerical solver, as represented in
Fig. 16(a). On the other hand, when the applied magnetic field is of the order
of or greater than the calculated critical magnetic field, the Rayleigh–Taylor
instability is anticipated to reach a stable regime. Consequently, the
interface between the two fluids oscillates, and the penetration of the
heavier fluid to the lighter one is prevented. The solver successfully
reproduced these expected results, as demonstrated in Fig.16(a). Additionally,
it is apparent from Fig.16(a) that for higher values of the applied magnetic
field, the interface oscillates at a faster rate, consistent with the
dispersion relation described by Eq. (59).

The magnetic permeability ratio between the two fluids,
$\hat{\mu}_{\mathrm{m}}$, also influences the evolution of the Rayleigh–Taylor
instability. By keeping the imposed transverse magnetic field constant and by
varying $\hat{\mu}_{\mathrm{m}}$, the critical magnetic field,
$\hat{H}_{\mathrm{c}}$, changes. Consequently, if the Rayleigh–Taylor
instability becomes stable for a given $\hat{\mu}_{\mathrm{m}}$ and
$\hat{H}_{0}$, altering the value of $\hat{\mu}_{\mathrm{m}}$ can render it
unstable. It is for this reason that the magnetic permeability ratio has a
dual role of destabilizing and stabilizing the system, as was mentioned by
Awatshi [75]. To assess the solver’s capability to accurately handle different
magnetic permeability ratios across the interface, three cases of
$\hat{\mu}_{\mathrm{m}}=0.01$, $\hat{\mu}_{\mathrm{m}}=1$, and
$\hat{\mu}_{\mathrm{m}}=100$ are investigated for a constant imposed
$\hat{H}_{0}$, shown in Fig. 16(b). It can be visually confirmed that the case
of $\hat{\mu}_{\mathrm{m}}=1$ is the most unstable one, a finding that is
consistent with Eq. (59). The Rayleigh–Taylor instability transitions to a
stable regime when $\hat{\mu}_{\mathrm{m}}=0.01$. However, by changing the
magnetic permeability ratio to $100$, the system becomes unstable,
demonstrating the destabilizing effect of the $\hat{\mu}_{\mathrm{m}}$.

(a)

(b)

Figure 16:  (a) The interface evolution of the Rayleigh–Taylor instability
with $\hat{\rho}=1/3$ and $\hat{\mu}_{\mathrm{m}}=0.01$, at
$t\sqrt{\mathrm{At}}=0.25,0.5,0.75,1.25,1.5,$ and 1.75, for three conditions
of $\hat{H}_{\mathrm{0}}<\hat{H}_{\mathrm{c}}$,
$\hat{H}_{\mathrm{0}}\approx\hat{H}_{\mathrm{c}}$, and
$\hat{H}_{\mathrm{0}}>\hat{H}_{\mathrm{c}}$, along with the nonmagnetic case.
(b) The interface evolution for three conditions of
$\hat{\mu}_{\mathrm{m}}=0.01$, $\hat{\mu}_{\mathrm{m}}=1$, and
$\hat{\mu}_{\mathrm{m}}=100$, with a constant $H_{0}=5$ and
$\mu_{\mathrm{m,1}}=0.001$, at $t\sqrt{\mathrm{At}}=0.5,0.75,1,1.25,1.75,$ and
2.25. The simulations are performed on $32\times 128$ mesh resolution and time
step $\Delta t=5\times 10^{-4}/\sqrt{\mathrm{At}}$.

##  7 Conclusion

This paper presents a numerical development of a two-phase incompressible
solver for magnetic flows. The proposed numerical toolkit couples the
Navier–Stokes equations of hydrodynamics with Maxwell’s equations of
electromagnetism, enabling us to model the behaviour of magnetic flows in the
presence of a static magnetic field. In order to achieve this goal, first, a
detailed implementation of a second-order two-phase solver for incompressible
nonmagnetic flows is introduced. This two-phase solver utilizes a fifth-order
conservative level set method within the finite-difference framework to
capture the interface evolution during the simulation. The accuracy,
robustness, and mass conservation properties of the level set solver are
verified by conducting three test cases, namely, a rotating circle, a circle
in a deformation field, and Zalesak’s disk. Additionally, the incompressible
Navier–Stokes equations are modelled based on the projection method and solved
in the conservative form. The order of accuracy for the implemented two-phase
nonmagnetic solver is verified through the capillary wave test case. The
Rayleigh–Taylor instability test is also performed to evaluate the performance
of the solver for more complicated interface evolution and higher density
ratios. For the magnetic case, the reduced Maxwell equations under the
magnetostatic assumption are solved based on the vector potential formulation.
Subsequently, the primary two-phase solver is extended to account for magnetic
flows by incorporating the Lorentz force in the momentum equation. The
developed solver demonstrated the capability to handle high magnetic
permeability jumps across the interface. Three benchmarks were conducted to
evaluate the performance and robustness of the implemented two-phase solver
for magnetic flows in the presence of a static magnetic field. The first test
case explored droplet deformation under the influence of an external magnetic
field in a quiescent flow, examining the interplay between the Lorentz force
and capillary force at the interface across varying magnetic field amplitudes.
These findings were validated by comparing them with analytical solutions and
existing numerical and experimental results from the literature. In the
sheared ferrofluid droplet test, the impact of magnetic susceptibility on
droplet dynamics was studied. The results revealed that an increase in the
susceptibility of the ferrofluid droplet could influence its deformation and
rotation in low capillary regimes. However, in higher capillary flows, it was
observed that, under constant capillary and magnetic Bond number values,
increasing magnetic susceptibility could lead to droplet breakup. Furthermore,
the solver was employed to simulate Rayleigh–Taylor instability in magnetic
fluids. Numerical results were compared with linear stability analysis across
different magnetic field intensities and magnetic permeability ratios at the
interface. This analysis facilitated a discussion regarding the effects of
these variables in either stabilizing or destabilizing the Rayleigh–Taylor
instability.

## Acknowledgments

The authors would like to thank Victor Boniou for his valuable comments on the
surface tension implementation and curvature calculations. Paria Makaremi-
Esfarjani further acknowledges the helpful suggestions of Khashayar F. Kohan
on the initial draft of the manuscript. The authors thank Mathias Larrouturou
for useful feedback on the manuscript. This study was supported by the Natural
Science and Engineering Research Council of Canada (NSERC) through a NSERC
Discovery Grant and Mitacs through the Mitacs Accelerate program.

##  Appendix A Fifth-order weighted essentially non-oscillatory (WENO)
interpolation

The computation of the level set convective flux requires the interpolation of
cell center $\psi$ values to cell faces using an upwind-based scheme. In this
study, a fifth-order WENO interpolation is implemented in our numerical
solver, explained in greater detail in this section. The concept of the WENO
interpolation is similar to the WENO reconstruction. This method involves
using a weighted combination of three sub-stencils to calculate the flux
values, resulting in a high-order scheme. The weights are chosen to minimize
the contribution of stencils with discontinuities, thereby avoiding numerical
instability and providing non-oscillatory interpolation.

Figure 17:  The stencil for the fifth-order WENO interpolation scheme to find
$\tilde{\psi}_{\mathrm{L},\,x_{i+1/2}}$ and
$\tilde{\psi}_{\mathrm{R},\,x_{i-1/2}}$ values at $x_{i+1/2}$ and $x_{i-1/2}$
cell faces, respectively.

In Fig. 17, the one-dimensional stencil for the fifth-order interpolation in
$x_{i}-$direction is shown. Subscripts $\mathrm{L}$ (left) and $\mathrm{R}$
(right) are employed to denote the interpolated values at each cell face based
on upwinding considerations. As a result, cell values at
$\tilde{\psi}_{\mathrm{L},x_{i+1/2}}$ and
$\tilde{\psi}_{\mathrm{R},x_{i-1/2}}$ have the same five-points stencil, shown
in Fig. 17, and can be divided into three sub-stencils

$\mathrm{S}_{1}=\left\\{x_{i-2},\,x_{i-1},\,x_{i}\right\\},$ (63a)
$\mathrm{S}_{2}=\left\\{x_{i-1},\,x_{i},\,x_{i+1}\right\\},$ (63b) and
$\mathrm{S}_{3}=\left\\{x_{i},\,x_{i+1},\,x_{i+2}\right\\}.$ (63c)

Using the third-order interpolation, the cell face values can be obtained for
each sub-stencil as

$\psi^{(1)}_{\mathrm{L},\,x_{i+1/2}}=\frac{3}{8}\psi_{x_{i-2}}-\frac{5}{4}\psi_{x_{i-1}}+\frac{15}{8}\psi_{x_{i}},$
(64a)
$\psi^{(2)}_{\mathrm{L},\,x_{i+1/2}}=-\frac{1}{8}\psi_{x_{i-1}}+\frac{3}{4}\psi_{x_{i}}+\frac{3}{8}\psi_{x_{i+1}},$
(64b)
$\psi^{(3)}_{\mathrm{L},\,x_{i+1/2}}=\frac{3}{8}\psi_{x_{i}}+\frac{3}{4}\psi_{x_{i+1}}-\frac{1}{8}\psi_{x_{i+2}},$
(64c)

and

$\psi^{(1)}_{\mathrm{R},\,x_{i-1/2}}=-\frac{1}{8}\psi_{x_{i-2}}+\frac{3}{4}\psi_{x_{i-1}}+\frac{3}{8}\psi_{x_{i}},$
(65a)
$\psi^{(2)}_{\mathrm{R},\,x_{i-1/2}}=\frac{3}{8}\psi_{x_{i-1}}+\frac{3}{4}\psi_{x_{i}}-\frac{1}{8}\psi_{x_{i+1}},$
(65b)
$\psi^{(3)}_{\mathrm{R},\,x_{i-1/2}}=\frac{15}{8}\psi_{x_{i}}-\frac{5}{4}\psi_{x_{i+1}}+\frac{3}{8}\psi_{x_{i+2}}.$
(65c)

The fifth-order interpolation is then obtained by combining the calculated
third-order interpolation for three sub-stencils using non-linear weights,
defined as

$\tilde{\psi}_{\mathrm{L},\,x_{i+1/2}}=\sum_{k=1}^{3}\omega_{\mathrm{L},\,k}\,\psi^{(k)}_{\mathrm{L},\,x_{i+1/2}},\quad\textrm{and}\quad\tilde{\psi}_{\mathrm{R},\,x_{i-1/2}}=\sum_{k=1}^{3}\omega_{\mathrm{R},\,k}\,\psi^{(k)}_{\mathrm{R},\,x_{i-1/2}},$
(66)

where $\omega$ variables are the nonlinear weights adopted from Jiang and Shu
[40], given as

$\omega_{\mathrm{L},\,k}=\frac{\beta_{\mathrm{L},\,k}}{\sum_{m=1}^{3}\beta_{\mathrm{L},m}},\quad\textrm{and}\quad\omega_{\mathrm{R},\,k}=\frac{\beta_{\mathrm{R},\,k}}{\sum_{m=1}^{3}\beta_{\mathrm{R},m}},$
(67)

where

$\beta_{\mathrm{L},\,k}=\frac{\bar{\omega}_{\mathrm{L},\,k}}{\left(\epsilon+\mathrm{IS}_{k}\right)^{2}},\quad\textrm{and}\quad\beta_{\mathrm{R},\,k}=\frac{\bar{\omega}_{\mathrm{R},\,k}}{\left(\epsilon+\mathrm{IS}_{k}\right)^{2}}.$
(68)

Here parameter $\epsilon$ is used to prevent the denominators to become zero.
The value of $\epsilon$ is typically chosen to be between $10^{-5}$ and
$10^{-7}$ [40]. In this study, we used a value of $\epsilon=10^{-6}$.
Variables $\mathrm{IS}_{k}$ are the smooth indicators for each sub-stencil,
defined as

$\displaystyle\mathrm{IS}_{1}=\frac{1}{3}\left[\psi_{x_{i-2}}\left(4\psi_{x_{i-2}}-19\psi_{x_{i-1}}+11\psi_{x_{i}}\right)+\psi_{x_{i-1}}\left(25\psi_{x_{i-1}}-31\psi_{x_{i}}\right)+10\psi^{2}_{x_{i}}\right],$
(69a)
$\displaystyle\mathrm{IS}_{2}=\frac{1}{3}\left[\psi_{x_{i-1}}\left(4\psi_{x_{i-1}}-13\psi_{x_{i}}+5\psi_{x_{i+1}}\right)+13\psi_{x_{i}}\left(25\psi_{x_{i}}-\psi_{x_{i+1}}\right)+4\psi^{2}_{x_{i+1}}\right],$
(69b) and
$\displaystyle\mathrm{IS}_{3}=\frac{1}{3}\left[\psi_{x_{i}}\left(10\psi_{x_{i}}-13\psi_{x_{i+1}}+11\psi_{x_{i+2}}\right)+\psi_{x_{i+1}}\left(25\psi_{x_{i+1}}-19\psi_{x_{i+2}}\right)+4\psi^{2}_{x_{i+2}}\right],$
(69c)

where $\bar{\omega}_{\mathrm{L},\,k}$ and $\bar{\omega}_{\mathrm{R},\,k}$ are
the linear optimal weights for the fifth-order interpolation, calculated using
the Lagrange interpolation, and are determined as

$\bar{\omega}_{\mathrm{L},\,k}=\left\\{\frac{1}{16},\frac{10}{16},\frac{5}{16}\right\\},\quad\textrm{and}\quad\bar{\omega}_{\mathrm{R},\,k}=\left\\{\frac{5}{16},\frac{10}{16},\frac{1}{16}\right\\}.$
(70)

The same procedure can be applied in each direction for multi-dimensional
problems as well.

##  Appendix B Third-order Runge–Kutta scheme

The third-order, total variation diminishing (TVD) Runge–Kutta scheme
introduced by Gottlieb and Shu [77] is an explicit temporal integration
scheme, which can attenuate spurious oscillations appearing in the solution.

Consider a system of differential equations given by

$\frac{\partial\mathbf{U}}{\partial t}=\mathbf{F}(\mathbf{U},t)\\\ ,$ (71)

where $\mathbf{U}$ is the vector of conservative variables and $\mathbf{F}$ is
a right-hand-side operator. For a known solution, $\mathbf{U}^{n}$, the third-
order Runge–Kutta scheme approximates the solution at the next time step,
$\mathbf{U}^{n+1}$, using an intermediate step

$\mathbf{U}^{1}=\mathbf{U}^{n}+\Delta t\mathbf{F}(\mathbf{U}^{n},t),$ (72)

followed by

$\mathbf{U}^{2}=\frac{3}{4}\mathbf{U}^{n}+\frac{1}{4}\mathbf{U}^{1}+\frac{1}{4}\Delta
t\mathbf{F}(\mathbf{U}^{1},t),$ (73)

which in turn is followed by the final step

$\mathbf{U}^{n+1}=\frac{1}{3}\mathbf{U}^{n}+\frac{2}{3}\mathbf{U}^{2}+\frac{2}{3}\Delta
t\mathbf{F}(\mathbf{U}^{2},t),\\\ $ (74)

where $\Delta t$ denotes the time step.

##  Appendix C  Evaluation of level set solver numerical accuracy: Rotating
circle test

The rotating circle test is examined to study the numerical order of accuracy
of the level set solver. In this test case, a circular interface of radius
$r=0.15$, with its center initially located at $(x_{0},y_{0})=(0.25,0.25)$ is
considered. The circle rotates in the two-dimensional computational domain,
$[x,y]\in[-0.8,0.8]\times[-0.8,0.8]$, under the constant velocity field
$\left(u,v\right)=(y,-x)$. According to the initial condition of the problem,
the circle should return to its initial place after one revolution at
$t=2\pi$. The initial velocity field and the evolution of the rotating circle
from $t=0$ to $t=2\pi$ is represented in Fig. 18(a). The initial thickness of
the level set profile is set to be $\epsilon=(\Delta x)^{0.8}/2$, with the
constant time step equal to $\Delta t=0.001$ and $\Delta\tau=0.0005$ for
transport and re-initialization equations, respectively. The re-initialization
process is executed every $20$ steps, with a maximum of two iterations per
step. The simulation is performed for five different grid resolutions, and
$L_{2}$ and $L_{\infty}$ errors are determined for each case. The values of
$L_{2}$ and $L_{\infty}$ errors are calculated as

$L_{2}\,\mathrm{error}=\sqrt{\frac{\Sigma_{i=1}^{N}\left(\psi_{\mathrm{anal,}i}-\psi_{\mathrm{num},i}\right)^{2}}{N}}\
\ \ \ \ \ \mathrm{and}\ \ \ \ \ \
L_{\infty}\,\mathrm{error}=\max\big{|}\psi_{\mathrm{anal,}i}-\psi_{\mathrm{num},i}\big{|},$
(75)

where $\psi_{\mathrm{anal,}i}$ and $\psi_{\mathrm{num},i}$ denote the
analytical and numerical results at grid point $i$, respectively, and $N$ is
the total number of the grid points. In Fig. 18(b), the error values are
plotted versus the number of nodes in the logarithmic scale, and, hence, the
slope of the plot depicts the estimated convergence rate of the implemented
numerical scheme. As can be observed from Fig. 18(b), the order of accuracy is
between $4$ and $4.5$, which is pretty close to the excepted fifth-order
accuracy of the implemented method, even for the defined sharp interface case
($\epsilon<\Delta x/2$). However, the re-initialization step is expected to
affect the global accuracy of the level set solver.

(a)  (b)

Figure 18:  (a) The velocity field of the rotating circle test case along with
the $\psi=0.5$ location at $t=0,\pi/4,\pi/2,3\pi/4,$ and $2\pi$. (b) Order of
accuracy analysis for the implemented level set solver using the rotating
circle test case. The error values are computed for five different grid
resolutions, $20\times 20$, $40\times 40$, $80\times 80$, $160\times 160$, and
$320\times 320$, with a constant time step for all simulations.

##  Appendix D  Incompressible solver test cases

In this section, two test cases are presented, the Taylor–Green vortex and the
lid-driven cavity, to examine the order of accuracy and robustness of the
implemented incompressible solver.

###  D.1 Decay of a Taylor–Green vortex

The Taylor–Green vortex is a well-known test case in the literature to examine
the ability of an incompressible Navier–Stokes numerical solver to simulate
transient problems. In this study, the Taylor–Green test is investigated to
evaluate the order of accuracy of the implemented incompressible solver. The
analytical solution of this benchmark is given as

$u=\sin\left(kx\right)\cos\left(ky\right)\exp\left(-2\mu k^{2}t/\rho\right),$
(76a) $v=-\cos\left(kx\right)\sin\left(ky\right)\exp\left(-2\mu
k^{2}t/\rho\right),$ (76b) and
$p=\frac{1}{4}\rho\left[\cos\left(2kx\right)+\cos\left(2ky\right)\right]\exp\left(-4\mu
k^{2}t/\rho\right),$ (76c)

where $k$ is the wavenumber, set to $k=1$. The magnitude of the initial
velocity field is represented in Fig. 19(a), which will be exponentially
damped due to the presence of the viscosity. The simulation is performed on
the computational domain $[x,y]\in[0,2\pi]\times[0,2\pi]$, with a periodic
boundary condition implemented at all boundaries. The values of density and
dynamic viscosity are set to $1$ and $0.01$, respectively. The simulation is
run for five different grid resolutions, $16\times 16,32\times 32,64\times
64,128\times 128$, and $256\times 256$, with the constant time step, $\Delta
t=0.001$. The maximum error of the velocity field is calculated at $t=0.2$ for
each case. In Fig. 19(b), the error values are plotted against the number of
nodes in the logarithmic scale, demonstrating the estimated convergence rate
of approximately three for the implemented incompressible solver. It should be
noted that the global accuracy of the velocity field is affected by the
approximation of the pressure gradient and the discretization of the diffusion
term, which are second-order accurate herein. In Table 1, second-order
convergence is observed for the pressure field as expected. Also, the maximum
calculated error for the velocity divergence is quite small,
$\mathcal{O}\left(10^{-13}\right)$, showing the proper divergence-free
property of the velocity field.

(a)  (b)

Figure 19:  (a) Initial velocity field of the Taylor–Green test case. (b) Order of accuracy analysis for the implemented incompressible solver using the decaying Taylor–Green vortex test case. The error values are computed for five different grid resolutions, $16\times 16$, $32\times 32$, $64\times 64$, $128\times 128$, and $256\times 256$, with the constant time step, $\Delta t=0.001$, for all the simulations. Table 1:  Calculated $L_{\infty}$ error of the pressure and velocity divergence for different grid resolutions of the Taylor–Green test case. Mesh resolution |  Error of $p$ | Rate |  Error of $\left(\nabla\cdot\mathbf{u}\right)$
---|---|---|---
$16\times 16$ | $5.28\times 10^{-2}$ | —— | $9.92\times 10^{-14}$
$32\times 32$ | $1.69\times 10^{-2}$ | 1.64 | $2.50\times 10^{-13}$
$64\times 64$ | $5.44\times 10^{-3}$ | 1.63 | $3.96\times 10^{-13}$
$128\times 128$ | $1.50\times 10^{-3}$ | 1.85 | $7.54\times 10^{-13}$
$256\times 256$ | $3.91\times 10^{-4}$ | 1.94 | $9.85\times 10^{-13}$

###  D.2 Lid-driven cavity

The lid-driven cavity flow test case is used to evaluate the performance of
the implemented solver for complicated flow fields and its capability to reach
a steady state solution. This test case includes the square cavity,
$[x,y]\in[0,1]\times[0,1]$, initially filled with an incompressible fluid,
with the top wall moving with the constant velocity, $u=1$, while the other
three walls are stationary. The simulation is performed for two cases of
Reynolds numbers of $100$ and $1000$, with a no-slip boundary condition
imposed at all boundaries. The two simulations are run for the grid resolution
of $100\times 100$ and $\Delta t=0.001$. Figure 20(a) depicts the streamlines
of the velocity fields for $\mathrm{Re}=100$ and $\mathrm{Re}=1000$ cases at
$t=10$ and $t=50$, respectively. This figure visually confirms that the
calculated velocity field follows the expected behavior, which is the
formation of a primary vortex towards the center of the cavity and the
generation of smaller corner vortices at the bottom corners. For a better
evaluation, the computed results for the $u$ and $v$ velocities at the
vertical and horizontal centerlines, respectively, are compared with the
corresponding results from the study by Huang et al. [78], shown in Fig. 20(b)
and (c). As can be observed in this figure, there is a good agreement between
the results of the two studies for both cases of $\mathrm{Re}=100$ and
$\mathrm{Re}=1000$.

(a)

(b)  (c)

Figure 20: (a) Velocity streamlines of the lid-driven cavity test case for
[left] $\mathrm{Re}=100$ and [right] $\mathrm{Re}=1000$ cases, at $t=10$ and
$t=50$, respectively. Calculated results of the lid-driven cavity test case
for (b) $u$ along the vertical line passing through the center, and (c) $v$
along the horizontal line passing through the center, for $\mathrm{Re}=100$
(black) and $\mathrm{Re}=1000$ (blue) cases, at $t=10$ and $t=50$,
respectively. Results obtained from this study and the study by Huang et al.
[78] are plotted using a solid line and a circle symbol, respectively.

##  Appendix E  Formation of spurious currents: static droplet test

In this test case, a two-dimensional droplet is considered in a stationary
velocity field, $\mathbf{u}=0$. By setting the velocity field to zero and
ignoring gravity, the equilibrium solution of the Navier–Stokes equation will
result in $\Delta p=\mathbf{F}_{\sigma}=\sigma\kappa$, which satisfies the
pressure jump condition, $[p]=\sigma\kappa$. This equilibrium condition is the
well-known Laplace’s relation between pressure and surface tension forces of a
droplet in an equilibrium condition [79]. Therefore, the droplet is expected
to remain at rest since the pressure force balances the capillary force. If
the numerical solver cannot accurately calculate the curvature and, hence,
fails to recover the equilibrium solution, quasi-stationary velocity patterns
known as parasitic or spurious currents will appear in the solution [48].
Thus, this test case has been widely used by researchers to investigate the
capability of their solver to properly balance pressure and surface tension
forces across the interface and its accuracy in computing the interface
curvature [48]. Additionally, the viscosity discontinuity between two flows
results in discontinuity of the velocity derivative across the interface,
which is more pronounced for larger viscosity ratios. This discontinuity can
cause numerical errors in the divergence calculation of the velocity field,
which leads to a source of error for pressure [78]. The inaccuracy in pressure
calculation will generate errors when balancing forces across the interface
and leads to the formation of unphysical spurious currents. Therefore, the
spurious velocity formation is also investigated for different viscosity
ratios across the interface in this benchmark.

In the initial condition, a circle of diameter ${D}=0.4$ is considered at the
center of the two-dimensional computational domain $[0,1]\times[0,1]$, where a
no-slip boundary condition is imposed at all boundaries. The thickness of the
interface is set to $\epsilon=\Delta x/2$. In order to study the magnitude of
the spurious flow that will be generated due to the inaccuracies while
computing the interface curvature, the proposed test case by Desjardins and
Pitsch [80] is followed. In their test case, the surface tension coefficient
is set to $\sigma=1$, with the same viscosity of $\mu=0.1$ for both fluids and
a constant density ratio of one, i.e., $\rho_{1}=\rho_{2}$. The Laplace
number, represented by $\mathrm{La}=\sigma\rho D/\mu^{2}$, is then altered by
changing the density of both fluids to examine the performance of the solver
for various ratios of surface tension and viscous forces. The intensity of the
unwanted spurious currents appearing in the solution can be quantified by
measuring the capillary number, $\mathrm{Ca}=|u_{\mathrm{max}}|\mu/\sigma$,
for different Laplace numbers. The mesh resolution is set to $32\times 32$,
and the calculated capillary numbers at a non-dimensional time
$t\sigma/\left(\mu D\right)=250$ for five different Laplace numbers are
reported in Table 2. For each Laplace number, the time step should be
determined according to the CFL constraint introduced in Eq. (40). According
to the calculated capillary numbers, the magnitude of spurious currents, even
after a long simulation time, is minimal and independent of the Laplace
number, confirming the well-balanced results between pressure and surface
tension forces across the interface in the solver. It is worth mentioning that
the order of calculated capillary numbers for different Laplace numbers is
close to the ones obtained by Desjardins and Pitsch [80], showing that the
magnitude of the spurious errors is consistent with other studies in the
literature.

Table 2:  Calculated Ca values for different La numbers for the static droplet test case on the mesh resolution of $32\times 32$, with density and viscosity ratios of one. $\rho_{1}(\rho_{2})$ | La | Ca
---|---|---
$0.3$ | $12$ | $5.029\times 10^{-5}$
$3$ | $120$ | $3.378\times 10^{-5}$
$30$ | $1200$ | $3.612\times 10^{-5}$
$300$ | $12000$ | $9.428\times 10^{-5}$
$3000$ | $120000$ | $2.657\times 10^{-4}$

In order to investigate the influence of the shear stress discontinuity across
the interface due to the viscosity jump on the formation of spurious velocity,
a benchmark akin to that of Huang et al. [78] is also adopted. In this test
case, a circle with the diameter of ${D}=2$ in a computational domain
$[0,8]\times[0,8]$ with the grid resolution of $101\times 101$ is considered.
For the initial condition of $\sigma=730$, $\rho_{1}=\rho_{2}=1$, and
$\mu_{1}=0.001$, the calculated $L_{2}$ and $L_{\infty}$ error norms of the
velocity field after $1$ and $1000$ time steps are given in Table 3 for four
different viscosity ratios of $1$, $10$, $100$, and $1000$, along with the
inviscid case. It is appreciated from Table 3 that even for large viscosity
ratios, the intensity of spurious currents is small and is close to the
inviscid case. In other words, the strength of the parasite current is
independent of the viscosity jump, and the solver accurately balances pressure
and surface tension forces across the interface in the presence of a shear
stress jump. Also, the calculated $L_{2}$ and $L_{\infty}$ errors for both
cases, after 1 and 1000 time steps, are in the same order as those reported by
Huang et al. study [78].

Table 3:  Calculated $L_{2}$ and $L_{\infty}$ norms of the velocity field for four different viscosity ratios along with the inviscid case for the static droplet test case. |  at $1^{\mathrm{st}}$ time step |  at $1000^{\mathrm{th}}$ time step
---|---|---
$\mu_{2}/\mu_{1}$ | $L_{2}$ | $L_{\infty}$ | $L_{2}$ | $L_{\infty}$
0 | $2.198\times 10^{-5}$ | $2.286\times 10^{-4}$ | $1.534\times 10^{-3}$ | $7.919\times 10^{-3}$
1 | $6.608\times 10^{-6}$ | $6.903\times 10^{-5}$ | $4.619\times 10^{-4}$ | $2.522\times 10^{-3}$
10 | $6.608\times 10^{-6}$ | $6.903\times 10^{-5}$ | $4.618\times 10^{-4}$ | $2.521\times 10^{-3}$
100 | $6.608\times 10^{-6}$ | $6.903\times 10^{-5}$ | $4.613\times 10^{-4}$ | $2.509\times 10^{-3}$
1000 | $6.608\times 10^{-6}$ | $6.903\times 10^{-5}$ | $4.562\times 10^{-4}$ | $2.394\times 10^{-3}$

## References

  * [1] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost fluid method for simulating turbulent atomization, Journal of Computational Physics 227 (2008) 8395–8416.
  * [2] R. Clift, J. R. Grace, M. E. Weber, Bubbles, Drops, and Particles, Courier Corporation, 2005.
  * [3] S. Radman, C. Fiorina, A. Pautz, Development of a novel two-phase flow solver for nuclear reactor analysis: algorithms, verification and implementation in OpenFOAM, Nuclear Engineering and Design 379 (2021) 111178.
  * [4] A. Prosperetti, G. Tryggvason, Computational Methods for Multiphase Flow, Cambridge University Press, 2009.
  * [5] M. Bussmann, D. B. Kothe, J. M. Sicilian, Modeling high density ratio incompressible interfacial flows, in: Fluids Engineering Division Summer Meeting, Vol. 36150, 2002, pp. 707–713.
  * [6] X. Li, Z.-Q. Dong, Y. Li, L.-P. Wang, X.-D. Niu, H. Yamaguchi, D.-C. Li, P. Yu, A fractional-step lattice Boltzmann method for multiphase flows with complex interfacial behavior and large density contrast, International Journal of Multiphase Flow 149 (2022) 103982.
  * [7] B. Lai, Z. Ke, Z. Wang, R. Zhu, R. Gao, Y. Mao, Y. Zhang, A lattice Boltzmann front-tracking interface capturing method based on neural network for gas-liquid two-phase flow, International Journal of Computational Fluid Dynamics 37 (1) (2023) 49–66.
  * [8] T. Fonty, M. Ferrand, A. Leroy, A. Joly, D. Violeau, Mixture model for two-phase flows with high density ratios: A conservative and realizable SPH formulation, International Journal of Multiphase Flow 111 (2019) 158–174.
  * [9] X. Cui, W. G. Habashi, V. Casseau, Multiphase SPH modelling of supercooled large droplets freezing on aircraft surfaces, International Journal of Computational Fluid Dynamics 35 (1-2) (2021) 79–92.
  * [10] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method), Journal of Computational Physics 152 (2) (1999) 457–492.
  * [11] J. U. Brackbill, D. B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of Computational Physics 100 (2) (1992) 335–354.
  * [12] S. Mirjalili, S. S. Jain, M. Dodd, Interface-capturing methods for two-phase flows: An overview and recent developments, Center for Turbulence Research Annual Research Briefs 2017 (117-135) (2017) 13.
  * [13] S. O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of Computational Physics 100 (1) (1992) 25–37.
  * [14] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annual Review of Fluid Mechanics 31 (1) (1999) 567–603.
  * [15] J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Vol. 3, Cambridge University Press, 1999.
  * [16] M. Sussman, E. G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows, Journal of Computational Physics 162 (2) (2000) 301–337.
  * [17] W. Meng, L. Liao, M. Chen, C.-h. Yu, J. Li, R. An, An enhanced CLSVOF method with an algebraic second-reconstruction step for simulating incompressible two-phase flows, International Journal of Multiphase Flow 154 (2022) 104151.
  * [18] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A hybrid particle level set method for improved interface capturing, Journal of Computational Physics 183 (1) (2002) 83–116.
  * [19] M. Herrmann, A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids, Journal of Computational Physics 227 (4) (2008) 2674–2706.
  * [20] F. Gibou, R. Fedkiw, S. Osher, A review of level-set methods and some recent applications, Journal of Computational Physics 353 (2018) 82–109.
  * [21] D. Chen, X. Tong, B. Xie, F. Xiao, Y. Li, An accurate and efficient multiphase solver based on THINC scheme and adaptive mesh refinement, International Journal of Multiphase Flow 162 (2023) 104409.
  * [22] E. Olsson, G. Kreiss, A conservative level set method for two phase flow, Journal of Computational Physics 210 (1) (2005) 225–246.
  * [23] E. Olsson, G. Kreiss, S. Zahedi, A conservative level set method for two phase flow II, Journal of Computational Physics 225 (1) (2007) 785–807.
  * [24] A. J. Chorin, A numerical method for solving incompressible viscous flow problems, Journal of Computational Physics 135 (2) (1997) 118–125.
  * [25] P. A. Davidson, An Introduction to Magnetohydrodynamics, Cambridge University Press, 2001.
  * [26] G.-S. Jiang, C.-C. Wu, A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics, Journal of Computational Physics 150 (2) (1999) 561–594.
  * [27] C.-C. Wu, A high order WENO finite difference scheme for incompressible fluids and magnetohydrodynamics, Geophysical and Astrophysical Fluid Dynamics 101 (1) (2007) 37–61.
  * [28] P. Makaremi-Esfarjani, A. Najafi-Yazdi, Characteristic boundary conditions for magnetohydrodynamic equations, Computers & Fluids 241 (2022) 105461.
  * [29] H. Huang, A. Ying, M. Abdou, 3D MHD free surface fluid flow simulation based on magnetic-field induction equations, Fusion Engineering and Design 63 (2002) 361–368.
  * [30] D. Gao, N. B. Morley, V. Dhir, Understanding magnetic field gradient effect from a liquid metal droplet movement, Journal of Fluids Engineering 126 (1) (2004) 120–124.
  * [31] T. Tagawa, Numerical simulation of two-phase flows in the presence of a magnetic field, Mathematics and Computers in Simulation 72 (2-6) (2006) 212–219.
  * [32] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, Journal of Computational Physics 114 (1) (1994) 146–159.
  * [33] J. Yang, S. Mao, X. He, X. Yang, Y. He, A diffuse interface model and semi-implicit energy stable finite element method for two-phase magnetohydrodynamic flows, Computer Methods in Applied Mechanics and Engineering 356 (2019) 435–464.
  * [34] S. Afkhami, Y. Renardy, M. Renardy, J. S. Riffle, T. St Pierre, Field-induced motion of ferrofluid droplets through immiscible viscous media, Journal of Fluid Mechanics 610 (2008) 363–380.
  * [35] S. Afkhami, A. Tyler, Y. Renardy, M. Renardy, T. S. Pierre, R. Woodward, J. S. Riffle, Deformation of a hydrophobic ferrofluid droplet suspended in a viscous medium under uniform magnetic fields, Journal of Fluid Mechanics 663 (2010) 358–384.
  * [36] M. Majidi, M. A. Bijarchi, A. G. Arani, M. H. Rahimian, M. B. Shafii, Magnetic field-induced control of a compound ferrofluid droplet deformation and breakup in shear flow using a hybrid lattice boltzmann-finite difference method, International Journal of Multiphase Flow 146 (2022) 103846.
  * [37] Y. Morinishi, T. S. Lund, O. V. Vasilyev, P. Moin, Fully conservative higher order finite difference schemes for incompressible flow, Journal of Computational Physics 143 (1) (1998) 90–124.
  * [38] O. Desjardins, G. Blanquart, G. Balarac, H. Pitsch, High order conservative finite difference scheme for variable density low mach number turbulent flows, Journal of Computational Physics 227 (2008) 7125–7159.
  * [39] X.-D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, Journal of Computational Physics 115 (1) (1994) 200–212.
  * [40] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, Journal of Computational Physics 126 (1) (1996) 202–228.
  * [41] Y. Su, S. H. Kim, An improved consistent, conservative, non-oscillatory and high order finite difference scheme for variable density low Mach number turbulent flow simulation, Journal of Computational Physics 372 (2018) 202–219.
  * [42] R. Nourgaliev, T. Theofanous, High-fidelity interface tracking in compressible flows: Unlimited anchored adaptive level set, Journal of Computational Physics 224 (2) (2007) 836–866.
  * [43] K. Mizuno, M. Asahara, T. Kamiya, T. Miyasaka, Finite difference and reinitialization methods with level set to interfacial area transport equations for gas–liquid two-phase flows, International Journal of Computational Fluid Dynamics 36 (5) (2022) 361–383.
  * [44] C. D. Pierce, Progress-variable approach for large-eddy simulation of turbulent combustion, Ph.D. thesis, Stanford University (2001).
  * [45] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, Journal of Computational Physics 59 (2) (1985) 308–323.
  * [46] S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith, Efficient management of parallelism in object oriented numerical software libraries, in: E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 1997, pp. 163–202.
  * [47] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows, Cambridge University Press, 2011.
  * [48] V. Boniou, T. Schmitt, A. Vié, Comparison of interface capturing methods for the simulation of two-phase flow in a unified low-Mach framework, International Journal of Multiphase Flow 149 (2022) 103957.
  * [49] M. Finster, H. Chaves, R. Schwarze, Numerical simulation of the accumulation process in modulated liquid jets of defined volume using the volume-of-fluid methodology, International Journal of Multiphase Flow 160 (2023) 104356.
  * [50] A. Prosperetti, Motion of two superposed viscous fluids, Physics of Fluids 24 (7) (1981) 1217–1223.
  * [51] Z. Huang, G. Lin, A. M. Ardekani, Consistent, essentially conservative and balanced-force phase-field method to model incompressible two-phase flows, Journal of Computational Physics 406 (2020) 109192.
  * [52] H. Ding, P. D. Spelt, C. Shu, Diffuse interface model for incompressible two-phase flows with large density ratios, Journal of Computational Physics 226 (2) (2007) 2078–2095.
  * [53] J.-L. Guermond, L. Quartapelle, A projection FEM for variable density incompressible flows, Journal of Computational Physics 165 (1) (2000) 167–188.
  * [54] R. E. Rosensweig, Ferrohydrodynamics, Courier Corporation, 2013.
  * [55] M. A. Bijarchi, A. Favakeh, S. Alborzi, M. B. Shafii, Experimental investigation of on-demand ferrofluid droplet generation in microfluidics using a pulse-width modulation magnetic field with proposed correlation, Sensors and Actuators B: Chemical 329 (2021) 129274.
  * [56] O. T. Mefford, R. C. Woodward, J. D. Goff, T. Vadala, T. G. S. Pierre, J. P. Dailey, J. S. Riffle, Field-induced motion of ferrofluids through immiscible viscous media: Testbed for restorative treatment of retinal detachment, Journal of Magnetism and Magnetic Materials 311 (1) (2007) 347–353.
  * [57] P. Voltairas, D. Fotiadis, C. Massalas, Elastic stability of silicone ferrofluid internal tamponade (SFIT) in retinal detachment surgery, Journal of Magnetism and Magnetic Materials 225 (1-2) (2001) 248–255.
  * [58] M. A. Bijarchi, M. B. Shafii, Experimental investigation on the dynamics of on-demand ferrofluid drop formation under a pulse-width-modulated nonuniform magnetic field, Langmuir 36 (26) (2020) 7724–7740.
  * [59] M. Zarei Saleh Abad, M. Ebrahimi-Dehshali, M. A. Bijarchi, M. B. Shafii, A. Moosavi, Visualization of pool boiling heat transfer of magnetic nanofluid, Heat Transfer—Asian Research 48 (7) (2019) 2700–2713.
  * [60] J.-C. Bacri, D. Salin, R. Massart, Study of the deformation of ferrofluid droplets in a magnetic field, Journal de Physique Lettres 43 (6) (1982) 179–184.
  * [61] X. Li, P. Yu, X.-D. Niu, D.-C. Li, H. Yamaguchi, A magnetic field coupling lattice boltzmann model and its application on the merging process of multiple-ferrofluid-droplet system, Applied Mathematics and Computation 393 (2021) 125769.
  * [62] J. Bacri, D. Salin, Instability of ferrofluid magnetic drops under magnetic field, Journal de Physique Lettres 43 (17) (1982) 649–654.
  * [63] L. H. Cunha, I. R. Siqueira, T. F. Oliveira, H. D. Ceniceros, Field-induced control of ferrofluid emulsion rheology and droplet break-up in shear flows, Physics of Fluids 30 (12) (2018).
  * [64] M. R. Hassan, J. Zhang, C. Wang, Deformation of a ferrofluid droplet in simple shear flows under uniform magnetic fields, Physics of Fluids 30 (9) (2018).
  * [65] M. R. Hassan, C. Wang, Magnetic field induced ferrofluid droplet breakup in a simple shear flow at a low reynolds number, Physics of Fluids 31 (12) (2019).
  * [66] M. R. Kennedy, C. Pozrikidis, R. Skalak, Motion and deformation of liquid drops, and the rheology of dilute emulsions in simple shear flow, Computers & fluids 23 (2) (1994) 251–278.
  * [67] S. Guido, M. Villone, Three-dimensional shape of a drop under simple shear flow, Journal of Rheology 42 (2) (1998) 395–415.
  * [68] G. I. Taylor, The viscosity of a fluid containing small drops of another fluid, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 138 (834) (1932) 41–48.
  * [69] G. I. Taylor, The formation of emulsions in definable fields of flow, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 146 (858) (1934) 501–523.
  * [70] W. C. Jesus, A. M. Roma, H. D. Ceniceros, Deformation of a sheared magnetic droplet in a viscous fluid, Communications in Computational Physics 24 (2) (2018) 332–355.
  * [71] Y. Hu, D. Li, X. Niu, Phase-field-based lattice Boltzmann model for multiphase ferrofluid flows, Physical Review E 98 (3) (2018) 033301.
  * [72] A. L. Velikovich, F. Cochran, J. Davis, Suppression of Rayleigh–Taylor instability in Z-pinch loads with tailored density profiles, Physical Review Letters 77 (5) (1996) 853.
  * [73] E. Harris, Rayleigh–Taylor instabilities of a collapsing cylindrical shell in a magnetic field, The Physics of Fluids 5 (9) (1962) 1057–1062.
  * [74] A. Velikovich, P. Schmit, Bell–Plesset effects in Rayleigh–Taylor instability of finite-thickness spherical and cylindrical shells, Physics of Plasmas 22 (12) (2015) 122711.
  * [75] M. K. Awasthi, Viscous potential flow analysis of magnetohydrodynamic Rayleigh–Taylor instability with heat and mass transfer, International Journal of Dynamics and Control 2 (2014) 254–261.
  * [76] Y. O. El-Dib, Nonlinear hydrodynamic Rayleigh–Taylor instability of viscous magnetic fluids: effect of a tangential magnetic field, Journal of Plasma Physics 51 (1) (1994) 1–11.
  * [77] S. Gottlieb, C.-W. Shu, Total variation diminishing Runge–Kutta schemes, Mathematics of Computation 67 (221) (1998) 73–85.
  * [78] Z. Huang, G. Lin, A. M. Ardekani, A mixed upwind/central WENO scheme for incompressible two-phase flows, Journal of Computational Physics 387 (2019) 455–480.
  * [79] S. Popinet, Numerical models of surface tension, Annual Review of Fluid Mechanics 50 (2018) 49–75.
  * [80] O. Desjardins, H. Pitsch, A spectrally refined interface approach for simulating multiphase flows, Journal of Computational Physics 228 (5) (2009) 1658–1677.