Scalable adaptive algorithms for next-generation multiphase simulations

Kumar Saurabh§

Iowa State University, USA

Hari Sundar

University of Utah, USA

Masado Ishii§

University of Utah, USA

Makrand A. Khanwale

Stanford University, USA

Baskar Ganapathysubramanian

Iowa State University, USA

Equal contribution

The accuracy of multiphysics simulations is strongly contingent up on the finest resolution of mesh used to resolve the interface. However, the increased resolution comes at a cost of inverting a larger matrix size. In this work, we propose algorithmic advances that aims to reduce the computational cost without compromising on the physics by selectively detecting the key regions of interest (droplets/filaments) that requires significantly higher resolution.

The overall framework uses an adaptive octree–based mesh generator, which is integrated with  's linear algebra solver. We demonstrate the scaling of the framework up to 114,688 processes on TACC  . Finally we deploy the framework to simulate primary jet atomization on an equivalent 35 trillion grid points – 64$\times$ greater than the state–of–the–art simulations.

multiphase simulations, octree, massively parallel algorithms.

§ INTRODUCTION

Understanding the mechanisms of phase interactions in two-phase flows is critical to answering fundamental and pressing questions in science and technology. In particular, insights into phase interactions may lead to accurate, low-cost coarse-scale models for large systems (e.g., chemical/biological reactors). Further, they can be used for the optimization-based design of micro-scale systems (e.g., bio-microfluidics and advanced manufacturing using multiphase flows). We need physically accurate models and simulations that capture the varying time and length scales of phase interactions to achieve this understanding.

Interface resolved two-phase modeling can be characterized into two main categories – sharp–interface methods and diffuse–interface methods. The sharp interface methods rely on representing the interface with a sharp, discontinuous function (e.g., Volume of Fluids (VOF)), whereas the diffuse interface method involves using a diffuse continuous representation of the interface (e.g., conservative diffuse interface, Cahn–Hillard Navier–Stokes (CHNS)).

The class of sharp interface methods has been state–of–the–art for high-fidelity simulations of turbulent multiphase flows  [Fuster et al., 2009, Tomar et al., 2010, Pairetti et al., 2020, Herrmann, 2010]. These methods involve solving a PDE for transporting the discontinuous sharp volume fraction function.  An interface reconstruction procedure is generally required to construct normals.
The sharp interface methods suffer from artificial breakup – also known as numerical surface tension, when the interfacial features (droplets/filaments) are comparable in length scale $r$ to the grid size $\Delta x$  [Gorokhovski and Herrmann, 2008]. There has been some progress in the development of computational techniques that reduces this effect involving interface reconstruction techniques [Lu and Tryggvason, 2018, Chiodi, 2020], Moment-of-Fluid methods [Jemison et al., 2015], Manifold death algorithms [Chirco et al., 2022]. However, these methods are limited to sharp interface methods and are still computationally complex.

Diffuse interface methods, particularly CHNS, provide a diffuse interface representation and avoid any interface reconstruction process. However, diffuse interface methods—under low mesh resolutions—also suffer from artificial breakup and mass loss, but for completely different reasons. The limitation arises because of the thickness of the diffuse interface $\epsilon$. In the case of CHNS, when the ratio of $\epsilon/r$ is $\mathcal{O}(1)$, it can lead to bound violations in addition to artificial breakup, and these structures might disappear and get absorbed into the larger structures (due to coarsening of Cahn-Hilliard)

So, both the sharp and the diffuse interface methods for insufficient resolution face the challenge of artificial breakup. The main contributions of this work are as follows:

* We propose a novel scalable algorithms to identify regions of interest, i.e., regions where $\epsilon/r  \approx \mathcal{O}(1)$ and selectively increase the mesh resolution in these identified regions.

Providing targeted resolution is an important feature for performing cost-effective simulations

* In the case of the multiphase simulations, the levels of the mesh can vary by several orders of magnitudes to resolve the interfacial and flow physics. We tailor existing octree refinement and coarsening algorithms to accelerate remeshing and decrease the associated overhead, especially for multi-level refinement, where the element sizes drop substantially. This contrasts existing approaches, where refinement or coarsening of the octrees is done level by level. [Saurabh et al., 2021, Saurabh et al., 2021, Fernando et al., 2018, Popinet, 2003, Popinet, 2009, Bangerth et al., 2012]

* Further, we analyze the data movement during matrix and vector assembly of FEM operators and propose a new data layout for local matrix and vector assembly that is more amenable to avoiding cache misses.

* We showcase the application of the proposed algorithm for the canonical problem of primary jet atomization, where we perform a simulation equivalent to $35$ trillion grid points on a uniform mesh – $64\times$ more resolved than the current state–of–the–art simulation.

§ METHODS

In this work, we present algorithmic developments that are applicable

for all interface capturing methods,
however, we specifically deploy a thermodynamically consistent Cahn–Hillard Navier–Stokes (CHNS) for the numerical simulations.

§.§ Numerical methods

Governing equations:

Let $\phi$ be the phase field variable that tracks the location of the phases and varies
smoothly between +1 to -1 with a characteristic diffuse interface thickness $\epsilon$, and $v_i$ be the $i^{th}$ component of the mixture velocity of the two phases. A thermodynamically consistent coupled CHNS equations [Khanwale et al., 2021] are written as follows:

* Momentum equation:

\begin{equation}
\begin{split}
&\pd{\left(\rho(\phi) v_i\right)}{t} + \pd{\left(\rho(\phi)v_iv_j\right)}{x_j} + \frac{1}{Pe}\pd{\left(J_jv_i\right)}{x_j}+\frac{Cn}{We} \pd{}{x_j}\left({\pd{\phi}{x_i}\pd{\phi}{x_j}}\right)
\\
&+\frac{1}{We}\pd{p}{x_i} - \frac{1}{Re}\pd{}{x_j}\left({\eta(\phi)\pd{v_i}{x_j}}\right) - \frac{\rho(\phi)\hat{{g_i}}}{Fr} = 0,\\
& \quad\quad \text{where,} \quad J_i = \frac{\left(\rho_- - \rho_+\right)}{2\;\rho_+ Cn} \, m(\phi),
\end{split}
\end{equation}

* Solenoidality and Continuity:

\begin{equation}
\pd{v_i}{x_i} = 0, \quad \pd{\rho(\phi)}{t} + \pd{\left(\rho(\phi)v_i\right)}{x_i}+
\frac{1}{Pe} \pd{J_i}{x_i} = 0,
\end{equation}

* Cahn–Hillard Equations

\begin{equation}
\begin{split}
&\pd{\phi}{t} + \pd{\left(v_i \phi\right)}{x_i} - \frac{1}{PeCn} \pd{}{x_i}\left(m(\phi){\pd{\mu}{x_i}}\right) \\
&\mu = \psi'(\phi) - Cn^2 \pd{}{x_i}\left({\pd{\phi}{x_i}}\right)
\end{split}
\end{equation}

where $\rho(\phi)$ is the non–dimensional mixture density given by $\left({\rho_+ - \;\rho_-}/{2\rho_+}\right) \phi + \left({\rho_+ + \;\rho_-}/{2\rho_+}\right)$, the non–dimensional mixture viscosity $\eta\left(\phi\right)$ is given by $\left({\eta_+ - \;\eta_-}/{2\eta_+}\right) \phi + \left({\eta_+ + \;\eta_-}/{2\eta_+}\right)$.

We use a degenerate mobility $m(\phi)$ given by $\sqrt{1-\phi^2}$.

Non-dimensional parameters are as follows: Peclet, $Pe = \frac{u_{r} L_{r}^2}{m_{r}\sigma}$; Reynolds, $Re = \frac{u_{r} L_{r}}{\nu_{r}}$; Weber, $We = \frac{\rho_{r}u_{r}^2 L_{r}}{\sigma}$; Cahn, $Cn = \frac{\varepsilon}{L_{r}}$; and Froude, $Fr = \frac{u_{r}^2}{gL_{r}}$, with $u_{r}$ and $L_r$ denoting the reference velocity and length, respectively.

Numerical discretization:

