A Simple Parallel Solution Method for the Navier–Stokes Cahn–Hilliard Equations

: We present a discretization method of the Navier–Stokes Cahn–Hilliard equations which offers an impressing simplicity, making it easy to implement a scalable parallel code from scratch. The method is based on a special pressure projection scheme with incomplete pressure iterations. The resulting scheme admits solution by an explicit Euler method. Hence, all unknowns decouple, which enables a very simple implementation. This goes along with the opportunity of a straightforward parallelization, for example, by few lines of Open Multi-Processing (OpenMP) or Message Passing Interface (MPI) routines. Using a standard benchmark case of a rising bubble, we show that the method provides accurate results and good parallel scalability.


Introduction
Phase-field (or diffuse-interface) models have become a valuable tool for the numerical simulation of two-phase flows. In many fields these methods have been proven advantageous compared to traditional sharp interface methods (e.g., Level-Set, Arbitrary Lagrangian-Eulerian, Volume-Of-Fluid). Phase field methods can easily handle (i) moving contact lines (e.g., [1]), (ii) topological transitions (e.g., [2,3]), and (iii) additional chemical species in the fluid and on the interface (e.g., [4]). Further, these methods are relatively robust and do not require additional stabilization techniques such as reinitialization or convection stabilization. It has been shown in recent years that phase field models can be computationally cheap [5] and accurate [6,7] depending on the simulated problem.
Mathematically, phase field methods for two-phase flow are typically described by the Navier-Stokes Cahn-Hilliard equations. Efficient techniques have been proposed for their discretization in space (e.g., [8]) and in time (e.g., [5,9,10]). Due to the necessary high resolution at the fluid-fluid interface, typical discretizations are based on adaptive finite element methods which comes along with a complexity that makes implementations from scratch almost impossible. In addition, the scalable parallelization of such methods is typically challenging and requires a solid knowledge on memory-distributed matrices and vectors and parallel preconditioning and solution techniques. However, given the current shift from serial towards parallel computing, parallelizability is the key component of future numerical algorithms. The goal of this paper is to combine parallelizability with simplicity in a single numerical method. Therefore, we propose a parallel finite difference scheme for the Navier-Stokes Cahn-Hilliard equations.
A finite difference scheme is arguably the simplest way to discretize partial differential equations in space. Efficient finite difference schemes have been proposed for phase field methods, e.g., in [2,3]. Recently, the method has been extended to feature mass conservation and energy stability [11]. All these approaches are rather difficult to understand for a novice and complicated to parallelize as they all involve solvers for coupled systems of equations. Parallel implementations from scratch usually require several thousand lines of code.
Here we aim to avoid the solution of (non-)linear systems of equations by employment of an explicit (forward) Euler method. In this case all unknowns decouple, no system of equations needs to be assembled, and all degrees of freedom can be directly computed from values of the previous time step. This enables straightforward parallelization by few lines of Open Multi-Processing (OpenMP) or Message Passing Interface (MPI) routines. We make this explicit time discretization applicable to these equations of saddle-point structure by a special pressure projection scheme with incomplete pressure iterations. It is known that explicit Euler methods introduce strong time step restrictions of the form dt < Ch 2 for a second order equation and dt < Ch 4 for a fourth order equation, such as the Cahn-Hilliard equation. However, we will show that in practice these restrictions are much less severe and that our proposed method is very efficient as the number of computations per time step is low.
Accordingly, we develop a parallel solution method for the Navier-Stokes equations which can be understood and implemented even for a novice within a few days. We share the (short) code attached to this publication and expect it to be valuable for PhD students and numerical analysts who want to have an easy-to-understand code basis at hand to solve two-phase flow problems with a fast parallel method. Using a standard benchmark case of a rising bubble, we show that our method can provide accurate results in a highly efficient, yet very simple, parallel implementation.
The structure of the remaining paper is as follows. In Section 2 we present the Navier-Stokes Cahn-Hilliard equations. In Sections 3 and 4, we present the time and space discretization, respectively. The resulting method is challenged in a standard two-phase flow benchmark in Section 5. Afterwards we give some comments on the attached implementation of the method in Section 6. Finally, conclusions are drawn in Section 7.

