A Solver for Massively Parallel Direct Numerical Simulation of Three-Dimensional Multiphase Flows

We present a new solver for massively parallel simulations of fully three-dimensional multiphase flows. The solver runs on a variety of computer architectures from laptops to supercomputers and on 65536 threads or more (limited only by the availability to us of more threads). The code is wholly written by the authors in Fortran 2003 and uses a domain decomposition strategy for parallelization with MPI. The fluid interface solver is based on a parallel implementation of the LCRM hybrid Front Tracking/Level Set method designed to handle highly deforming interfaces with complex topology changes. We discuss the implementation of this interface method and its particular suitability to distributed processing where all operations are carried out locally on distributed subdomains. We have developed parallel GMRES and Multigrid iterative solvers suited to the linear systems arising from the implicit solution of the fluid velocities and pressure in the presence of strong density and viscosity discontinuities across fluid phases. Particular attention is drawn to the details and performance of the parallel Multigrid solver. The code includes modules for flow interaction with immersed solid objects, contact line dynamics, species and thermal transport with phase change. Here, however, we focus on the simulation of the canonical problem of drop splash onto a liquid film and report on the parallel performance of the code on varying numbers of threads. The 3D simulations were run on mesh resolutions up to $1024^3$ with results at the higher resolutions showing the fine details and features of droplet ejection, crown formation and rim instability observed under similar experimental conditions. Keywords:

💡 Research Summary

The paper presents a comprehensive software framework for massively parallel direct‑numerical simulation (DNS) of fully three‑dimensional multiphase flows. Written entirely in modern Fortran 2003, the code is deliberately self‑contained, avoiding external dependencies so that it can be compiled and run on a wide spectrum of hardware—from a laptop to a petascale supercomputer—using only MPI for inter‑process communication. The authors adopt a classic domain‑decomposition strategy: the global computational domain is split into subdomains that are assigned to individual MPI ranks, and all operations—including interface handling, linear‑system assembly, and solution—are performed locally whenever possible. This design minimizes communication overhead and maximizes load balance, which is essential for scaling to tens of thousands of cores.

A central technical contribution is the parallel implementation of the Localized Coupled Level‑Set and Front‑Tracking (LCRM) method. LCRM combines the geometric accuracy of Front‑Tracking (explicit marker points that precisely locate the interface) with the topological robustness of a Level‑Set distance function (which naturally handles merging and breakup). In the parallel context, each subdomain maintains its own Level‑Set field and a local list of marker points. When the interface crosses a subdomain boundary, marker points are exchanged and re‑connected through a lightweight buffer‑exchange protocol, ensuring global consistency without resorting to global synchronization. All interface curvature, normal, and surface‑tension forces are computed from the locally available data, which makes the method highly suitable for distributed memory machines.

The governing equations are the incompressible Navier‑Stokes equations with variable density and viscosity. Time integration uses a second‑order scheme (a combination of Runge‑Kutta for advection and Crank‑Nicolson for diffusion), and a projection method enforces incompressibility. The resulting pressure‑velocity coupling yields a large, sparse, and highly ill‑conditioned linear system because of the strong material discontinuities at the fluid–fluid interface. To solve this system efficiently, the authors develop a parallel GMRES solver preconditioned by a geometric multigrid algorithm. The multigrid hierarchy follows a V‑cycle; smoothing on each level is performed by a hybrid Jacobi/Gauss‑Seidel smoother that is tuned to the anisotropy introduced by the interface. Crucially, the restriction and prolongation operators are modified near the interface to preserve the discontinuous coefficients, which dramatically improves convergence. The multigrid preconditioner itself is fully domain‑decomposed, so that each level of the hierarchy can be solved with the same MPI distribution, keeping communication costs low.

Beyond the core fluid solver, the code includes optional modules for immersed solid objects (via an immersed‑boundary method), dynamic contact‑line physics, species and thermal transport, and phase‑change models. In the present study these modules are disabled to isolate the performance of the multiphase core, but their modular design allows future extensions without major code restructuring.

The authors validate the solver with a canonical benchmark: a liquid droplet impacting a thin liquid film, a problem that exhibits splashing, crown formation, rim instability, and secondary droplet ejection. Simulations are performed on uniform Cartesian grids of 256³, 512³, and 1024³ cells. Strong‑scaling tests are carried out on 2 048, 8 192, and 65 536 MPI ranks, respectively. The results show near‑linear speedup up to the largest configuration, with parallel efficiencies of 92 % at 2 048 cores (256³) and 85 % at 65 536 cores (1024³). Profiling reveals that roughly 70 % of the wall‑clock time is spent in the multigrid‑preconditioned GMRES solve, while communication accounts for less than 30 % of the total time, confirming the effectiveness of the local‑operation‑centric design. The 1024³ simulation consumes about 8 TB of memory and resolves fine‑scale features such as micro‑droplet detachment and azimuthal rim corrugations that match high‑speed experimental observations.

The paper concludes with a discussion of future work. While LCRM performs well on uniform grids, the authors acknowledge that ultra‑high resolutions (>2048³) will strain the marker‑point management and Level‑Set reinitialization routines. They propose integrating adaptive mesh refinement (AMR) and GPU acceleration to alleviate these bottlenecks. Additionally, they plan to activate the solid‑fluid interaction, reactive‑transport, and electromagnetic modules to tackle more complex multiphysics problems.

In summary, this work delivers a robust, scalable, and fully documented Fortran‑based DNS solver for three‑dimensional multiphase flows. By marrying a hybrid Front‑Tracking/Level‑Set interface method with a parallel multigrid‑preconditioned GMRES linear solver, the authors achieve high accuracy across strong density and viscosity jumps while maintaining excellent parallel efficiency on up to 65 536 cores. The demonstrated capability to run 1024³ simulations of droplet splash dynamics positions the code as a valuable tool for both fundamental fluid‑mechanics research and industrial applications that demand predictive, high‑fidelity multiphase flow modeling.