In this work, we utilize the two-block projection scheme proposed in  Khanwale et al., 2021. The proposed scheme is based on projection-based pressure splitting of Navier–Stokes (NS). Overall, the scheme involves four solutions per block:

* CH-Solve: As the first step, we solve the advective Cahn–Hillard equation. This is solved as a non–linear solver in a fully implicit manner.

* NS-Solve:

We compute mixture density and viscosity as a function of $\phi$
from the CH-Solve in the previous step and use it to solve the momentum equations.
We use a semi-implicit discretization based on the Crank-Nicolson scheme that results in linearized momentum equations.

* PP-Solve:
The projection-based splitting results in a variable density pressure Poisson equation that allows us to compute pressure at the current time level.

* VU-solve: As the last step for each block, we correct the velocity predicted in the NS-solve using the Helmholtz equations to get the solenoidal velocity.

To minimize the memory footprint, we split the final VU-solve to update velocity once per direction. This reduces the assembled matrix size from $N\times \mathrm{DIM} \times k$ to $N \times k$, where $N$ is the total number of independent DOFs, $k$ is the number of non–zeros per row, and $\mathrm{DIM}$ is the spatial dimension (2D/3D). We note that the mass matrix, as a result of VU-solve, does not need to be recomputed for each of the $\mathrm{DIM}$ separately and is reused till the mesh does not change. Once the matrix is assembled, no subsequent call to  and   is made, reducing the preconditioner computation cost.

The above second-order semi-implicit formulations have been demonstrated to have a milder time-stepping restriction $\Delta t$ compared to the fully explicit counterparts. The explicit parts of the discretization allow for the linearization of NS, which avoids an expensive setup of Newton iteration for NS.

Although the proposed numerical scheme can be integrated with arbitrary order basis functions with Continuous Galerkin (CG) elements, we have restricted the deployment of the numerical algorithm to the linear basis functions (spatially second-order convergence).

Local Cahn:

The diffuse interface methods, including CHNS,  under low mesh resolution, suffers from artificial breakup and mass loss. Particularly for CHNS, this limitation arises when the interface thickness $\epsilon$ has a comparable length scale to the length scale of fluid features (droplets/filaments/sheets), leading to the artificial breakup. Our simulations use the $Cn$ number as a non-dimensional $\epsilon$.

A smaller $Cn$ represents a thinner diffuse interface which allows the capture of more accurate physics. A smaller $Cn$ compared to the other interface features would prevent an artificial breakup.

However, a decrease in $Cn$ everywhere comes at an additional cost of increased mesh resolution. A rough estimate indicates that a $2\times$ decrease in $Cn$ leads to approximately $8\times$ increase in mesh resolution. In this work, we develop novel scalable algorithms to tackle this issue, where we selectively decide the regions where the interface thickness has become comparable to the length scale of interfacial features and reduce $Cn$ only in those regions. We refer to this approach as “local Cahn". The current state–of–the–art approaches for CHNS simulation relies on a constant Cahn number throughout the domain

This approach balances the accurate capture of relevant physics and the controlled increase in mesh resolution, which is of utmost importance in tackling engineering problems.

§.§ Identification of key features

As we have discussed earlier, identifying key features such as thin filaments, sheets, and small drops is quintessential for accurate physics modeling.   sec:LocalCahn discusses the importance of accurately detecting these key features and decreasing the Cahn locally.

Decreasing $Cn$ at all interfacial regions can lead to intractable problems due to increased mesh resolution.

We now present algorithms to detect such regions of interest.
We first highlight the salient features of the algorithm on a uniform mesh before discussing the parallel deployment on octree meshes.

§.§.§ Uniform mesh case:

Our algorithm is inspired by the classic image processing idea of erosion and dilation.  fig:schematic briefly shows the algorithmic idea behind the approach for a drop ( fig:Drop) and an elongated filament ( fig:filament).

We perform a series of following steps to achieve the desired goal:

* The phases in CHNS are described by a continuously changing variable $\phi$ ranging from -1 to 1. As a first step, we convert the image from continuous $\phi$ to a binary 0/1 representation $\phi_{BW}$ by thresholding it at a certain value ($\delta$): (denoted by $\mathcal{T}(\phi)$ in  fig:schematic).

\begin{equation*}
\phi_{BW} = \begin{cases}
1 & \quad \phi \leq  \delta, \\
0 & \quad \phi > \delta.
\end{cases}
\end{equation*}

The value of $\delta$ is decided such that the immersed phase gets the value of 1 and the bulk phase gets 0. We

used a value of $\delta$ of 0.8 or -0.8 depending on the value of immersed phase.

* We then perform a series of erosion steps to shrink the region of immersed ($\phi_{BW} = 1$) phase. (denoted by $\mathcal{E}(\phi)$ in  fig:schematic).

* Subsequently,  a series of dilation steps are performed to expand the region of $\phi_{BW}$ (denoted by $\mathcal{D}(\phi)$ in  fig:schematic). We note that we need to perform more dilation than erosion to incorporate the shrinking of the interface encountered in the first step. Our experiments showed that 3 – 4 extra dilation steps are typically sufficient.

* Finally we compute the difference of the pixels after the dilation ($\phi_{BWd}$) and $\phi_{BW}$ (denoted by $\mathcal{S}(\phi)$ in  fig:schematic). The pixels marked as 1 in $\phi_{BW}$ and 0 on $\phi_{BWd}$ constitute the region of interest.

* As the final step, we decrease the Cahn in the regions identified in the previous step and refine the interface region ($|\phi| < \delta^{*}$) with the appropriate resolution.

§.§.§ Challenges with octree mesh:

The above steps work well on a uniform mesh with accessible neighbors; however, performing these steps in a parallel distributed setting on adaptive meshes poses the following challenges:

* The octree–based mesh has elements with non– uniform sizes.

* CG-based finite element (FEM) data–structures may not have neighbor information. This is an issue with unstructured mesh, as elements can have a varying number of neighbors with no plausible upper limits.

* Hanging nodes involve the interpolation of values from parent to child elements. So, the $\phi_{BW}$ value will not have only two distinct values, 0 and 1, but can have arbitrary values between 0 and 1, depending on the type of hanging nodes – edge hanging or face hanging in 3D.

§.§.§ Parallel deployment on octree meshes

We extend the idea of the approach illustrated in  sec:UniformDilation to the octree meshes in parallel settings, albeit with some modifications.

* The first modification involves changing the limits of $\phi_{BW}$. We modify $\phi_{BW}$ as:

\begin{equation}
\phi_{BW,o} = \begin{cases}
1 & \quad \phi \leq  \delta, \\
-1 & \quad \phi > \delta.
\end{cases}
\label{eq:thresholdOctree}
\end{equation}

The above modification in the definition of $\phi_{BW,o}$ is purely a mathematical convenience in detecting the interface elements and helps in subsequent erosion and dilation algorithms.

* Interface elements: An element contains interface when:

\begin{equation}
\lvert \sum_{i=1}^{nodes} \phi_{BW,o} \rvert \neq nodes
\label{eq:interface}
\end{equation}

* Erosion: Erosion step involves visiting the interfacial element ( eq:interface) and making all the nodal values to be -1.

* Dilation: Dilation step involves the reverse of the erosion step. It involves visiting the interface element and making all the nodal values +1.

We note that erosion and dilation involve the call to  rather than . So, the algorithm can proceed without any issue, even if multiple elements try to write to the same nodal location.