Navier-Stokes Cahn-Hilliard Equations for Two-Phase Flow
The phase field method was originally developed to model solid-liquid phase transitions, see, e.g., [12][13][14]. The interface is represented as a thin layer of finite thickness ε and an auxiliary function, the phase field φ, is used to indicate the phases. The phase field function varies smoothly between distinct values in both phases and the interface can be associated with an intermediate level set of the phase field function. Diffuse interface approaches for mixtures of two immiscible, incompressible fluids lead to the Navier-Stokes Cahn-Hilliard equations and have been considered by several authors, see, e.g., [8,[15][16][17]. If variable (phase-dependent) densities are involved, different formulations of the model equations have been proposed [17][18][19][20]. However, in real scenarios away from the critical point, the involved interface thickness ε is small such that all these models yield comparable results [6]. Therefore, we use the simplest model which reads [17] in the domain Ω. Here u, p, φ, and µ are the velocity, pressure, phase field variable, and chemical potential, respectively. The parameter m is a mobility constant, ε defines a length scale over which the interface is smeared out. Furthermore, D(u) = ∇u + ∇u T is the (doubled) viscous rate of strain tensor, ρ(φ), η(φ), and f are the (phase dependent) density, viscosity, and body force. The polynomial φ 3 − φ arises as the first derivative of the double well potential W(φ) = 1/4(φ 2 − 1) 2 with respect to φ. This choice of the polynomial leads to the values φ ≈ ±1 in two distinct fluid phases. The parameter σ is a scaled surface tension which is a multiple of the physical surface tension σ. For the given double-well potential, the relationσ = σ 3 2 √ 2 holds [21].
Since the solution of the Cahn-Hilliard equation does not satisfy a maximum principle, we follow the approach in [22] and use a truncated phase field function φ = φ, |φ| ≤ 1, sign(φ), |φ| > 1, (5) to construct linear interpolations of the viscosity η and density ρ, where η ±1 and ρ ±1 denote the viscosity and density in the fluid phases where the phase field is ±1, respectively. Similar truncated phase field functions have played an essential role in the simulation of multi-agent systems [23].

Time Discretization
In this section we adopt the explicit Euler scheme for the time discretization of the Navier-Stokes Cahn-Hilliard equations. Let the time interval [0, T] be divided in N subintervals of size dt, such that t n := n · dt, where the upper index denotes the time step number.
Before applying the time discretization to the individual systems we use an operator splitting approach to decouple the Navier-Stokes equations from the Cahn-Hilliard equations. To be more precise, both problems are solved subsequently in every time step. First, we solve the Navier-Stokes system using the values of φ and µ from the previous time step. Then we solve the Cahn-Hilliard system using the just computed velocity field.

Time Discretization of the Navier-Stokes Equations
The main difficulty in numerically solving the Navier-Stokes equations is that the pressure field is coupled indirectly to the velocity by the incompressibility condition. Chorin [24] and Temam [25] proposed an effective way to decouple the pressure from the velocity by a projection type scheme. As their original scheme has poor accuracy in time many improved versions of these schemes have been introduced since, for details we refer to the review article [26]. Keeping in mind that we aim for an explicit time discretization which enforces small time steps, the time stepping error of the original scheme becomes less important.
In terms of simplicity, the original scheme of Chorin [24] and Temam [25] surpassed all its succeeding versions. The idea is to decouple the momentum and mass balance equations. To this end, an intermediate velocity fieldũ n+1 is computed, which satisfies the momentum balance (without the pressure term), but is not divergence-free. Then this field is projected to the space of divergence free functions to obtain the velocity u n+1 . The latter step involves to first solve an additional equation for the pressure field p n+1 which is then used to update the velocity field. Combined with an explicit time discretization the scheme reads: where ∂/∂n denotes the normal derivative on the domain boundary ∂Ω.
Equations (8)-(10) are decoupled and can be solved subsequently. The first equation is used to compute an intermediate (not divergence-free) velocity fieldũ n+1 . The second equation computes the pressure which is then used to obtain a divergence-free velocity field u n+1 from the last equation. Hence, Equations (9) and (10) can be seen as the projection ofũ n+1 onto the space of divergence-free vector fields.
The first and third equation are completely explicit, i.e., the sought variablesũ n+1 and u n+1 can be directly computed. The second equation is a Poisson equation and is usually handled by assembling a matrix which can be solved with the conjugate gradient method or by shifting the problem to Fourier space and solving with spectral methods.
Again, as we aim for simplicity we follow a different approach here. The idea is to make Equation (9) amenable to an explicit Euler method by including a time derivative of the pressure. In fact, the solution p n+1 from Equation (9) can be seen as the stationary solution of the convergent heat equation Instead of solving this equation exactly, we propose to solve only K steps of the time evolution by an explicit Euler scheme. Accordingly, we introduce K − 1 intermediate values p n+k/K for k = 1, . . . , K−1 which can be computed successively from with a different time step sizedt. As we aim for the stationary solution it is desired to choosedt as large as possible. From the explicit Euler scheme, the time step restrictiondt < ρh 2 /4 emerges. To be well below this barrier, we use half of this time step size,dt = ρh 2 /8, throughout this work.
As the time step size of the overall projection method, dt, is small, we can expect that the velocity and pressure change only slightly in every time step. Hence, very few iterations K of Equation (11) should suffice to update the pressure. In the simplest case, K = 1, Equation (11) reduces to a direct computation of p n+1 by This choice can be sufficient if the time step size is small enough.

Time Discretization of the Cahn-Hilliard Equations
The explicit Euler scheme for the Cahn-Hilliard equation reads Due to the explicit presence of a fourth order space derivative, a time step restriction arises: dt < Ch 4 /(mε 2 ). Luckily, this restriction turns out to be less severe after the following considerations. Firstly, the interface length ε is very small. The real thickness of a fluid-fluid interface is in the order of nanometers, unless the temperature is close to the critical point. This length scale cannot be resolved in numerical methods as soon as the length scale of the problem exceeds the micrometer range. Hence, in any useful numerical approximation that aims to be close to the real physical solution, ε is chosen as small as possible. In practice, that means ε is chosen similar to the grid size, ε ≈ h. Consequently, the above time step restriction effectively reduces to dt < Ch 2 /m. This condition is weakened further, as m is typically very small. The interface moves physically with the fluid velocity, thus m needs to be so small that the right-hand side of Equation (12) hardly moves the interface but only conserves the interface shape, similar to reinitialization in Level-Set methods. In addition, matched asymptotic analysis shows the convergence of diffuse-interface methods toward the sharp interface equations only for small CH mobility [20].
The remaining time step restriction can be tolerated, given that the number of computations in each time step is relatively low in comparison to implicit time discretization methods. While the latter reduces time step restrictions, they typically involve several orders more computations per time step, due to the necessary solution of a (linear) system of equations, which is avoided in our method.

Space Discretization
The time discrete scheme presented in the previous section is discretized in space by a finite difference method. We use an equidistant rectangular grid with equal spacing h in all directions. For simplicity the following discussion is restricted to two dimensions, but we note that the scheme can be easily extended to three dimensions. The rectangular computational domain Scalar variables are defined in the centers of these cells. One additional layer of ghost points is included outside of each boundary facet. To stay close to the implementation we use only integer indices to describe the values of all variables on the grid points. Hence, the discrete grid values of a scalar field q(x, y) are defined as where q is a placeholder for φ, µ, p.
To avoid spatial oscillations from the odd-even decoupling between pressure and velocity, a staggered Cartesian grid is employed [27]. On a staggered grid the velocity variables (u 0 , u 1 ) =: u are shifted by h/2 in their respective direction to the cell faces. Figure 1 shows the arrangement of all the different variables on the grid. To stay close to the implementation we use only integer indices to describe the values of all variables on the grid points: Note, that ghost points are also included for the velocity field. While this might not be necessary, depending on the boundary condition, it simplifies and unifies indexing and implementation at negligible additional memory consumption.
We introduce for some quantity q the discrete Laplacian ∆ h , all-sided average A, horizontal and vertical averages A x , A y , backward differences D x , D y , and central differences D 2x , D 2y as follows:

Matched Densities and Viscosities
In the beginning let us consider the simple case of matched viscosity (η = η −1 = η +1 ) and matched density (ρ = ρ −1 = ρ +1 ). In this case the viscous stress contribution reduces to η∆u and the space discretization of Equation (8) reads u n+1 Boundary conditions to the velocity are included in this step as explained in the Appendix A. After computingũ n+1 directly from these equations, the pressure p n+1 is calculated from the discretized version of Equation (11), which is solved for k = 0, . . . K−1 to successively compute p n+ 1 K i,j , p n+ 2 K i,j , ..., p n+1 i,j , for a small constant K.

Finally, the new velocity field is computed from
After calculation of the new velocity, the Cahn-Hilliard system is solved. Using central differences for the advective term, a second order in space discretization of the Cahn-Hilliard Equations (12) and (13) reads The homogeneous Neumann boundary conditions are included as explained in the Appendix A.
Equations (14)- (19) provide a way to directly solve the Navier-Stokes Cahn-Hilliard system. We again note the equations can be solved individually for each i and j, no coupled system of equations is involved. Hence, the solution algorithm can be easily implemented with a few dozen lines of code. Due to the direct calculation, parallelization (e.g., with OpenMP or MPI) is straightforward, can be done with a few lines of code and can be expected to scale very well.

Non-Matched Density and Viscosity
In the general case of non-matched density and viscosity the discretization becomes more elaborate. We define η n i,j = η(φ n i,j ) and ρ n i,j = ρ(φ n i,j ) and replace ρ in Equations (14) and (17) by A x ρ n i,j and in Equations (15) and (18) by A y ρ n i,j . Additionally, in Equation (16) we replace ρ by ρ n i,j and ∆ h p To deal with non-constant viscosity we introduce the fields α, β, γ to represent the viscous stress, In the discrete case, α and β are naturally defined in the cell centers, γ is defined on the cell vertices. This leads to α n i,j = 2η n i,j D x u n 0 i+1,j , β n i,j = 2η n i,j D y u n 1 i,j+1 , γ n i,j = D y u n 0 i,j + D x u n 1 i,j Aη n i,j .
These terms are then used to replace the viscous terms η∆ h u n 0 i,j and η∆ h u n 1 i,j in Equations (14) and (15) by D x α n i,j + D y γ n i,j+1 and D y β n i,j + D x γ n i+1,j , respectively.

Benchmark Test Results
To validate the accuracy of the model, a benchmark test as described by test case 1 in [28] is conducted. This test can be seen as the standard benchmark for two-phase flows and has been conducted for many different numerical methods (e.g., Level-Set, Volume-of-Fluid, Arbitrary Lagrangian-Eulerian). The benchmark was applied to various phase field methods in [6].

Test Scenario
The test specifies a circular bubble of one fluid (here φ ≈ 1), which rises in a rectangular domain Ω = [1, 0] × [2, 0] filled with another fluid (here φ ≈ −1). The initial position of the center of the bubble is at (0.5, 0.5) and its radius is 0.25, which results in the initial condition Test case 1 of [28] prescribes for the outer fluid the physical quantities ρ −1 = 1000 and η −1 = 10. The inner fluid of the bubble is defined by ρ +1 = 100 and η +1 = 1. Furthermore, the gravitational force vector is set to f = (0, −0.98(ρ(φ) − ρ −1 )) and the surface tension σ to 24.5. No slip boundary conditions are specified at the top and bottom, a free slip condition is imposed at the lateral boundaries. Due to the different densities of the inner and outer fluids, the bubble will rise in the vertical tube over time. The simulation is run for 3 time units during which some benchmark quantities are computed.
The chosen four combinations of numerical parameters are listed in Table 1. Their choice is guided by reported values in the phase field benchmark [6]: The interface thickness ε is refined proportional to grid size h, which keeps the number of elements along the interface fixed. The mobility is set to a small constant, m = 0.01, which is by a factor ten higher than the value chosen in [6] (acknowledging the different scaling of mu there). We conducted a preliminary study to determine that this value of m is large enough to keep a good shape of the phase field interface profile while introducing only small (undesired) Cahn-Hilliard diffusion, which is confirmed by the accuracy of the results in Section 5.3. In addition, note that results are essentially independent of m when → 0, see [20]. Given the previously discussed time step restrictions, the time step size is scaled quadratically with grid size.
Similar to [6], the total number of degrees of freedom per time step is reported for the four different simulations. Additionally, the CPU time of each parameter case is listed in Table 1 for an OpenMP implementation running on four cores.

Measured Quantities
In order to compare our results with the literature some of the benchmark quantities presented in [28] are used. Their calculation is described below: (1) Center of Mass Rise Velocity . (

3) Error Norms and Orders of Convergence
The following three error norms are also specified in [28] and will be used for the evaluation of our results: where q n is the temporal evolution of quantity q. We use the results of the simulation with the finest grid and smallest time steps as the reference solution q n re f . Spline interpolation is used to match the time steps of the reference solution with our solution. Using the so defined errors, a rate of convergence (ROC) can be computed as follows: where error h denotes the error as defined above on a grid of grid size h.

Bubble Shape
As already described in [28], the initial circular bubble changes its shape over time to an ellipsoidal shape. In Figure 2 one can see the different bubble shapes at the end time t = 3 compared to the reference results from group 1 in [28]. It shows that the shapes of our model approximate the reference solution as grid resolution and interface thickness are refined. Especially the final bubble shapes of the two finest grids with h = 1/128 and h = 1/256 got very close to the reference results, as they are wider than the bubble shapes produced by the other two simulations.

Benchmark Quantities
In addition to the bubble shapes, we evaluated different benchmark quantities as defined before. Figure 3 shows the center of mass and the rise velocity of the bubble over time in comparison to the reference solution from [28]. Both measured quantities approach the reference values with increasing grid resolutions. In the "eye norm" the center of mass is nearly identical for the two finest grids, while the simulation with h = 1/32 is more off. Differences between the parameter cases are more pronounced when looking at the rise velocity. Here, one can see more clearly that all curves approach the reference solution with increasing grid resolutions. This behavior is also confirmed in Table 2 which lists the maximum rise velocity and the final center of mass. A comparison with the phase field results from [6] reveals similarly small deviations from the reference solutions in [28] which indicates that these differences stem mostly from the modeling error due to non-zero interface thickness. Table 2. Benchmark quantities: Maximum rise velocity V c,max with corresponding time t| V c =V c,max and final position of the center of mass y c (t = 3) with reference results from group 1 in [28].

Error Norms
Finally, the calculated error norms and the resulting orders of convergence are presented in Table 3. Overall, the previous observations are confirmed, since all the calculated errors become smaller with finer grid resolutions and shorter time step lengths. The ROC values indicate clear second order convergence for all error norms. Note that this implies second order convergence in space and first order in time as time step size is scaled quadratically with the grid size.

Parallelization
Finally, we tested the model for different numbers of cores to examine the parallel efficiency. The resulting time measurements and calculated speedup can be found in Figure 4. Here, the dimensions of the field were 512 × 1024 grid points and the number of time steps was set to 3000. It can be seen that the implementation of the proposed discretization method is highly parallelizable leading to a nearly perfect linear speedup. It is noteworthy that parallelization has been achieved by a few simple OpenMP 'parallel for' directives. Given the problem size is large enough, a similarly good linear speedup can also be expected for much larger number of cores, for OpenMP as well as MPI parallelization.

Implementation
We share the code in the supporting material attached to this publication and additionally on GitHub: https://github.com/soranad-a/parallel-navier-stokes-cahn-hilliard. It contains the discretized Navier-Stokes Cahn-Hilliard system implemented in C++, whereby the variable names were kept analogous to those in the equations above. Furthermore, the code includes OpenMP for parallelization on multiple cores. It consists of the following three files: • helpers.h-handles output and data structure of 2 dimensional arrays. • config.conf-self explanatory file defining all constants and simulation parameters. • main.cc-contains core code of the numerical algorithm, but most of the lines are needed for reading parameters and managing output.
The code can be compiled by: g++ -std=c++17 -fopenmp -O3 main.cc. To enable the correct output of the simulation files, the program needs a subfolder named paraview/into which it writes its output files, which can be visualized in the open-source data analysis software ParaView (www.paraview.org). Furthermore, during runtime, a file named benchmark.csv is created, which contains the measured benchmark quantities.

Conclusions
We have presented a discretization method of the Navier-Stokes Cahn-Hilliard equations. The method offers an impressing simplicity, which is based on an explicit Euler time discretization of a special pressure projection scheme with incomplete pressure iterations. Hence, it is not necessary to solve a linear system of equations as all unknowns decouple, and computations can be done separately for each grid point. Consequently, the method can be easily implemented from scratch and parallelization is straightforward by few lines of OpenMP or MPI routines.
A standard benchmark case of a rising bubble is used to evaluate the method in practice. We find that the arising time step restrictions are not critical as the number of compute instructions per time step is relatively low. We show that the method provides accurate solutions comparable to previous results from finite-element methods [6,28]. The numerical results approach reference values as grid size and interface thickness are decreased. The computed ROC indicates second order convergence. Finally, we measure computation time for up to 24 CPU cores and find a nearly perfect parallel speedup.
We share the (short) code attached to this publication and expect it to be valuable for novices in phase field simulations who want to have an easy-to-understand code basis at hand to solve two-phase flow problems with a fast parallel method. The method is flexible and can easily be augmented to account for additional physical processes. In the future we plan to extend the method to describe flows in complex evolving geometries [29] and fluid-structure interaction using the approach in [30]. We further address the large-scale parallelization by use of GPUs.