With these modifications, we are ready to define the key steps associated with identifying the features of interest in octree-based meshes.  alg:localCahnIdentifier briefs the major steps involved in identifying the local regions of interest and locally decreasing the $Cn$ number. As the first step, we begin by converting the continuous phase field variable $\phi$ to $\phi_{BW,o}$. We then impose erosion and dilation as a series of   operations.( alg:erodeDilate) [  operations are at the heart of FEM computations and has been shown to have excellent scaling  [Ishii et al., 2019, Saurabh et al., 2021]. Each   operation involves a single pass over all the local elements with associated ghost exchange. We note that these ghost exchange communication is typically overlapped with computation  [Ishii et al., 2019, Saurabh et al., 2021, Fernando et al., 2017, Burstedde et al., 2011].] To perform each erosion step, we make a single pass over the elements in a distributed fashion, with each processor looping over its own set of local elements. When the element containing the interface is visited, the erosion step is triggered, which involves converting all the nodal values of that element to -1. As previously mentioned, the octree mesh contains elements at various levels of refinement. To address this challenge, we define a base level $b_l$, typically the finest level in the octree mesh. Additionally, we maintain a counter on how often the element has been visited for erosion. An erosion is triggered when the counter equals the $(b_l - l)$, where $l$ is the current level of the octree mesh. For example, an element two times coarser than the finest mesh, i.e., $b_l - l = 2$, will have to wait for two iterations before erosion can finally be triggered. This gives a way to balance the difference in levels of octree mesh. Similarly, as before, we then perform a series of dilation steps to follow the erosion steps. The dilation is also triggered similarly to account for the difference in octree levels. We perform a larger number of dilation steps compared to the erosion steps to compensate for the thresholding encountered during $\phi_{BW,o}$, thus ensuring that no region that is a part of a larger region gets marked for high refinement, leading to a drastic increase in the problem size. Finally, we make a last pass over the elements and mark it with reduced Cahn based on the following condition:

\begin{equation}
% \begin{split}
\sum_{i=1}^{nodes}\phi_{BW,o}   = nodes   \quad \& \quad
-\sum_{i=1}^{nodes}\phi_{BW,d} = nodes
% \end{split}
\label{eq:condnCahn}
\end{equation}

eq:condnCahn states that we need to reduce Cahn in the element, which has all the nodes marked as +1 after thresholding and -1 after extra dilation. This is analogous to the view in the uniform mesh, except we are looking at element by element rather than nodes. ( alg:elemCahn)

Our approach relies heavily on the number of erosion and dilation steps, a hyper-parameter that controls the identified regions. To prevent a very small region of local $Cn$, we further perform erosion and dilation on the elemental $Cn$ vector ( alg:dilateCahn). This is done by first converting the elemental $Cn$ vector to the nodal $Cn$ vector and performing the subsequent erosion and dilation steps on the nodal vectors, as previously discussed. The erosion steps ensure that all the very small zones of reduced $Cn$ are removed. Similar to erosion, we perform more dilation steps than erosion. The dilation steps pad the region of local $Cn$ elements. These padded elements serve as a potential candidate that would contain the key features of interest till the next level of local $Cn$ identification is performed. We note that padding the elements does not increase refinement as the refinement is only performed at the interface. Even if the pure phase region (away from the interface) gets the decreased $Cn$, no refinement is triggered as we refine only near the interface. Once we have identified the elements for reducing the local $Cn$, we then proceed to refine the regions of the interface marked by $|\phi|<\delta$. We note that the mesh is refined only in the region near the interface, even on the elements with reduced Local $Cn$.

The overall complexity of the proposed approach still scales as $\mathcal{O}(N)$ per dilation and erosion step, which is of the same complexity as with the uniform case. Thus, even without neighbor information, we are able to achieve the same asymptotic complexity. Although we only demonstrate the algorithm deployment for octree-based meshes, this can also be extended for any other class of unstructured meshes. The algorithm proposed relies on efficient   operations, which also form the core kernel for FEM operation. Thus, it can be applied to any unstructured FEM code without any major changes to the framework and does not rely on any additional data–structure like neighbors information, which can be difficult to get in the case of unstructured meshes. Furthermore, although we have demonstrated the applicability of the algorithm to identify the regions with only two levels of local $Cn$, it can be easily extended to multi-level $Cn$. Each level of $Cn$ will have its own set of numbers of erosion and dilation steps.

§.§ Refinement & Intergrid transfer

§.§.§ Adaptive refinement algorithms

As the interface moves and becomes more complex over time, certain
domain regions require greater resolution to maintain accuracy.
The needed level of detail can also change suddenly due to a
large timestep.

We enforce the adaptive mesh resolution in a scalable manner with an
octree-based approach. Adaptive refinement in octrees amounts to
substituting sets of leaves by their ancestors or descendants. We adopt
similar principles as state-of-the-art octree meshing libraries
[Burstedde et al., 2011, Isaac et al., 2015, Ishii et al., 2019, Saurabh et al., 2021].
Namely, linearization, 2:1-balancing, supporting incomplete octrees and
deducing element-to-nodal relations hierarchically on the fly. This
paper focuses on the aspects of octrees relevant to refining by
multiple levels at once. Due to the drastic changes in the
interface, the transformed octree may have leaves that are not direct
children of leaves in the input tree but descendants at an arbitrarily
greater depth.

Refinement in serial

alg:multi-level-refine takes a sorted leaf set and replaces each leaf with its
descendants at the level specified, in sorted order. This algorithm is
elegant in that, even without pointers marking the leaves of subtrees,
the octree traversal passes through the input leaves only once.

The traversal begins at the root. If the current input octant needs
refining below the current traversal depth, then child subtrees are
traversed; otherwise, the current subtree root is emitted. After
processing the descendant subtree of the current octant, the input
pointers are advanced. The next subtree is traversed only if it overlaps
the next input octant.

Note that $\textrm{oct}_\textrm{in}$, $\textrm{level}_\textrm{in}$, and $\textrm{oct}_\textrm{out}$
are passed by reference,
while $R$ and $\textrm{sfc}$ are passed by value.

Refine: Replace leaf octants by sorted descendants

[1]
List of leaf octants and desired levels, with pointers $\textrm{oct}_\textrm{in}$ and  $\textrm{level}_\textrm{in}$.
Mutable list to capture output octants, $\textrm{oct}_\textrm{out}$.
Root of the subtree to traverse, $R$, and its SFC orientation, $\textrm{sfc}$.
$\textrm{oct}_\textrm{out}$ contains the descendants of
$\textrm{oct}_\textrm{in}$ at $\textrm{level}_\textrm{in}$.

not $R$ overlaps $\textrm{oct}_\textrm{in}[0]$
$\textrm{level}(R) < \textrm{level}_\textrm{in}[0]$
$0 \leq c < 2^\textrm{dim}$
refine($\textrm{oct}_\textrm{in}$, $\textrm{level}_\textrm{in}$, $\textrm{oct}_\textrm{out}$, $\textrm{child}(R, \textrm{sfc}, c)$, $\textrm{curve}(\textrm{sfc}, c)$)
push($\textrm{oct}_\textrm{out}, R$)
$\textrm{oct}_\textrm{in} < \textrm{end}$ and $\textrm{oct}_\textrm{in} = R$

Obtaining a refined set of octants can be done independently on each process,
so we do not discuss a parallel version.
However, once the refinement is completed, the 2:1-balance condition must be restored after the refinement.

Refining octants near the domain boundary can create void descendants,
either directly due to subdivision or as an aftereffect of 2:1-balancing.
Void descendants of boundary-intercepted octants need to be discarded.

Coarsening in serial

Analogously to refinement, the octree can be coarsened
by various levels throughout the domain.
An integer is attributed for the coarsest acceptable level to which each octant can be promoted.
Unlike refinement, the decision to coarsen a subtree depends on the inputs of all descendants.
Consensus is defined based on two requirements:
An ancestor $A$ of an input leaf is output if (i) no descendants of $A$ in the input require $A$ to be refined and if (ii) the same cannot be said of the parent of $A$.

alg:multi-level-coarsen traverses the tree while iterating in order over the input—just
as for refinement—but here, the output is pushed and popped at every subtree,
rather than the leaves only.
Per requirement (i), the inputs within the current subtree must be read before deciding
whether the root of the subtree is too coarse.
Therefore, the root is emitted in the post-order action only after checking all descendants.
On the other hand, by (ii), the root of the current subtree could be too fine,
depending on inputs outside the subtree.
Therefore, we leave it to the parent to retract emitted roots of child subtrees if needed.

Coarsen: Replace leaf octants by sorted ancestors

[1]
List of leaf octants and desired levels, with pointers $\textrm{oct}_\textrm{in}$ and  $\textrm{level}_\textrm{in}$.
Mutable list to capture output octants, $\textrm{oct}_\textrm{out}$.
Root of the subtree to traverse, $R$, and its SFC orientation, $\textrm{sfc}$.
$\textrm{oct}_\textrm{out}$ contains the ancestors of
$\textrm{oct}_\textrm{in}$ at $\textrm{level}_\textrm{in}$,
subject to the least amount of coarsening voted by any descendant.

coarsen_to $\leftarrow$ 0
not $R$ overlaps $\textrm{oct}_\textrm{in}$[0]
level($R$) $<$ level($\textrm{oct}_\textrm{in}$[0])
pre_size $\leftarrow$ size($\textrm{oct}_\textrm{out}$)
Remember where to emit parent
$c$ in range(0,pow2(dim))
$L_c$ $\leftarrow$
coarsen($\textrm{oct}_\textrm{in}$, $\textrm{level}_\textrm{in}$, $\textrm{oct}_\textrm{out}$, child($R$, sfc, $c$), curve(sfc, $c$))
coarsen_to $\leftarrow$ max(coarsen_to, $L_c$)
coarsen_to $\leq$ level($R$)
repeat(pop($\textrm{oct}_\textrm{out}$), size($\textrm{oct}_\textrm{out}$) - pre_size)
Undo child emit
push($\textrm{oct}_\textrm{out}$, $R$)
Emit parent octant
push($\textrm{oct}_\textrm{out}$, $R$)
coarsen_to $\leftarrow$ $\textrm{level}_\textrm{in}$[0]
$\textrm{oct}_\textrm{in}$.nonempty() and $\textrm{oct}_\textrm{in}$[0] = $R$

In [Burstedde et al., 2011] serial coarsening and refinement
algorithms are discussed, including multi-level variants. Compared to our algorithm
(<ref>), theirs accepts a slightly more general
coarsening criterion, expressed as a Boolean callback function
that is dynamically evaluated on every subtree. In
contrast, our algorithm is tailored to the case that the application
declares up front the maximum change in level that should be allowed for
each octant. Our algorithm, however, removes the requirement that
subtrees be complete.
Otherwise, both algorithms have an equivalent structure, featuring a post-order tree traversal,
a single pass through the input, and writing outputs through a pure stack interface (push/pop).

Next, we present algorithms to coarsen in parallel and
transfer functions across multiple levels of refinement in parallel.

Coarsening in parallel

For the same reasons that coarsening in serial is more complicated than refinement,
coarsening in parallel requires extra communication to guarantee consensus.

ParCoarsen: Replace distributed leafs by partitioned and sorted ancestors

[1]
List of leaf octants and desired levels, $\textrm{oct}_\textrm{in}$ and  $\textrm{level}_\textrm{in}$.
MPI communicator $\textrm{comm}$.
Returns the ancestors of
$\textrm{oct}_\textrm{in}$ at the consensus of $\textrm{level}_\textrm{in}$, with duplicates removed.

first coarsening pass
tentative = coarsen($\textrm{oct}_\textrm{in}$,  $\textrm{level}_\textrm{in}$)

exchange tentative coarse octants at partition endpoints
$r \leftarrow$ rank(comm)
head_tail[1] $\leftarrow$ $\{$ tentative[0], tentative[-1] $\}$
head_tail[2] $\leftarrow$ send_recv(head_tail[1], to:$r-1$, from:$r+1$, comm)
head_tail[0] $\leftarrow$ send_recv(head_tail[1], to:$r+1$, from:$r-1$, comm)

find overlapped inputs
overlap_prev = [],  overlap_next = []
level(tail[0]) $\leq$ level(head[1])
overlap_prev = [$x$ if overlaps($x$, tail[0]) for $x$ in $\textrm{oct}_\textrm{in}$]
level(head[2]) $<$ level(tail[1])
overlap_next = [$x$ if overlaps($x$, head[2]) for $x$ in $\textrm{oct}_\textrm{in}$]

repartition overlapped inputs (omitted: also repartition levels)
$\textrm{oct}_\textrm{pre} \leftarrow$ send_recv(overlap_next, to:$r+1$, from:$r-1$, comm)
$\textrm{oct}_\textrm{post} \leftarrow$ send_recv(overlap_prev, to:$r-1$, from:$r+1$, comm)
$\textrm{oct}_\textrm{in} \leftarrow \textrm{oct}_\textrm{in}$ - overlap_prev - overlap_next + $\textrm{oct}_\textrm{pre}$ + $\textrm{oct}_\textrm{post}$

second coarsening pass

coarsen($\textrm{oct}_\textrm{in}$,  $\textrm{level}_\textrm{in}$)

alg:par-coarsen proceeds as follows.
The globally sorted input octree undergoes a local coarsening pass
resulting in a locally linearized incomplete tree.
The results are tentative, as the output octants achieve consensus locally on a process,
but may not represent the global consensus.
Furthermore, coarse octants might be duplicated even if the results agree across processes.
The input must be slightly repartitioned before the second coarsening pass to prevent conflicts.
Note that the required partitioning resolution
for multi-level coarsening is known only after the first coarsening attempt.

After the first local coarsening pass, the next step is finding the overlap region between consecutive processes.
Ideally, the tentative coarse octants do not overlap between two consecutive processes $p_i$ and $p_{i+1}$.
If so, the inputs on processes $p_{j \leq i}$ are independent from those on processes $p_{j > i}$.
Otherwise, the worst case for two processes is for a single coarse octant
on one process to overlap one or more octants on the other process.

Note that incomplete trees require special handling compared to complete trees.
For complete trees, conflicts can be resolved by globally sorting the tentative coarse octants,
then locally splitting any coarse octants until all overlaps are removed.
For incomplete trees, the same splitting procedure could produce octants in void regions,
leading to unnecessary fill-in.
There are several possible solutions.

Option one is to reapply the domain test after splitting coarse octants.
The advantage here is that no further communication is needed after the partitioning step.
But, the cost of the domain test may not be optimal for pinpoint tests at partition endpoints;
rather we treat the cost as amortized over bulk queries during startup and refinement.
Also, the domain test is completely redundant for coarsening
because the input octree already contains the needed information.

Option two is to send each split octant back to the process(es) whose input could overlap
that octant and have the originating process(es) figure out whether the octant is void or not.
Equivalently, the overlapped input could be requested, with void octants discarded on the same process that performed the splits.
(Multiple processes might need to be consulted, but only if the input grain size is very small).
Either way, the disadvantage is that a round of ping-pong messages
is needed after the partitioning step.

Option three is to repartition the input directly toward the coarsest contender of each conflict.
This solution cuts out the lag time needed to fully partition and split tentative coarse octants
before resolving conflicts.
As a result, communication is limited to an exchange of tentative coarse octants at the partition endpoints and a repartition of overlapped inputs.
The new partition makes it just as if there were no overlaps to begin with.
Coarsening can be finished independently of other processes,
and only the exchanged regions need to be recomputed after the exchange.

We have taken the approach of option three.
Although option two can perform partitioning at a finer granularity,
it does not promise to restore load balance any better than option three.
We consider proper load balancing a separate step in any case,
especially since the subsequent restoration of 2:1-balance tends to affect the load.
Therefore, we choose option three, which both avoids redundant domain tests
and streamlines communication.

In the rare case that a tentative octant is so aggressively coarsened
that it overlaps multiple remote partitions, the initial coarse octant exchange
can be considered the first step of a distributed exponential search
to find the processes containing the region of overlap.
The details of this adjusted algorithm are omitted for brevity.

Proof of worst-case overlap:
It is claimed above that, after local coarsening, the worst case overlap for two processes
is that one or more octants are overlapped by a single remote coarse octant.
The proof hinges on the assumption that the input to local coarsening is globally sorted:
Suppose process $p$ precedes process $q$, and last two octants of $p$ are $s_p < t_p$ while the first two octants of $q$ are $s_q < t_q$.
Now suppose, contrary to the claim, that $s_p$ and $s_q$ overlap, and $t_p$ and $t_q$ overlap.
It must be that $s_q$ does not overlap $t_p$; for otherwise $s_q$ would be an ancestor of $t_p$ and thus would overlap $t_q$, which is impossible because $s_q$ and $t_q$ belong to the same locally linearized tree.
Since $s_q$ overlaps $s_p$ and $s_p < t_p$, the pre- or post-order implies $s_q$ precedes $t_p$ (still without overlapping).
But, this implies that the inputs descending from $s_q$ precede the inputs descending from $t_p$.
In other words, the input is not globally sorted, contradicting the assumption that it is.

§.§.§ Distributed multi-level inter-grid transfer

After updating the resolution of the octree mesh on a new time step
(and restoring 2:1-balance),
data from the previous grid must be transferred to the resolution and partition of the new grid.

Octree-based meshing frameworks typically support inter-grid transfers only between successive levels [Sundar et al., 2012].
This constraint simplifies the implementation because it means overlapping cells differ in size by a factor of at most two.
However, for the problems with moving interfaces and suddenly changing refinement patterns,
constructing intermediate grids one level at a time creates significant overhead.
This section describes our methods to transfer data after refining or coarsening an arbitrary number of levels.

Serial algorithms

Transferring coarse-to-fine cell-centered data is the simplest case. We
traverse the octants of both grids in the SFC order simultaneously. Once
a leaf is reached in the coarse grid, all cells overlapped in the fine
grid are accessible. The value in the coarse cell is imposed on all the
overlapped fine cells.

Coarse-to-fine interpolation of node-centered data begins similarly.
Octants of both grids are traversed together. Then, at a coarse leaf,
nodal values on all overlapped fine cells can be interpolated from the
coarse cell nodes. Using mesh-free techniques [Ishii et al., 2019],
the coordinates of the relevant nodes are available in a contiguous list
at this point in the traversal. The transfer is done by evaluating the field
in the coarse element at each of the fine nodes. If neighboring subtrees redundantly interpolate a node, the results are consistent in each instance; the final value is arbitrarily picked from one of the instances.

Fine-to-coarse transfers follow essentially the same structure, except
mirrored. As multiple fine unknowns are aggregated into fewer coarse
unknowns at a coarse leaf, the exact aggregation function depends on the
differential operator and the type of transfer. Cell-centered values
might be averaged. Node-centered values might be simply injected or
otherwise weighted by the transpose of the interpolation operator.

The advantage of using mesh-free grid traversals here is that the nodes
need not be initially arranged in any particular order.
The only requirement is that nodal values are tagged by their unique location code key.
The keys make it possible to route the appropriate nodes to each subtree during the traversal.
This property is needed by our distributed memory approach
of sending nodes—with their keys—to interested processes in a detached manner.

Parallel algorithm

The serial grid transfer algorithms can be applied in distributed memory
after the partitions are somehow matched between the two grids.

Our inter-grid transfer algorithms are structured in four steps:
First, search for grid-grid overlaps in the splitter table.
Second, detach and send coarse element nodes to the fine grid.
Third, execute serial inter-grid interpolation or aggregation.
Fourth, if coarsening, receive coarse element data from the fine grid and attach it.

It is possible to communicate either before or after the local inter-grid transfer.
We choose the order that communicates the coarse data.
This helps reduce the total message size.
Also, since the number of fine elements dominates the local inter-grid transfers,
this strategy uses the fine partition to keep the computational workload balanced.

Extension of rank search from octants to overlaps

The function \(\textrm{rank}_{X^\sim}(e) \equiv |x \in X : x \sim e|\)
on a totally ordered set \((Y,\sim), X \subset Y\) guarantees that rank
is non-decreasing on a sorted subset of \(Y\). This property suggests the desired key can be found by narrowing down an interval, one of the main ideas behind the binary search. To apply binary search algorithms to find overlaps between two arbitrary leaf sets, extending the notion of rank from octants to a kind of overlapping region is useful.
Specifically, from the total ordering of octants, this section constructs a new total ordering over equivalence classes whose members share a common ancestor.

This section also uses the properties of hierarchical ordering and leaf

* Hierarchical octant ordering: Let \(a, x, y\) be octants such
that \(a\) is an ancestor of \(x\) but not of \(y\). Then
\(y<a \Leftrightarrow y<x\). Equivalently,
\(a<y \Leftrightarrow x<y\). (Both preorders and postorders have this

* Leaf set: If \(G\) is a grid and \(L=\textrm{leafs}(G)\), then
\(\forall x \in L\), the only ancestor of \(x\) in \(L\) is \(x\)

Let \(G\) and \(H\) be two leaf sets. Define \(L=G \cup H\).

For any \(x,y \in L\), define the equivalence relation \(x \frown y\) to
mean there is \(a \in L\) that is an ancestor of both \(x\) and \(y\).
Reflexive: \(\forall x \in L\), \(x\) is an ancestor of \(x\).
Symmetric: The definition only uses the expression of both
\(x\) and \(y\), which is symmetric. Transitive: Let
\(x,y,z \in L\) with \(x \frown y\) via ancestor \(u \in L\), while
\(y \frown z\) via ancestor \(v \in L\). Since \(u\) and \(v\) are both
ancestors of \(y\), one of \(u\) or \(v\) is an ancestor of the other,
and thus of \(x\) and \(z\) as well.

Each equivalence class of \(\frown\) has a coarsest representative. This
can be seen by cumulatively repeating the transitivity argument over the
whole class, producing an ancestor shared by all members of the class.

The picture is especially clear-cut for two leaf sets, as \(L\) is
defined above. For two leaf sets, the coarsest element belongs to one
grid, and any overlapped finer elements are disjoint from each other and
belong to the other grid. Additionally, for two leaf sets, it is
immediately apparent whether elements from different grids belong to the
same class: If so, they overlap. The common ancestor of two class
members comes from one of the leaf sets; thus it must equal the member
from the same leaf set.

Next is a proof that the equivalence classes are totally ordered. Take a
total ordering \(\leq\) that has the hierarchical property (such as a
pre- or post-order SFC traversal) and define a new relation: For any
\(x,y \in L\), define \(x \sqsubseteq y\) to mean that \(x < y\) or
\(x \frown y\). This is a quasiorder whose equivalence relation is
\(\frown\) (see antisymmetry below). Thus the quotient of
\(\sqsubseteq\) partially orders the equivalence classes of \(\frown\)
by \(\leq\), which is a total order. Reflexive: \(\frown\) is
reflexive. Antisymmetric: If \(x \sqsubseteq y\) and
\(y \sqsubseteq x\), then \(x \frown y\). Otherwise, using the symmetry of
\(\frown\), it must be that \(x<y\) and also \(y<x\), which is
impossible. Transitive: Let \(x,y,z \in L\) with
\(x \sqsubseteq y\) and \(y \sqsubseteq z\). There are four cases. (i)
\(x<y\) and \(y<z\). Then \(x<z\). (ii) \(x \frown y\) and
\(y \frown z\). Then \(x \frown z\), by transitivity of \(\frown\).
(iii) \(x<y\) and \(y \frown z\), via ancestor \(v \in L\). If \(v\) is
an ancestor of \(x\), then \(x \frown z\) via \(v\). Else, by the
hierarchy property, \(x<y \leftrightarrow x<v \leftrightarrow x<z\).
(iv) \(x \frown y\), via ancestor \(u \in L\), and \(y < z\). If \(u\)
is an ancestor of \(z\), then \(x \frown z\) via \(u\). Else, by the
hierarchy property, \(y<z \leftrightarrow u<z \leftrightarrow x<z\).

The inclusive rank function, \(\textrm{rank}_{X \sqsubseteq}(e)\), maps
all members of a class to the same rank. This also holds for the
exclusive rank, \(\textrm{rank}_{X \sqsubset}(e)\). (The irreflexive
kernel of \(\sqsubseteq\) is \(x \sqsubset y\), which means that
\(x < y\) without \(x \frown y\)). The inclusive rank of \(e\) is the
exclusive rank of the next class after that of \(e\).

Finally, it is apparent that searching for grid overlaps is nearly the
same as for \(<\) rank. A leaf set sorted by \(\leq\) is also
sorted by \(\sqsubseteq\). Double-rank search then only requires
evaluating \(x \sqsubseteq y\) for \(x\) and \(y\) from distinct leaf
sets. This is possible since, for two grids, \(x \frown y\) can be
evaluated between distinct leaf sets.

Intervals of components

Denote every partition's first and last leaves in grid \(G\) by the
arrays \(G^-\) and \(G^+\), such that on process \(r\), any leaf
\(g_r\) would fall in the interval \(G^-_r \leq g_r \leq G^+_r\). Similarly,
denote arrays \(H^-\) and \(H^+\) for grid \(H\). Each interval of
octants maps to an interval of overlap regions with the total ordering
\(\sqsubseteq\). A \(\sqsubseteq\)-interval \(G^-_p \dots G^+_p\)
intersects \(H^-_q \dots H^+_q\) if, and only if, both
\(G^-_p \sqsubseteq H^+_q\) and \(H^-_q \sqsubseteq G^+_p\).

These observations show that the range of processes \(q\) overlapped in
\(H\) by the interval in \(G\) owned by process \(p\) can be found with
two binary searches:
\(\textrm{rank}_{H^+ \sqsubset}(G^-_p) \leq q < \textrm{rank}_{H^- \sqsubseteq}(G^+_p)\)
. Because the searches only involve partition endpoints, not
process-local data, the intersections are detected consistently across

Given one of the processes \(q\) and its partition endpoints in \(H\),
the range of overlapping local octants \(G_p[i]\) is found by searching:
\(\textrm{rank}_{G_p \sqsubset}(H^-_q) \leq i < \textrm{rank}_{G_p \sqsubseteq}(H^+_q)\)
. In practice, one range starts no earlier than one octant before the
end of the preceding range. Then, the coarse elements are selected from each range of octants.

Cell-centered data is already in the correct order and can be sent and received
once the partition overlap search is complete.
Node-centered data, on the other hand, must be clipped out from the nodal ordering
separately for each destination process before it can be communicated.

Detaching and attaching coarse node-centered data

The process-to-element overlap range is combined with the element-to-node relation
to find the remote interpolation source nodes for each process.
The nodes are not detached element-by-element.
That would inflate the message size more than necessary since many elements destined for the same process neighbor each other and share nodes.
Instead, for each process, we loop just over the range of elements destined for that process, setting a flag for all nodes on those elements.
After the loop, the flagged nodes are gathered and ready to be sent in a contiguous message.
This is repeated for all destination processes.

The process of attaching node-centered data onto coarse cells after aggregation is similar.
The flags from the elemental loops are reused from the detachment phase.
The earlier gather movement is reversed, as the packed nodes are scattered to the positions of the flags.

§.§.§ Scalability Checkup

Preliminary runs on up to 30K cores alerted us to some suboptimal usage
of MPI collective communication in the previous versions of the code.

Sorting octree keys (Allreduce and Alltoallv)

Sorting octree keys in distributed memory is a building block of our
meshing routines, such as for repartitioning, 2:1-balancing, and
enumeration of nonhanging nodes. In our old implementation of the
distributed octree sorting algorithm, we experienced poor scalability of
splitter search, Allreduce, and Alltoall exchange. In the new version,
we adopt a hierarchical \(k\)-way staged communication pattern, similar
to other hypercube exchange sorting algorithms
[Sundar et al., 2013, Fernando et al., 2017]. In the \(k\)-way
hierarchical scheme, the number of (super)partitions is kept below a
constant \(k\) for each of \(\mathcal{O}(\log_k(p))\) stages. (By
default \(k=128\), so at most, three stages are required up to 2M processes).
Dissociating the partitions from process count reduces storage cost of
splitter selection from \(\mathcal{O}(p)\) to \(\mathcal{O}(k)\), and
data transfer of Allreduce from \(\mathcal{O}(p)\) to
\(\mathcal{O}(k \log_k(p))\). Performing the Alltoall exchange in stages
is also a standard defense against network congestion.

Hierarchical communication (Comm_split)

During the setup phase of the \(k\)-way hierarchical exchange, the
global communicator is subdivided recursively with
, such that the final stage involves \(k\) or
fewer processes. Splitting a communicator is a costly global operation.
Doing it repeatedly poses a serious threat to scalability. Since the
arguments to  are not dependent on the data to
be sorted, we memoize the sequence of communicators in an MPI user
cache attribute attached to the root communicator. On later executions
of distributed octree sort, the saved sequence of communicators is
recalled without extra splitting.

Remote processing (Sparse Alltoallv)

Our nodal enumeration algorithm contains an outsourcing pattern, whereby
candidate nodes are sorted globally, checked for duplication and
hangingness on remote processes, and sent back to the owners of the
originating elements. This algorithm, although not necessarily optimal,
is simple and has scaled adequately, following a communication fixup.

The last step of the algorithm uses a concept of return address
to return nodes to originating processes. Due to the high-locality
heuristic of SFC sorted orders, the number of destination processes is
actually sparse. At the same time, the return addresses are treated as
if arbitrary. Our first implementation used a raw
collective to obtain the receive counts at the destinations.
tends to produce heavy network congestion and
requires time to populate the send count array, which is linear in
process count. We noticed that this step of the algorithm had low
overhead up to 28K cores, at which point scaling halted and the overhead
blew up 15x from 28K to 56K cores. To fix this we adopted the NBX sparse
exchange algorithm from [Hoefler et al., 2010], which uses
non-blocking collectives and point-to-point messages to eliminate any
\(\Omega(p)\) primitives like .

§.§ Matrix & Vector assembly

In this work, we interface our framework with  data–structure to perform linear algebra operations. Specifically, we store the matrix in the form of block storage MATMPIBAIJ. This format has been demonstrated to be much more efficient than the non-block version MATMPIAIJ, specifically for the multi–dof system. The block size, here, refers to the number of DOFs per node.

This storage format, although ideal for the global format, is not suitable for local assembly. This is primarily because of the fact that any operator of the form $\mathcal{L}(dof_{i},dof_{j})$ tends to have the strided memory access. The top row  fig:vector shows the memory view for the global vector. [We note that vector assembly will have only one loop over basis function, i.e., they will be of the form , whereas the matrix assembly involves two loops and are of form .] Any loop over $(dof_{i})$ in vector assembly will write to memory in strided manner. For example, ($i = 0$) will visit 0,2,4,6 and  ($i = 1$) will write to 1,3,5,7. This is not ideal for optimal performance.

In order to circumvent this issue, we perform a sequence of the following operations:

* $zip$: We perform the $zip$ operations for the DOFs, where similar DOFs are arranged contiguously in memory. This is done by making just a single pass over the elemental vector.

* Perform elemental vector assembly.

* $unzip$: Finally, we perform the $unzip$ operations to revert back to the original global format.

Similar to the vector, matrix assembly faces the same issue. The assembly of each DOF writes in the strided assembly. For example, in case of 2D with 2 dof system, $\mathcal{L}(dof_0,dof_0)$ would write to $0,2,4,6,16,18,20,22,32,34,36,38,42,44,46,48$ [considering flattened row major addressing]. Clearly, we can see that there is a lot of strided access while writing the values in the local vector. We perform a similar $zip$ operation to pack all the DOFs together. A good way to look the $zip$ operation is that, after the $zip$ operation, the matrix size is $16\times4$, where ($dof_0,dof_0$) writes to $(0\cdots15)$, ($dof_0,dof_1$) to $(16\cdots31)$, ($dof_1,dof_0$) to $(32\cdots47)$ and ($dof_1,dof_1$) to $(48\cdots63)$.

In the case of matrix assembly, an explicit $zip$ operation is never performed. We begin each elemental assembly with a matrix filled with 0 entries. So, we perform only $unzip$ operations once the local matrix is completely assembled and written back to the global matrix.

Additionally, we extended the idea proposed by Saurabh et al., 2021 and exposed each matrix assembly operator as a matrix–matrix multiplication and each vector assembly operator as matrix–vector multiplication. This way of expressing the FEM operator allows us to leverage the benefits of vendor–optimized GEMM and GEMV kernels. We note the $zip$ and $unzip$ operation is quintessential to express the operator as matrix–matrix or matrix–vector multiplication. This also ensures the portability of the code across different hardware platforms. It forms a middle ground between the portability and explicit vectorization of the FEM operator.

§.§ Challenges with parallel settings

The parallel deployment of multiphase simulation comes up with its own set of challenges. As the simulation progress in time, the interface blooms, break into droplets, and forms complex regions and shapes. This leads to a significant increase in the mesh resolution as the time and simulation progress. We showcase one such simulation of primary jet atomization later in  sec:Atomization. With the increase in the mesh resolution, it is desirable to increase the number of processes accordingly.

We dump the files checkpoint files at frequent intervals. Once the checkpoint is reloaded, we can load it either from the same or larger number of processes. Specific care must be taken while loading with increased processor count. When the processor count is increased,  is split into  and  communicators. The checkpoint is loaded with  communicators, and the mesh is created only within the  communicator. The  communicator has the size of the processes that dump the checkpoint files. Once a repartition or remeshing is triggered, the new mesh is then redistributed among all the available processes. The  communicators become  once the repartition or remeshing function is called.

§ RESULTS & DISCUSSION

Show the performance improvement for a simple bubble rise case.

§.§ Region Identification:   scaling

The proposed algorithm for erosion and dilation scheme relies solely on the   operations, which involve performing the local elemental traversal with associated ghost exchange. Therefore, we show the performance of   operations of our framework for linear basis functions.  fig:mvecscaling shows the strong and weak scaling of   for our framework. For strong scaling, we consider an adaptive mesh of around 13M elements and 13.7 M DOFs. fig:mavecStrongScaling shows the   execution time as a function of increasing number of cores.   execution time decreased from 2.87 s on 224 processes to 0.027 s on 28K processes, resulting in 81% parallel efficiency for 128 fold increase in processor count. For the weak scaling runs, we created grids with a fixed grain size of about 35K elements per core and timed   execution time. The coarsest mesh consists of 981K elements on 28 processes  with 1.02M DOFs for linear and 8.01M DOFs.  fig:mavecWeakScaling shows the    averaged over 100 iterations as a function of the number of cores. A constant execution time would imply ideal weak scaling efficiency. We observed a slowly growing weak-scaled execution time. Overall the time increased from about 1.58 s on 28 cores to 1.9 s on 14 K cores for linear elements, resulting in 82% weak scaling efficiency.

Single processor improvements.

§.§ Application performance scaling

fig:applicationScaling shows the scaling performance of the full framework on TACC  [All simulation was performed using 56 cores per node]. The overall mesh consists of around 700 M mesh elements. The simulation was run for 11 timesteps. We have timed each of the solvers and remeshing separately. Overall, we see a good scaling behavior in all the solvers with NS-solve, resulting in $6.6 \times$ speedup for $8 \times$ increase in the processor from $\mathcal{O}(14K)$ processor to $\mathcal{O}(112K)$ processor. Similarly, PP–solve showed a $5.3 \times$ reduction, whereas VU–solve showed $5.5\times$ and CH–solve showed a $4 \times$ reduction in total solve time. The remeshing cost showed $2.5 \times$ reduction for $4 \times$ increase in processes till $\mathcal{O}(57K)$ processes but then began to grow with further increase in processor. This increased cost in the remeshing needs further investigation.

The variable density PP–solve is the most time-consuming step till the point where remeshing cost begins to dominate. Solving pressure Poisson efficiently, especially with variable coefficients, is still a current area of research. Scalable solvers, like Geometric multigrid (GMG), promise to yield a better solve time but relies on optimized algorithms for creating different mesh hierarchy and   operation. [We have tested the framework with Algebraic multigrid via   and Hypre interface. Our study showed that iterative solver yields a better solve time at a large processor count due to the setup cost associated with AMG. Therefore, we used iterative solvers for the current study.] This is left as future work.

Weak scaling: We do not produce any weak scaling results as it is hard to create an adaptive mesh with a fixed grain size per processor. Arbitrarily refining the elements lead to a much-increased iteration count in during KSP solve and is therefore not a correct metric for the scaling study for complete application problem.

§ APPLICATION TO PRIMARY JET ATOMIZATION

We finally demonstrate the ability of our algorithm to simulate the most resolved simulation of primary jet atomization. The finest resolution in our application problem consists of octree at level 15. This is equivalent to solving it on a 35 trillion grid point – a 64$\times$ larger than the current state–of–the–art.  fig:jet shows the simulation result at 6.3 $\mu s$. We can see that the simulation framework can capture tiny droplets by selectively identifying the key reasons of interest and selectively increasing the resolution.

Simulation of primary jet atomization at 6.3 $\mu s$.

fig:zoomview shows the progressive refinement of the mesh. We can see that the algorithm proposed in  sec:FeatueDetection is able to detect complicated structures like filaments and drops and selectively refine those regions. The overall interface is resolved at level 13 [An octree of level $l$ implies $\Delta x = \Delta y = \Delta z = L/2^{l}$, where L is the length of the domain.] with the key features resolved at level 15.  fig:elemFraction shows the fraction  of the element at different levels. We see that a significant portion of the elements is at levels 13 and 14 ($\approx 25 \%)$, with a maximum fraction at level 15. A rough estimate indicates that resolving the complete interface at level 15 would result in at least an $8 - 10 \times$ increase in the overall mesh count, resulting in approximately $20-25\times$ higher time to solve,  [assuming $\mathcal{O}(NlogN)$ scaling of KSP solve.] yielding the simulation impractical in a reasonable time. We note that level 15 has the maximum element fraction but covers only $0.01\%$ of total volume. This illustrates the importance of adaptivity in resolving multiphase simulations. We defer the more detailed physics discussion for a future paper.

§ RELATED WORK

In this work, we have shown the importance of identifying the key regions of interest and selectively refining those specific regions. The alternative approaches for accompanying the task can be through the route of determining the connected components  [Harrison et al., 2015]. However, the connected component is computationally more expensive and non–trivial to implement compared to the current approach. Additionally, connected component labeling is not enough to identify all the regions of interest, such as the thin filament, which is a part of the big structure (for example, fig:filament). Furthermore, the proposed approach by  Harrison et al., 2015 is limited to octree-based meshes, whereas the current approach in this work can be easily integrated into other numerical libraries with structured and unstructured meshes, as it solely relies on efficient   computations – a key component of scientific computations.

In the context of sharp interface capturing methods, there are advances in techniques, for example, in interface reconstruction [Lu and Tryggvason, 2018, Chiodi, 2020], moment-of-fluid [Jemison et al., 2015], manifold death [Chirco et al., 2022] that also resolve try to resolve small structures. These techniques are, however, limited to sharp interface methods at the moment.

§ CONCLUSION & FUTURE WORK

In this work, we present algorithmic advances in the field of multiphase modeling. We proposed a novel scalable algorithm to identify the key regions of interest that need sufficiently high resolution for capturing the physics accurately. This approach is quintessential for efficient simulation of multiphase simulation in a reasonable time. Further, we modify the octree refinement algorithm to refine to an arbitrary level. We also accelerate the matrix and vector assembly by carefully analyzing the data movement. We showcase the application of these algorithmic advances by simulating the canonical problem of primary jet atomization simulation and push the current state–of–the–art by a factor of 64$\times$. The approach proposed in this work can serve as a reference for the next–generation multiphase modeling for diffuse as well as sharp interface modeling.

In future, we plan to extend these algorithms, by integrating the framework with GPU accelerated linear algebra solvers like Amat [Tran et al., 2022], Magma [Tomov et al., 2009]. Further, we plan to improve the DA scaling and utilize GMG to improve the solve–time, specifically for the variable coefficient pressure Poisson problem. We also plan to integrate the setup with Domain Specific Language (DSL) like Finch [Heisler et al., 2022], Fenics [Alnæs et al., 2015] and extend FEM framework to Finite Volume and Discontinuous Galerkin schemes.

[Fuster et al., 2009]
D. Fuster, A. Bagué, T. Boeck, L. Le Moyne, A. Leboissetier, S. Popinet,
P. Ray, R. Scardovelli, and S. Zaleski, “Simulation of primary atomization
with an octree adaptive mesh refinement and vof method,” International
Journal of Multiphase Flow, vol. 35, no. 6, pp. 550–565, 2009.

[Tomar et al., 2010]
G. Tomar, D. Fuster, S. Zaleski, and S. Popinet, “Multiscale simulations of
primary atomization,” Computers & Fluids, vol. 39, no. 10, pp.
1864–1874, 2010.

[Pairetti et al., 2020]
C. I. Pairetti, S. M. Damián, N. M. Nigro, S. Popinet, and S. Zaleski,
“Mesh resolution effects on primary atomization simulations,”
Atomization and Sprays, vol. 30, no. 12, 2020.

[Herrmann, 2010]
M. Herrmann, “Detailed numerical simulations of the primary atomization of a
turbulent liquid jet in crossflow,” Journal of Engineering for Gas
Turbines and Power, vol. 132, no. 6, 2010.

[Gorokhovski and Herrmann, 2008]
M. Gorokhovski and M. Herrmann, “Modeling primary atomization,” Annu.
Rev. Fluid Mech., vol. 40, pp. 343–366, 2008.

[Lu and Tryggvason, 2018]
J. Lu and G. Tryggvason, “Direct numerical simulations of multifluid flows in
a vertical channel undergoing topology changes,” Physical Review
Fluids, vol. 3, no. 8, p. 084401, 2018.

[Chiodi, 2020]
R. M. Chiodi, Advancement of numerical methods for simulating primary
atomization.1em plus 0.5em minus 0.4emCornell University,

[Jemison et al., 2015]
M. Jemison, M. Sussman, and M. Shashkov, “Filament capturing with the
multimaterial moment-of-fluid method,” Journal of Computational
Physics, vol. 285, pp. 149–172, 2015.

[Chirco et al., 2022]
L. Chirco, J. Maarek, S. Popinet, and S. Zaleski, “Manifold death: a volume of
fluid implementation of controlled topological changes in thin sheets by the
signature method,” Journal of Computational Physics, vol. 467, p.
111468, 2022.

[Saurabh et al., 2021]
K. Saurabh, B. Gao, M. Fernando, S. Xu, M. A. Khanwale, B. Khara, M.-C. Hsu,
A. Krishnamurthy, H. Sundar, and B. Ganapathysubramanian, “Industrial scale
large eddy simulations with adaptive octree meshes using immersogeometric
analysis,” Computers & Mathematics with Applications, vol. 97, pp.
28–44, 2021.

[Saurabh et al., 2021]
K. Saurabh, M. Ishii, M. Fernando, B. Gao, K. Tan, M.-C. Hsu, A. Krishnamurthy,
H. Sundar, and B. Ganapathysubramanian, “Scalable adaptive pde solvers in
arbitrary domains,” in Proceedings of the International Conference for
High Performance Computing, Networking, Storage and Analysis, 2021, pp.

[Fernando et al., 2018]
M. Fernando, D. Neilsen, H. Lim, E. Hirschmann, and H. Sundar, “Massively
parallel simulations of binary black hole intermediate-mass-ratio
inspirals,” arXiv preprint arXiv:1807.06128, 2018.

[Popinet, 2003]
S. Popinet, “Gerris: a tree-based adaptive solver for the incompressible euler
equations in complex geometries,” Journal of computational physics,
vol. 190, no. 2, pp. 572–600, 2003.

[Popinet, 2009]
——, “An accurate adaptive solver for surface-tension-driven interfacial
flows,” Journal of Computational Physics, vol. 228, no. 16, pp.
5838–5866, 2009.

[Bangerth et al., 2012]
W. Bangerth, C. Burstedde, T. Heister, and M. Kronbichler, “Algorithms and
data structures for massively parallel generic adaptive finite element
codes,” ACM Transactions on Mathematical Software (TOMS), vol. 38,
no. 2, pp. 1–28, 2012.

[Khanwale et al., 2021]
M. A. Khanwale, K. Saurabh, M. Ishii, H. Sundar, J. A. Rossmanith
et al., “A projection-based, semi-implicit time-stepping approach for
the cahn-hilliard navier-stokes equations on adaptive octree meshes,”
arXiv preprint arXiv:2107.05123, 2021.

[Ishii et al., 2019]
M. Ishii, M. Fernando, K. Saurabh, B. Khara, B. Ganapathysubramanian, and
H. Sundar, “Solving pdes in space-time: 4d tree-based adaptivity, mesh-free
and matrix-free approaches,” in Proceedings of the International
Conference for High Performance Computing, Networking, Storage and Analysis,
2019, pp. 1–61.

[Fernando et al., 2017]
M. Fernando, D. Duplyakin, and H. Sundar, “Machine and application aware
partitioning for adaptive mesh refinement applications,” in
Proceedings of the 26th International Symposium on High-Performance
Parallel and Distributed Computing, 2017, pp. 231–242.

[Burstedde et al., 2011]
C. Burstedde, L. C. Wilcox, and O. Ghattas, “p4est: Scalable algorithms for
parallel adaptive mesh refinement on forests of octrees,” SIAM Journal
on Scientific Computing, vol. 33, no. 3, pp. 1103–1133, 2011.

[Isaac et al., 2015]
T. Isaac, C. Burstedde, L. C. Wilcox, and O. Ghattas, “Recursive algorithms
for distributed forests of octrees,” SIAM Journal on Scientific
Computing, vol. 37, no. 5, pp. C497–C531, 2015.

[Sundar et al., 2012]
H. Sundar, G. Biros, C. Burstedde, J. Rudi, O. Ghattas, and G. Stadler,
“Parallel geometric-algebraic multigrid on unstructured forests of
octrees,” in SC'12: Proceedings of the International Conference on
High Performance Computing, Networking, Storage and Analysis.1em plus
0.5em minus 0.4emIEEE, 2012, pp. 1–11.

[Sundar et al., 2013]
H. Sundar, D. Malhotra, and G. Biros, “HykSort: a new variant of hypercube
quicksort on distributed memory architectures,” in Proceedings of the
27th international ACM conference on international conference on
supercomputing, 2013, pp. 293–302.

[Hoefler et al., 2010]
T. Hoefler, C. Siebert, and A. Lumsdaine, “Scalable communication protocols
for dynamic sparse data exchange,” ACM Sigplan Notices, vol. 45,
no. 5, pp. 159–168, 2010.

[Harrison et al., 2015]
C. Harrison, J. Weiler, R. Bleile, K. Gaither, and H. Childs, “A
distributed-memory algorithm for connected components labeling of simulation
data,” in Topological and Statistical Methods for Complex Data.1em plus 0.5em minus 0.4emSpringer, 2015, pp. 3–19.

[Lu and Tryggvason, 2018]
J. Lu and G. Tryggvason, “Direct numerical simulations of multifluid flows in
a vertical channel undergoing topology changes,” Physical Review
Fluids, vol. 3, no. 8, p. 084401, 2018.

[Chiodi, 2020]
R. M. Chiodi, Advancement of numerical methods for simulating primary
atomization.1em plus 0.5em minus 0.4emCornell University,

[Jemison et al., 2015]
M. Jemison, M. Sussman, and M. Shashkov, “Filament capturing with the
multimaterial moment-of-fluid method,” Journal of Computational
Physics, vol. 285, pp. 149–172, 2015.

[Chirco et al., 2022]
L. Chirco, J. Maarek, S. Popinet, and S. Zaleski, “Manifold death: a volume of
fluid implementation of controlled topological changes in thin sheets by the
signature method,” Journal of Computational Physics, vol. 467, p.
111468, 2022.

[Tran et al., 2022]
H. D. Tran, M. Fernando, K. Saurabh, B. Ganapathysubramanian, R. M. Kirby, and
H. Sundar, “A scalable adaptive-matrix spmv for heterogeneous
architectures,” in 2022 IEEE International Parallel and Distributed
Processing Symposium (IPDPS).1em plus 0.5em minus 0.4emIEEE,
2022, pp. 13–24.

[Tomov et al., 2009]
S. Tomov, J. Dongarra, V. Volkov, and J. Demmel, “Magma library,” Univ.
of Tennessee and Univ. of California, Knoxville, TN, and Berkeley, CA, 2009.

[Heisler et al., 2022]
E. Heisler, A. Deshmukh, and H. Sundar, “: Domain specific language and code
generation for finite element and finite volume in julia,” in
International Conference on Computational Science.1em plus
0.5em minus 0.4emSpringer, 2022, pp. 118–132.

[Alnæs et al., 2015]
M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,
C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, “The fenics project
version 1.5,” Archive of Numerical Software, vol. 3, no. 100, 2015.