Werner Verdier
email: werner.verdier@cea.fr

Pierre Kestener
email: pierre.kestener@cea.fr

Alain Cartalade
email: alain.cartalade@cea.fr

Performance portability of lattice Boltzmann methods for two-phase flows with phase change

Keywords: Lattice Boltzmann method, phase-field model, two-phase flows with phase change, performance portability, Kokkos library, LBM saclay, conservative Allen-Cahn model

Numerical codes using the lattice Boltzmann methods (LBM) for simulating one-or two-phase flows are widely compiled and run on graphical process units. However, those computational units necessitate to re-write the program by using a lowlevel language which is suited to those architectures (e.g. CUDA for GPU NVIDIA or OpenCL). In this paper we focus our effort on the performance portability of LBM i.e. the possibility of writing LB algorithms with a high-level of abstraction while remaining efficient on a wide range of architectures such as multicores x86, GPU NVIDIA , ARM, and so on. For such a purpose, implementation of LBM is carried out by developing a unique code, LBM saclay written in the C++ language, coupled with the Kokkos library for performance portability in the context of High Performance Computing. In this paper, the LBM is used to simulate a phase-field model for two-phase flow problems with phase change. The mathematical model is composed of the incompressible Navier-Stokes equations coupled with the conservative Allen-Cahn model. Initially developed in the literature for immiscible binary fluids, the model is extended here to simulate phase change occurring at the interface between liquid and gas. For that purpose, a heat equation is added with a source term involving the time derivative of the phase field. In the phase-field equation a source term is added to approximate the mass production rate at the interface. Several validations are carried out to check step-by-step the implementation of the full model. Finally, computational times are compared on CPU and GPU platforms for the physical problem of film boiling.

Introduction

The Lattice Boltzmann Method (LBM) [START_REF] Krueger | The Lattice Boltzmann Method: Principles and Practice[END_REF][START_REF] Guo | Lattice Boltzmann Method and Its Applications in Engineering[END_REF] is a very attractive method to simulate problems involving fluid flows. Since more than ten years, numerical codes using that method are widely compiled and run on Graphical Process Units (GPU) [START_REF] Li | Implementing lattice Boltzmann computation on graphics hardware[END_REF][START_REF] Tölke | Implementation of a Lattice Boltzmann kernel using the Compute Unified Device Architecture developed by nVIDIA[END_REF][START_REF] Kuznik | LBM based flow simulation using GPU computing processor[END_REF][START_REF] Zhou | GPU implementation of lattice Boltzmann method for flows with curved boundaries[END_REF]. The GPUs allow for a very high calculation throughput and they are particularly efficient for repetitive workloads with simple memory access patterns. These units were initially designed for image processing or graphics rendering, but LBM simulations can also benefit from their use, because the stages of streaming and collide are two simple (stencil-like) computational operations. Numerous works have demonstrated the efficiency of LBM on single GPU (e.g. [START_REF] Kuznik | LBM based flow simulation using GPU computing processor[END_REF]) and later on clusters of GPUs (e.g. [START_REF] Obrecht | Multi-GPU implementation of the lattice Boltzmann method[END_REF][START_REF] Januszewski | Sailfish: A flexible multi-gpu implementation of the lattice boltzmann method[END_REF]). However, those computational units necessitate to re-write the code by using a low-level language which is suited to their specific architectures (e.g. CUDA for GPU NVIDIA or OpenCL). In this paper we focus our effort on the performance portability of LBM i.e. the possibility of writing LBM algorithms with a high-level of abstraction, but by remaining efficient on a wide range of architectures such as multicores x86, GPU NVIDIA , ARM, and so on.

The issue of performance portability has already been studied and implementation of numerical algorithms running on various architectures (GPU and so on) can be done by directive approaches (mostly OpenMP or OpenACC). Directivebased parallel programming solutions consist in decorating source code with comments that are interpreted by the compiler to derive the actual parallel code. They are useful when porting a legacy simulation code with a large number of lines, for which it is not reasonable to rewrite it from scratch. However, those programming models deal with computational patterns (for loops, reduction loops, ...) and do not provide tools for data or memory containers. Here we present an application of a more promising approach that uses a library-based solution which offers high-level abstract programming concepts and hardware agnostic solution for a better integration into C++ codes. Among libraries sharing the same goal of performance portability (like RAJA or SYCL), the Kokkos library [START_REF] Edwards | Kokkos: Enabling manycore performance portability through polymorphic memory access patterns[END_REF] is used for simulating two-phase flows with LBM. Kokkos implements a programming model in C++ for writing performance portable applications targeting all major High Performance Computing (HPC) platforms. Programming tools provide abstractions for both parallel execution of code and data management, i.e. they provide memory containers (multidimensional arrays) where the actual memory layout will be chosen by the library during compilation. Directive-based solution does not provide such advanced features regarding memory. The Kokkos library can currently use OpenMP, Pthreads and CUDA as backend programming models. The library has already been applied to accelerate high-order mesh optimization in [START_REF] Eichstädt | Accelerating high-order mesh optimisation with an architecture-independent programming model[END_REF].

Because of its explicit scheme and local interactions, the LBM ideally exploits the massively parallel supercomputers based on either CPUs or GPUs or heterogeneous architectures. In this paper, we take advantage of those benefits to study two-phase flows. Several topical reviews exist in the literature for modeling two-phase flows in LBM framework [START_REF] Huang | Multiphase Lattice Boltzmann Methods. Theory and Application[END_REF][START_REF] Li | Lattice Boltzmann methods for multiphase flow and phase-change heat transfer[END_REF]. The main families of methods are the color-gradient method [START_REF] Gunstensen | Lattice Boltzmann model of immiscible fluids[END_REF], the pseudo-potential method [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF], the free-energy method [START_REF] Swift | Lattice boltzmann simulations of liquid-gas and binary fluid systems[END_REF], and the phase-field method [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]. Most of approaches consider the interface as a diffuse zone (characterized by a thickness and a surface tension) which can be seen as a small region of transition between bulk phases. In pseudo-potential methods [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF][START_REF] Shan | Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation[END_REF] an additional force term is added in the Navier-Stokes equations to take into account an equation of state which is not the classical law of perfect gases [START_REF] Yuan | Equations of state in a lattice Boltzmann model[END_REF]. In that case, the density plays the role of a phase index varying smoothly between densities of gas and liquid. Several recent applications use that method for simulating liquid-gas phase change [START_REF] Li | Lattice Boltzmann simulations for transition from dropwise to filmwise condensation on hydrophobic surfaces with hydrophilic spots[END_REF][START_REF] Li | Lattice Boltzmann simulation of condensation in the presence of noncondensable gas[END_REF]. Another class of diffuse interface methods is the color-gradient model [START_REF] Gunstensen | Lattice Boltzmann model of immiscible fluids[END_REF] for which two distribution functions are introduced for computation of each phase (red and blue). In those approaches, surface tension is derived from a recoloring step involving both distribution functions [START_REF] Leclaire | Numerical evaluation of two recoloring operators for an immiscible two-phase flow lattice boltzmann model[END_REF][START_REF] Leclaire | Enhanced equilibrium distribution functions for simulating immiscible multiphase flows with variable density ratios in a class of lattice boltzmann models[END_REF]. The final approach that is commonly applied in the LBM literature is to capture the interfacial behavior through a phase-field equation. In this paper, we follow this latter method: the phase-field theory for two-phase flows [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF]. The phase-field method is quite similar to the free-energy lattice Boltzmann method [START_REF] Swift | Lattice boltzmann simulations of liquid-gas and binary fluid systems[END_REF] in the sense that both models are thermodynamically consistent and can be derived from a free-energy functional. However, in the free-energy LB approach, the density gradient appears explicitly in the free-energy functional and the phase separation is described by a non-ideal equation of state. For that purpose, the equilibrium distribution function is modified to include a non-ideal thermodynamic pressure tensor. In this paper, both fluids are considered as quasi-incompressible, i.e. we assume that the incompressibility condition holds in the bulk phases except in the interfacial zone where the mass production rate ṁ acts. That mass production rate comes from the phase change between the gas and the liquid. A new function φ is introduced to track the interfacial zone where the density varies.

Two main phase-field models for interface tracking between two immiscible fluids exist in literature: the first one is the Cahn-Hilliard (CH) model [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF][START_REF] Cahn | Free Energy of a Nonuniform System. I. Interfacial Free Energy[END_REF][START_REF] Jasnow | Coarse-grained description of thermo-capillary flow[END_REF] which was extensively applied in LBM literature for simulating spinodal decomposition [START_REF] Kendon | Inertial effects in three-dimensional spinodal decomposition of a symmetric binary fluid mixture: a lattice Boltzmann study[END_REF], buoyancy of bubbles [START_REF] Zheng | A lattice boltzmann model for multiphase flows with large density ratio[END_REF], drop impact [START_REF] Lee | Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces[END_REF], Rayleigh-Taylor instability [START_REF] Zu | Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts[END_REF] and so on. The second one is a more recent model, called the conservative Allen-Cahn (CAC) model, which was first developed in [START_REF] Sun | Sharp interface tracking using the phase-field equation[END_REF] and derived in conservative form in [START_REF] Chiu | A conservative phase field method for solving incompressible two-phase flows[END_REF]. The model became popular in the LB community [START_REF] Geier | Conservative phase-field lattice Boltzmann model for interface tracking equation[END_REF][START_REF] Fakhari | A weighted multiple-relaxation-time lattice Boltzmann method for multiphase flows and its application to partial coalescence cascades[END_REF][START_REF] Mitchell | Development of a three-dimensional phase-field lattice Boltzmann method for the study of immiscible fluids at high density ratios[END_REF] and several papers compare the Cahn-Hilliard and conservative Allen-Cahn models, e.g. [START_REF] Lee | Comparison study of the conservative Allen-Cahn and the Cahn-Hilliard equations[END_REF] without LBM and [START_REF] Wang | Comparative study of the lattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations[END_REF] with LBM. In this work the CAC model is chosen for interface tracking in order to eliminate the curvature-driven interface motion which is implicitly contained in the CH equation (see Section 2). Moreover, the CAC model involves only a second-order derivative and does not require to compute the fourth-order derivative (Laplacian of chemical potential) which appears in the CH equation.

In this paper, we take advantage of the simplicity of LBM to develop a new portable code for simulating two phase flows with the coupled Navier-Stokes/conservative Allen-Cahn (NS/CAC) model. The new code, called LBM saclay, targets all major HPC platforms such as multi-GPUs and multi-CPUs. In this paper, we also check the capability of the NS/CAC model to simulate phase change problems in the vicinity of the critical temperature. Near the critical temperature, properties of each phase vary smoothly and the range of variation of those parameters remains small. Several fluid flow models of phase change have already been proposed in the literature with the Cahn-Hilliard equation [START_REF] Dong | Lattice Boltzmann Simulation of Growth and Deformation for a Rising Vapor Bubble Through Superheated Liquid[END_REF][START_REF] Safari | Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow[END_REF]. Following those references, the NS/CAC model is extended here by adding a source term in both the mass balance and the CAC equations. The source term involves the mass production rate ṁ occurring at the interface. In references [START_REF] Dong | Lattice Boltzmann Simulation of Growth and Deformation for a Rising Vapor Bubble Through Superheated Liquid[END_REF][START_REF] Safari | Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow[END_REF], the liquid is often considered at saturation temperature and its thermal conductivity is neglected. Under those assumptions, ṁ is calculated by a gradient operator (Fourier's law) involving only the thermal conductivity of gas. Moreover, in order to avoid computing the temperature equation in liquid phase (because the thermal conductivity is neglected), a cut-off value of the phase-field is introduced beyond which the temperature equation will not be computed [START_REF] Safari | Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow[END_REF]. Here we propose an alternative way to calculate ṁ that avoids computing this gradient and avoids introducing this cut-off value. For that purpose, ṁ will be related to the normal interface velocity and expressed as a source term close to what is done in solidification models (section 2.3). Implementation of lattice Boltzmann methods will be checked step-by-step by considering separately solutions of the phase-field equation, the phase-field coupled successively with a fluid flow, and the phase-field coupled only with temperature for which the ratio of physical properties remain low. Finally, the aspects of two-phase flow, phase change and heat transfer are coupled to simulate the phenomena of film boiling [START_REF] Kharangate | Review of computational studies on boiling and condensation[END_REF].

This paper is organized as follows. Section 2 presents the continuous mathematical model based on the conservative Allen-Cahn equation which is extended to handle phase change. The model derivation will be reminded, as well as definition of the chemical potential and interpolation methods for kinematics viscosities and densities. Section 3 presents the Lattice Boltzmann schemes based on the Bhatnagar-Gross-Krook (BGK) collision operator for each equation. That collision operator is chosen because of its simplicity of implementation. Several improvements exist such as the two-relaxation-times (TRT) and the multiple-relaxation-times (MRT). Their benefits will be quickly discussed in that section. Computation of gradient and Laplacian operators that are involved in equations of phase-field and fluid flow will also be specified. Details on numerical implementation with the Kokkos library and various optimizations of LBM kernel will be discussed in Section 3.5. In Section 4, several basic code verifications are presented to check the implementation of each equation step-by-step. In Section 5, two purely qualitative simulations will be presented on the two-dimensional test case of film boiling. The first one will illustrate the capability of the model to simulate the detachment of bubbles on nodes and antinodes. The second one will illustrate the influence of the Jacob number on their detachment and shape. Here, we give a comparison of the code performance running on two architectures (CPU Intel and GPU NVIDIA ). Finally, Section 6 and three appendices will conclude this paper.

Two-phase flow with mass transfer

Phase change model

A single component fluid is considered, which can be either in a liquid (l) or gas (g) phase. The system is then composed with two incompressible fluids with constant densities ρ l and ρ g . A phase index φ ≡ φ(x, t) is introduced which can vary between 0 and 1 with φ = 0 (respectively φ = 1) corresponding to fluid l (resp. g) which is characterized by its density ρ l (resp. ρ g ) and its kinematic viscosity ν l (resp. ν g ). All other values of φ represent the interfacial zone or a mixture of both fluids l and g. When 0 < φ < 1, the densities ρ(φ) and the kinematic viscosities ν(φ) are respectively interpolated by

ρ(φ) = φ(x, t)ρ g + (1 -φ(x, t))ρ l , (1a)
ν(φ) = ν l ν g φ(x, t)ν l + (1 -φ(x, t))ν g . ( 1b
)
Local densities depending on position and time are noted ρχ (for χ = g, l) and write ρg (x, t) = ρ g φ(x, t) and ρl (x, t) = (1φ(x, t))ρ l . The total density writes ρ(x, t) = ρ g φ(x, t) + (1φ(x, t))ρ l . The method of harmonic mean is used in this work to interpolate the viscosity (Eq. ( 1b)) for simulating flows with viscosity contrast ([29, Eq. (29c)]). A comparison of both interpolation methods (linear and harmonic mean) is presented on the double-Poiseuille flow in Section 4.1. The local velocity u χ of each component χ is related to the volume averaged velocity u, the constant bulk density value ρ χ , and the volume diffusive flow rate j χ by [START_REF] Lee | Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces[END_REF] ρ χ j χ = ρχ (u χ -u) i.e. ρχ u χ = ρχ u + ρ χ j χ . The mass balance equations for each phase g and l writes

∂ ρg ∂t +∇ ∇ ∇ • ( ρg u + ρ g j g ) = + ṁ , (2a)
∂ ρl ∂t +∇ ∇ ∇ • ( ρl u + ρ l j l ) = -ṁ , ( 2b
)
where ṁ is the volumic production term (+) or sink term (-) due to phase change. Its physical dimension is M.L -3 .T -1 and its computation will be discussed in Section 2.3. In Eqs. (2a) and (2b), signs are chosen such as the phase change produces gas phase g to the detriment of liquid phase l. The mass flux relative to advection in each phase is ρχ u. In interfacial region, the mass flux ρ χ j χ has a diffusive origin and results of a regular transition of composition between two phases. By expressing Eqs. (2a) and (2b) with respect to φ(x, t) and assuming that the fluxes j g and j l are identical and opposite, j = j g = -j l , the following equations are obtained:

∂φ ∂t +∇ ∇ ∇ • (uφ + j) = + ṁ ρ g , (3a)
∂(1 -φ) ∂t +∇ ∇ ∇ • (u(1 -φ) -j) = - ṁ ρ l , (3b)
which after summing yield

∇ ∇ ∇ • u = ṁ 1 ρ g - 1 ρ l . (4)
To derive the interface tracking equation, in references [START_REF] Lee | Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces[END_REF][START_REF] Safari | Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow[END_REF] the flux j is assumed to be given by the Cahn-Hilliard flux defined by j = -M φ ∇ ∇ ∇µ φ where µ φ is the chemical potential. In that case Eq. (3a) becomes the CH equation with a source term of production in the second member. The Navier-Stokes/Cahn-Hilliard (NS/CH) model is very popular for simulations of two-phase flow since more than twenty years (e.g. without LBM [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF][START_REF] Jasnow | Coarse-grained description of thermo-capillary flow[END_REF] and [START_REF] Kendon | Inertial effects in three-dimensional spinodal decomposition of a symmetric binary fluid mixture: a lattice Boltzmann study[END_REF][START_REF] Zheng | A lattice boltzmann model for multiphase flows with large density ratio[END_REF][START_REF] Lee | Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces[END_REF][START_REF] Zu | Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts[END_REF] with LBM). However the chemical potential can be interpreted as the product of surface tension σ and curvature κ (see details in Section 2.2), and the CH equation imposes in its formulation a motion due to σ and κ even without coupling with a fluid flow. Here, in order to eliminate the curvature-driven interface motion inside the phase-field equation, we assume that the flux is defined by [START_REF] Sun | Sharp interface tracking using the phase-field equation[END_REF][START_REF] Chiu | A conservative phase field method for solving incompressible two-phase flows[END_REF] j = -M φ (∇ ∇ ∇φ -4φ(1φ)n/W ) and Eq. (3a) becomes the conservative Allen-Cahn (CAC) model with a source term:

∂φ ∂t +∇ ∇ ∇ • (uφ) = ∇ ∇ ∇ • M φ ∇ ∇ ∇φ - 4 W φ(1 -φ)n + ṁ ρ g . (5)
In Eq. ( 5), M φ is the interface mobility, W is the diffuse interface width and

n = ∇ ∇ ∇φ ∇ ∇ ∇φ (6)
is the unit normal vector at the interface directed from liquid toward gas. Eq. ( 5) is the Conservative version of Allen-Cahn (CAC) equation with a source term for modeling interface tracking with phase change. The accuracy of the phase-field simulations depends on two parameters: the interface thickness W and the mobility M φ . In reference [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]Sec. 5], a discussion is given regarding the numerical convergence of the phase-field method and the choice of those parameters in relation to the discretization step δx. For the Cahn-Hilliard equation, the mobility affects the thickness and perturbation magnitude of the chemical potential boundary layers. Here, for simulations of film boiling, preliminary sensitivity tests are performed on M φ and some details of its effects will be given in Section 5. The choice of ṁ will be discussed in Section 2.3. In the original paper [START_REF] Sun | Sharp interface tracking using the phase-field equation[END_REF], this equation is derived by assuming that the total advection velocity is composed of two terms: the external advective velocity u, plus the normal velocity to the interface u n n. That velocity u n is also defined as the sum of one term depending on the curvature κ, plus one independent of κ: u n n = ( ṽ -M φ κ)n. In the right-hand side of Eq. ( 5), the first term ∇ ∇ ∇ • j is an equivalent expression to the curvature term that is corrected with a "counter term" -M φ κ ∇ ∇ ∇φ [START_REF] Folch | Phase-field model for Hele-Shaw flows with arbitrary viscosity contrast. I. Theoretical approach[END_REF], in order to cancel the curvature-driven interface motion. The derivation is reminded in AppendixA by using the usual definition of curvature κ = ∇ ∇ ∇ • n with n defined by Eq. ( 6), and introducing the kernel function

φ = 1 2 1 + tanh 2ζ W (7)
in order to give an expression of ∇ ∇ ∇φ (see Eq. (A.7) in AppendixA):

∇ ∇ ∇φ = 4 W φ(1 -φ). (8)
That choice of kernel function imposes bulk phases for φ = 0 and φ = 1. Similar reasoning that cancels the curvature term can be found in [START_REF] Jamet | Thermodynamically consistent picture of the phase-field model of vesicles: Elimination of the surface tension[END_REF] in order to eliminate effects of surface tension (inherent in phase-field models) for membranes embedded in a Newtonian fluid. Let us notice that in this work the standard convention 0 ≤ φ ≤ 1 is used. Other conventions are possible, particularly when studying two-phase flow with high density ratio e.g. -φ ≤ φ ≤ φ where φ is defined by ρ g and ρ l (e.g. [27, Eq. (31)]). More generally, the inequality φ l ≤ φ ≤ φ g can be chosen. In that case the kernel function (Eq. ( 7)) and the expression of ∇ ∇ ∇φ must change. Moreover the source term in Eq. ( 5) must be modified by (see [12, Eq. (188)]): ṁ (φ g /ρ gφ l /ρ l ). Here, that expression is simplified to ṁ /ρ g with the standard choice φ g = 1 and φ l = 0.

The temperature equation is derived from the conservation law of total enthalpy ρH where H is the enthalpy (physical dimension E.M -1 where E is used for Energy) as carried out in crystal growth simulations [START_REF] Kobayashi | Modeling and numerical simulations of dendritic crystal growth[END_REF]:

∂(ρH ) ∂t +∇ ∇ ∇ • (uρH ) = ∇ ∇ ∇ • (K ∇ ∇ ∇T ) (9)
where the diffusive flux is given by the Fourier's law j T = -K ∇ ∇ ∇T with T being the temperature and K the thermal conductivity (physical dimension E.T -1 .L -1 .Θ -1 ). The enthalpy is defined by H = C p T + φL where C p is the specific heat (E.M -1 .Θ -1 ) and L is the latent heat of phase change (E.M -1 ). With this relation, enthalpies of liquid and gas are respectively equal to H l = C p T for φ = 0 and H g = C p T + L for φ = 1. With those notations and definitions the heat equation for temperature writes

∂T ∂t +∇ ∇ ∇ • (uT ) = α∇ ∇ ∇ 2 T - L C p ∂φ ∂t +∇ ∇ ∇ • (uφ) , (10)
where α = K /(ρC p ) is the thermal diffusivity, the second term in the right-hand side of Eq. ( 10) is interpreted as the release (or production) of latent heat during the displacement of the interface. When u = 0 the movement of the interface is only due to phase change between liquid and gas. Solving only Eq. ( 5) and [START_REF] Eichstädt | Accelerating high-order mesh optimisation with an architecture-independent programming model[END_REF] must be equivalent to solve the Stefan problem of phase change (see validation of Section 4).

Finally, the complete model of two-phase flows with phase change writes:

∇ ∇ ∇ • u = ṁ 1 ρ g - 1 ρ l , (11a)
∂(ρu) ∂t +∇ ∇ ∇ • (ρuu) = -∇ ∇ ∇p +∇ ∇ ∇ • η ∇ ∇ ∇u +∇ ∇ ∇u T + F tot , (11b)
∂φ ∂t +∇ ∇ ∇ • (uφ) = ∇ ∇ ∇ • M φ ∇ ∇ ∇φ - 4 W φ(1 -φ)n + ṁ ρ g , (11c)
∂T ∂t +∇ ∇ ∇ • (uT ) = α∇ ∇ ∇ 2 T - L C p ∂φ ∂t +∇ ∇ ∇ • (uφ) . (11d)
Eqs. (11a) and (11b) are the Navier-Stokes equations for modeling two Newtonian and incompressible fluids. In those equations p is the pressure, ρ(φ) is the density depending on the phase-field φ and η(φ) is the dynamic viscosity. F tot is the total force term defined as:

F tot = F s + F v ( 12
)
where F s is the surface tension force that is defined in the next subsection. The volumic force F v is the buoyancy force. Among different formulations of that force [START_REF] Amaya-Bower | Single bubble rising dynamics for moderate Reynolds number using Lattice Boltzmann Method[END_REF]Sec. 3.7], in this work the buoyancy is defined such as

F v = (ρ l -ρ(φ))g
where g is the constant acceleration due to the gravity. With that formulation, the gravity acts only on the gas phase for simulations of film boiling in Section 5.

Chemical potential and Cahn-Hilliard equation

The surface tension force F s is expressed here in its potential form [START_REF] Jacqmin | Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling[END_REF]:

F s = µ φ ∇ ∇ ∇φ (13)
where µ φ is the chemical potential which is defined as the change of free energy for a small variation of local composition of mixture: µ φ = δF /δφ. When the free energy is defined such as

F (φ) = ´v[V (φ) + K |∇ ∇ ∇φ| 2 /2]dv with V (φ) = Hφ 2 (1 -φ) 2 ,
the chemical potential writes

µ φ = 4Hφ(φ -1) φ - 1 2 -K∇ ∇ ∇ 2 φ. ( 14
)
The first term of the right-hand side of Eq. ( 14) is the derivative of V (φ) with respect to φ and the second term comes from the gradient energy term. The double-well ensures minima at φ = 0 and φ = 1. Coefficient H is the height of double-well and K is the gradient energy coefficient. It is well-known that the one-dimensional solution at equilibrium (i.e. µ φ = 0) of Eq. ( 14) is the hyperbolic tangent function defined by Eq. [START_REF] Obrecht | Multi-GPU implementation of the lattice Boltzmann method[END_REF]. A dimensional analysis of F (φ) indicates that H has the dimension of energy per volume unit, whereas K has the dimension of energy per length unit. In this formalism, the surface tension σ and the diffuse interface width W are proportional to the product and the ratio of both coefficients:

σ = 1 6 √ 2KH and W = 8K H (15a)
We also note that √ KH is homogeneous to an energy per surface unit which corresponds to the physical dimension of surface tension. The term K/H is homogeneous to a length as expected for the interface thickness. For the simulations of section 4, values of σ and W will be set and K and H will be derived by inverting those two relationships:

K = 3 2 W σ and H = 12 σ W . ( 15b
)
Let us notice that, if we use Eqs. ( 14) and (15b), the surface tension force F s = µ φ ∇ ∇ ∇φ can be written as

µ φ ∇ ∇ ∇φ = -(3/2)W σ ∇ ∇ ∇ 2 φ -16φ(1 -φ)(1 -2φ)/W 2 ∇ ∇ ∇φ.
The term inside the brackets is the curvature term κ ∇ ∇ ∇φ provided that the kernel function Eq. ( 7) is used for the second term (see Eq. (A.9) in AppendixA). In that case, the surface tension σ and the curvature κ appear explicitly in the definition of the chemical potential µ φ and the surface tension force is

F s = µ φ ∇ ∇ ∇φ = -(3/2)W σκ ∇ ∇ ∇φ ∇ ∇ ∇φ. Besides, if we set K = ε 2 and H = 1/4 in Eq. (15a), then we find (3/2)W = 6 √ 2ε. The surface tension force is F s = -σ(6 √ 2ε)(∇ ∇ ∇ • n) ∇ ∇ ∇φ ∇ ∇ ∇φ
which is the same relation in [44, Eq. ( 13)] provided that the kernel function Eq. ( 7) is applied for κ. As mentioned earlier, when the diffusive flux is proportional to the gradient of the chemical potential, then the evolution of φ follows the Cahn-Hilliard equation:

∂φ ∂t +∇ ∇ ∇ • (uφ) = ∇ ∇ ∇ • (M φ ∇ ∇ ∇µ φ ), (16)
with µ φ defined by Eq. ( 14). Compared to the standard CH equation, the main advantage of the conservative Allen-Cahn model lies in the computation of the right-hand side term. Indeed, the CH equation involves a fourth-order derivative because the flux is assumed to be proportional to gradient of chemical potential. A first Laplacian appears in Eq. ( 14) and a second one appears in the conservative equation Eq. ( 16). In the conservative Allen-Cahn equation (Eq. (11c)), only the second-order derivative is involved in its definition.

Production rate ṁ

Interface velocity of phase change

In sharp interface methods, the surface production rate ṁ (physical dimension M.L -2 .T -1 ) occurs on the separation area between liquid and gas. It is usually defined by [START_REF] Delhaye | Jump conditions and entropy sources in two-phase systems. Local instant formulation[END_REF][START_REF] Juric | Computations of boiling flows[END_REF]

ṁ = ρ g (u g -V I ) • n = ρ l (u l -V I )
• n where V I is the velocity of the interface, and u l and u g are respectively the velocities on liquid and gas sides. This relation is derived by integrating the mass conservation across the interface. Integration of the energy conservation yields an additional relation on ṁ which can be calculated in its simplest form by the difference of heat fluxes, ṁ

= ( K ∇ ∇ ∇T | l -K ∇ ∇ ∇T | g ) • n/L.
The driving force of evaporation is the heat quantity which is transferred at the interface. In [START_REF] Safari | Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow[END_REF], the liquid is assumed to be at saturation temperature T sat and in that case, only the heat quantity of the gas is considered and the temperature equation is solved only in the gas phase. Because of the diffuse interface, the rate ṁ is transformed to a volumic quantity ṁ by

ṁ = ṁ ∇ ∇ ∇φ = K ∇ ∇ ∇T • ∇ ∇ ∇φ/L
where φ follows the Cahn-Hilliard equation. The model was extended in [START_REF] Safari | Consistent simulation of droplet evaporation based on the phase-field multiphase lattice boltzmann method[END_REF] to include the gradient of the vapor concentration at the liquid-vapor interface as the driving force for vaporization. The model [START_REF] Safari | Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow[END_REF] was also applied in [START_REF] Begmohammadi | Simulation of pool boiling and periodic bubble release at high density ratio using lattice Boltzmann method[END_REF] to simulate nucleate pool boiling, including the bubble growth on and periodic departure from a superheated wall. Several other popular mass transfer models are reviewed in [39, Section 4.2] for phase change simulations.

Here, we notice that the source term ṁ /ρ g in Eq. (11c) can be identified as the normal velocity of the interfaceṽ ∇ ∇ ∇φ (see Eq. (A.8) in AppendixA) i.e. ṁ /ρ g =ṽ (because ṁ = ṁ ∇ ∇ ∇φ ). In Eq. (11c), the total velocity is the sum of an external velocity u plus the interface normal velocity. The latter has also been separated into one velocity depending on the curvature -M φ κ (which has been canceled) plus one velocity ṽ independent of the curvature. That velocity is responsible for the displacement of the interface because of the phase change. Its expression can be approximated by [30, Eq. (A.5)]:

ṽ = α A θ I -θ W , (17)
where θ is the dimensionless temperature defined as θ = (C p /L)(T -T sat ), θ I is the dimensionless interface temperature and A is a constant of proportionality that will be specified in section 2.3.2. Finally, if the kernel function

∇ ∇ ∇φ = (4/W )φ(1 -φ)
is used (see Eq. ( 8)), the source term ṁ /ρ g in Eq. (11c) takes the form ṁ

ρ g = -ṽ ∇ ∇ ∇φ = - 4α A W 2 (θ I -θ)φ(1 -φ). (18)

Value of coefficient A

In order to derive the value of A in Eq. ( 18), we proceed by analogy with the model of phase change for solidification and crystallization [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF]. First, Eq. (11c) with Eq. ( 18) are re-written in order to make appear the derivatives of the doublewell potential f (φ) and the interpolation function p(φ). Those functions are used in the solidification models derived from variational formulation based on the minimization of free energy [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF]. The interface is tracked by Eq. (11c) by assuming that the movement due to curvature is cancelled. That equation can be re-written (see AppendixA):

∂φ ∂t +∇ ∇ ∇ • (uφ) = M φ ∇ ∇ ∇ 2 φ - ∇ ∇ ∇φ •∇ ∇ ∇ ∇ ∇ ∇φ ∇ ∇ ∇φ -M φ κ ∇ ∇ ∇φ - 4α A W 2 (θ I -θ)φ(1 -φ). (19)
If the interface temperature is considered at saturation (i.e. θ I = 0), the source term is simplified to

(4α/A W 2 )θφ(1 -φ).
With the kernel function Eq. ( 7), the second term in the brackets writes (see Eq. (A.9

)) ∇ ∇ ∇φ • ∇ ∇ ∇ ∇ ∇ ∇φ / ∇ ∇ ∇φ = (16/W 2 )φ(1 - φ)(1 -2φ).
That term is proportional to the derivative (with respect to φ) of a double-well potential defined by

f (φ) = Hφ 2 (1 -φ) 2 with H = 1, hence ∇ ∇ ∇φ •∇ ∇ ∇ ∇ ∇ ∇φ / ∇ ∇ ∇φ = (8/W 2 )∂ f /∂φ. Besides if we set K ≡ ε 2 ,
then the two relationships Eqs.

(15b) with H = 1 yields ε 2 = W 2 /8. We also set M φ = ε 2 /T where T is the kinetic time, then Eq. ( 19) becomes

T ∂φ ∂t +∇ ∇ ∇ • (uφ) = ε 2 ∇ ∇ ∇ 2 φ - ∂ f ∂φ -ε 2 κ ∇ ∇ ∇φ - 4T α A W 2 (θ I -θ) ∂p ∂φ . ( 20
)
In the right-hand side of Eq. ( 20), the second term is the derivative of the double-well and the third term is the counter term. The last term is the coupling with temperature which involves the derivative (with respect to φ) of an interpolation function defined as p(φ) = φ 2 /2φ 3 /3. The factor 4 comes from the choice a = 1/2 in the kernel function (Eq. (A.6)) and we set W 0 = W /2. If we compare the coupling term of reference [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF] with the last term of Eq. ( 20), we can identify

λ = T α A W 2 0 , (21)
where W 2 0 = W 2 /4 and λ is the coupling coefficient in solidification/crystallization phase-field models. The star of λ means it is the particular value of λ that cancels the kinetic coefficient in the Gibbs-Thomson condition recovered by the matched asymptotic analysis of the phase-field model. Hence, that coupling term (Eq. ( 21)) means this is the particular model of phase change which cancels the kinetic coefficient in the Gibbs-Thomson equation. Besides, the curvature term is also removed by the counter term -ε 2 κ ∇ ∇ ∇φ . Finally, the coefficient A is identified to the coefficient a 2 in reference [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF]. Its value is a 2 = 0.6267 when the phase-field varies between -1 ≤ φ ≤ +1 and when the derivative of the interpolating function of temperature is p φ (φ) = 1φ 2 (the index φ indicates the derivative with respect to φ). In the present paper, the phase-field φ varies between 0 and 1 and the derivative of the polynomial function is p φ = φ(1φ). Because of those differences, the value of A must be computed from integrals obtained from the matched asymptotic expansion of the phase-field model. In AppendixB, details are given to obtain A = 10/48 ≈ 0.21, value that will be used for all simulations of this paper.

Lattice Boltzmann schemes

In this Section, we detail the lattice Boltzmann methods that are used to simulate the phase change model of Section 2 composed of Eqs (11a)-(11d) with Eq (13) for surface tension force and Eq. ( 18) for mass production rate. Simulations are performed by using three distribution functions ϑ i (x, t) ≡ ϑ i for ϑ = f , h, s where i = 0, ..., N pop and N pop is the total number of moving directions e i on a lattice (defined below). The first distribution function f i is used to recover the Navier-Stokes model (subsection 3.1); the second one g i is used for the phase-field equation (subsection 3.2) and the last one s i is used for the temperature equation (subsection 3.3). Each distribution function follows its own discrete lattice Boltzmann equation in which the collision term is considered with the Bhatnagar-Gross-Krook (BGK) approximation. That collision operator uses a unique relaxation parameter that is related to the diffusive parameter of the PDE (kinematic viscosity, mobility or diffusion coefficient). Several improvements exist such as the TRT [START_REF] Ginzburg | Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation[END_REF] or MRT [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF][START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] collision operators. They both use additional relaxation parameters (only one for TRT). With MRT, some of them can be related to physical parameters (e.g. anisotropic diffusion coefficient for transport equation) and the other ones control the stability of the algorithm when increasing the Reynolds number or Péclet number. Hence a wider range of parameters can be reached when simulations are performed with TRT and MRT. Let us mention that other alternatives exist in the literature (entropic, central moments, cumulants, ...) but an in-depth discussion of their benefits and drawbacks is out of the scope of this work. In Eq. ( 25), each discrete Boltzmann equation is expressed in terms of new variables f i , g i and s i , each one of them being defined by an appropriate variable change [START_REF] He | A Novel Thermal Model for the Lattice Boltzmann Method in Incompressible Limit[END_REF] (see details in AppendixC):

ϑ i = ϑ i + δt 2τ ϑ ϑ i -ϑ eq i - δt 2 S ϑ i for ϑ = f , h, s, (22)
where τ ϑ and S ϑ i are respectively the collision time and the source term relative to the distribution function ϑ; δt is the time step and ϑ eq i is the equilibrium distribution function. Two other notations are introduced: τ ϑ and ϑ i . The first one is the dimensionless collision rate that is defined by τ ϑ = τ ϑ /δt for each ϑ. The second one is the distribution function that is obtained after the stages of collision and streaming: ϑ i ≡ ϑ i (x + c i δt, t + δt). The use of this variable change (Eq. ( 22)) modifies the calculation of the zeroth-order moment M ϑ 0 of the distribution function ϑ i by (see AppendixC)

M ϑ 0 = ∑ i ϑ i + δt 2 S ϑ i for ϑ = f , h, s. (23)
It is also useful to introduce the variable change for the equilibrium function (see AppendixC.1)

ϑ eq i = ϑ eq i - δt 2 S ϑ i for ϑ = f , h, s, (24)
so that, with all those notations, the lattice Boltzmann equation writes

ϑ i = ϑ i - 1 τ ϑ + 1/2 ϑ i -ϑ eq i + S ϑ i δt (25)
for each distribution function ϑ = f , h, s. Before defining the equilibrium distribution functions and source terms, several lattices are introduced. In this work, the D2Q9 lattice and three 3D lattices are used: D3Q7, D3Q15 and D3Q19 (Fig. 1). For D2Q9 the moving vectors are defined by e 0 = (0, 0), e 1,3 = (±1, 0), e 2,4 = (0, ±1), e 5,6 = (±1, 1) and e 7,8 = (∓1, -1). for 3D lattices, the moving vectors e i are defined such as e 1 = (1, 0, 0) T , e 2 = (0, 1, 0) T , . . ., e 6 = (0, 0, -1) T for D3Q7 (Fig. 1a). For D3Q15, additional diagonal vectors are defined such as (see Fig. 1b) e 7 = (1, 1, 1) T , e 8 = (-1, 1, 1) T , . . ., e 14 = (1, -1, -1) T . Finally for D3Q19 (Fig. 1c):

Incompressible Navier-Stokes

Several lattice Boltzmann schemes exist for incompressible version of Navier-Stokes equations. The fully incompressible condition has already been proposed in literature but necessitates to solve an additional Poisson equation [START_REF] Inamuro | A lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF] or an additional predictor-corrector step [START_REF] Zu | Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts[END_REF]. Here we prefer to apply the artificial compressibility method [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] for which the solenoidal condition ∇ ∇ ∇ • u = 0 is approximated by (1/β)∂p/∂t + ∇ ∇ ∇ • u = 0 where β is the artificial compressibility coefficient. In LB framework, the method was derived in [START_REF] He | Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation[END_REF] with β = ρ 0 c 2 s where ρ 0 is the constant density of bulk phase. The LB scheme writes

f i = f i - 1 τ f + 1/2 f i -f eq i + S f i δt, (26a)
f eq i = w i p + ρ(φ)c 2 s c i • u c 2 s + (c i • u) 2 2c 4 s - u • u 2c 2 s , (26b)
with

f eq i = f eq i -S f
i δt/2 and ρ(φ) is given by Eq. (1a). In Eq. (26a) τ f is the collision rate which is related to the kinematic viscosity by ν = τ f c 2 s δt. Hence, the collision rate is obtained by τ f (φ) = 3ν(φ)(δt/δx 2 ) with the kinematic viscosity ν(φ) interpolated by Eq. (1b). In Eq. (26a), the source term S f i contains contributions of external forces (involving F tot ) plus the production term in mass conservation (involving ṁ ):

S f i = F f i + P f i (27a)
with [START_REF] Mohammadi-Shad | Phase-field lattice Boltzmann modeling of boiling using a sharp-interface energy solver[END_REF] F

f i = (c i -u) • (Γ i -w i )∇ ∇ ∇ρ(φ)c 2 s + Γ i F tot , (27b)
P f i = w i ρc 2 s ṁ 1 ρ g - 1 ρ l . (27c)
In Eq. (27b), F tot is the external force defined by Eq. ( 12) and the function Γ i ≡ Γ i (u) is defined by:

Γ i = w i 1 + c i • u c 2 s + (c i • u) 2 2c 4 s - u • u 2c 2 s . ( 27d
)
After the stages of collision and streaming, the first-order moment (momentum) and the zeroth-order moment (pressure) are updated by [START_REF] Mohammadi-Shad | Phase-field lattice Boltzmann modeling of boiling using a sharp-interface energy solver[END_REF]

ρu = 1 c 2 s ∑ i f i c i + δt 2 F tot , (28a)
p = ∑ i f i + δt 2 u •∇ ∇ ∇ρc 2 s + ρc 2 s ṁ 1 ρ g - 1 ρ l . (28b)

Conservative Allen-Cahn model

The lattice Boltzmann equation for the conservative Allen-Cahn model acts on the distribution function g i . The evolution equation is

g i = g i - 1 τ g + 1/2 g i -g eq i + S g i δt, (29a)
g eq = φΓ i , (29b)
with the variable change g eq i = g eq -δtS g i /2. The mobility coefficient is related to the collision rate by M φ = τ g c 2 s δt. The source term S g i contains two contributions:

S g i = F g i + P g i , (30a)
where the first one F g i involves the counter term with the normal vector n [START_REF] Mitchell | Development of a three-dimensional phase-field lattice Boltzmann method for the study of immiscible fluids at high density ratios[END_REF], and the second one P g i involves the mass production term ṁ :

F g i = 4 W φ(1 -φ)w i c i • n and P g i = w i ṁ ρ g . ( 30b
)
Let us notice that the scheme is equivalent (see AppendixC.2) to the lattice Boltzmann equation

g i = g i - 1 τ g + 1/2 g i -g eq,CAC i + P g i δt (31a)
where only the source term P g i appears in the source term and the equilibrium distribution function is redefined as [START_REF] Fakhari | A weighted multiple-relaxation-time lattice Boltzmann method for multiphase flows and its application to partial coalescence cascades[END_REF] g eq,CAC i

= φΓ i + M φ 4 W φ(1 -φ)w i c i • n c 2 s (31b) with g eq,CAC i = g eq,CAC i
-δtP g i /2. After the stages of collision and streaming, the new phase-field is obtained by the zeroth-order moment of g i which must be corrected with the production term:

φ(x, t) = ∑ i g i + δt 2 ∑ i P g i . (32)
This relation holds for both formulations that use g eq i and g

eq,CAC i

because ∑ i F g i δt/2 = 0.

Temperature equation

The lattice Boltzmann scheme for temperature equation writes:

s i = s i - 1 τ s + 1/2 s i -s eq i + S s i δt (33a)
s eq i = T Γ i (33b)
where the thermal diffusivity α is related to the collision rate by α = τ s c 2 s δt. The source term S s i is defined such as:

S s i = F s i + P s i ( 33c
)
where

F s i = w i L C p ∇ ∇ ∇ • (uφ) and P s i = w i L C p ∂φ ∂t (33d)
Finally, the new temperature is computed by

T = ∑ i s i - δt 2 L C p ∂φ ∂t +∇ ∇ ∇ • (uφ) . (34)
In Sections 4 and 5, simulations will be carried out with Dirichlet boundary conditions applied on temperature T and phase-field φ. In order to impose such a condition, for example on temperature T w on left boundary of a D2Q9 lattice, the unknown distribution functions s i | unknown are updated with the anti bounce-back method [START_REF] Ginzburg | Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations[END_REF]:

s i | unknown = -s i + 2w i T w
where i is the opposite direction of i.

Computations of gradients and Laplacian

The unit normal vector n and force term F s require computation of gradients. Moreover the chemical potential µ φ necessitates to calculate the Laplacian of φ. Gradients and Laplacian that are involved in definitions of n (Eq. ( 6)) and µ φ (Eq. ( 14)) are discretized by using the directional derivatives methods. The method has already demonstrated its performance for hydrodynamics problem in order to reduce parasitic currents for two-phase flow problem [START_REF] Lee | Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces[END_REF][START_REF] Lee | Eliminating parasitic currents in the lattice Boltzmann equation method for non ideal gases[END_REF][START_REF] Lee | Effects of incompressibility on the elimination of parasitic currents in the lattice Boltzmann equation method for binary fluids[END_REF]. The directional derivative is the derivative along each moving direction on the lattice. Taylor's expansion at second-order of a differentiable scalar function φ(x) at x + e i δx and xe i δx yields the following approximation of directional derivatives:

e i •∇ ∇ ∇φ x = 1 2δx [φ(x + e i δx) -φ(x -e i δx)] (35a)
The number of directional derivatives is equal to the number of moving direction e i on the lattice i.e. N pop . The gradient is obtained by

∇ ∇ ∇φ x = 3 N pop ∑ i=1 w i e i e i •∇ ∇ ∇φ x . (35b)
The three components of the gradient ∂ x φ, ∂ y φ and ∂ z φ are obtained by calculating each directional derivative e i • ∇ ∇ ∇φ x and next, by calculating the moment of first-order ∇ ∇ ∇φ x . For the calculation of ∇ ∇ ∇ 2 φ, all directions of propagation are taken into account by

(e i •∇ ∇ ∇) 2 φ x = 1 δx 2 [φ(x + e i δx) -2φ(x) + φ(x -e i δx)] . (36a)
The Laplacian is obtained by summing and weighting each term with

∇ ∇ ∇ 2 φ x = 3 N pop ∑ i=1 w i (e i •∇ ∇ ∇) 2 φ x . (36b)
Other approximations exist [START_REF] Lee | Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces[END_REF][START_REF] Lee | A stable discretization of the lattice boltzmann equation for simulation of incompressible two-phase flows at high density ratio[END_REF] such as the first-order and second-order upwind schemes (or biased differences) respectively defined by e i •∇ ∇ ∇ up 1 φ x = [φ(x + e i δx)φ(x)]/δx and e i •∇ ∇ ∇ up 2 φ x = [-φ(x + 2e i δx) + 4φ(x + e i δx) -3φ(x)] /(2δx). Here, by simplicity, the central difference approximation is applied for all simulations even though that approximation fails to capture the velocity profiles in low density regions [START_REF] Fakhari | Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios[END_REF] and biased directional derivatives can fix that issue [START_REF] Guo | Lattice Boltzmann Method and Its Applications in Engineering[END_REF]. Those biased differences could be tested in future works with LBM saclay.

Numerical implementation and kernel optimization

All LBM schemes of this Section were implemented in a new code called LBM saclay written in C++. The main advantage of this new code is its portability targeting all major HPC platforms and especially those based on GPU-and CPU-architectures. Actually, LBM saclay can run without modification on any architecture that Kokkos supports. The current compatibilities are indicated in [START_REF]Compatibilities of Kokkos library[END_REF] and summarized in Tab. 1. For more information, the reader can refer to the Kokkos documentation. Let us mention that the current support for AMD GPU is experimental through the C++ library HIP (Heterogeneous-Compute Interface for Portability) and it is planned to be supported at the end of 2020.

Two levels of parallelism are implemented in the code. The first one is the intra-node parallelism (shared memory) with the Kokkos library, an opensource C++ library with parallel algorithmic patterns and data containers. Specific commands of the Kokkos library optimize loops with OpenMP, Pthreads or CUDA during compilation. An example of using Kokkos' functionalities is presented on Fig. 2 to compute at each time-step the zeroth-order moment of a distribution function. The second level of parallelism is a standard domain decomposition performed with MPI: the full computational domain is cut into several sub-domains associated with each computational node (distributed memory).

When developing the code, several optimizations were implemented and compared in particular to enhance its performance on each architecture. The first way to consider the stages of collision and streaming is to "fuse" those two steps inside a single kernel, i.e. both stages are simply done in one single "for-loop" performed on the lattice nodes. The "fused" version does not require an intermediate memory load contrary to standard implementation for which both stages are well separated. However, the fused kernel contains more floating point operations per iteration of the "for-loop". This is a drawback if the number of floating point operations becomes large enough to exhaust the amount of register memory available on the architecture (this number is significantly lower on GPU than on CPU). So, if the register memory is full, additional variables will be allocated in the external DRAM memory, generating additional traffic on the memory bus and degrading performance. For GPUs NVIDIA and CPUs Intel Skylake, best performance is obtained with the fused version.

Alternatively, two optimizations were tested which are well suited for Intel KNL (KNights Landing) processors [START_REF] Calore | Early Experience on Using Knights Landing Processors for Lattice Boltzmann Applications[END_REF]: the first one is the "CSoA" optimization (Cluster of Structure of Array) i.e. for each line of the lattice, LBM nodes are stored in memory modulo M where typically M = 8 and each line is padded to be a multiple of M. The access of data container is done with data(iMem,j,k,ipop) where iMem is computed from the physical node location i. The CSoA optimization improves vectorization and memory alignment for streaming stage but performance decreases for large domain on D2Q9 lattice. The second optimization for KNL is "CSoA2", i.e. the population index ipop of data(i,j,k,ipop) is interverted to data(i,ipop,j,k), where i,j,k are indices of position. With this permutation, the memory locality is restored for the collision stage.

Comparisons were performed on a simplified diffusive problem. The CSoA2 optimization enhances performance on KNL processors, but on Fig. 3, we can see that it remains far below to that obtained on GPUs, even older generation GPUs (K80). Computational times are expressed in Million Lattice Updates per Second (MLUPS) as an effective metric measuring the number of millions of node lattice update per seconds. That performance metric is used by node so that it is independent of the type of lattice (e.g. D3Q7, D3Q15, D3Q19). With that metric, a larger lattice will give a smaller MLUPS. In the rest of this paper, most of validations and simulations of Sections 4 and 5 are carried out on GPUs. In Section 5.3, comparisons of computational times on GPU and CPU will be presented on the test case of film boiling for two mesh sizes.

Finally, let us mention that all kernels (Navier-Stokes, phase-field and temperature equations) have been developed in 2D and in 3D. They all run in 3D separately. However, all coupling terms, i.e. the surface tension force (Eq. ( 13)), the chemical potential (Eq. ( 14)) and the advective term in Eq. (33d)) were developed and checked only in 2D. Hence, verification of couplings in Section 4 and film boiling simulations of Section 5 will be presented only in 2D. The three-dimensional extension of coupling terms is planned for future works.

Code verifications

In this section, the numerical implementation of the LBM schemes of Section 3 is checked by comparison with wellknown solutions. Validations are gathered into two parts in order to check implementations step-by-step. In subsection 4.1, verifications are done without phase change, i.e. by neglecting the temperature equation and by assuming that the mass transfer is zero ( ṁ = 0 in Eq. (11a) and (11c)). The conservative Allen-Cahn model, and the coupling with fluid flow are verified successively. In subsection 4.2, the phase change model is checked by considering the phase-field equation coupled with temperature. The LBM code is compared with an analytical solution of Stefan's problem with two different diffusivities.

Verifications without phase change

We first compare implementation of the conservative Allen-Cahn model on two test cases: Zalesak's slotted disk and interface deformation inside a vortex. Next the coupling with Navier-Stokes model will be considered with the layered Poiseuille flow and the Laplace law.

Verifications of the phase-field model

Two verifications of phase-field implementation are presented. In the first one, we check that the contour of a slotted disk is well conserved inside a rotating fluid [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF]. In the second one, we check that the simulation retrieves a circle when an initial disk is deformed inside a vortex that changes its direction of rotation over time. For both simulations, the mesh is composed of 201 × 201 × 3 nodes with periodic boundary conditions applied on all faces, the time-step is δt = 10 -4 and the space-step δx = 5 × 10 -3 .

Zalesak's slotted disk. Inside a domain of lengths L x = L y = 1, and L z = 0.01, a disk is initialized at the center of the domain

x c = (100, 100, 1) T by φ(x, 0) = 1 + tanh (R -d c )/ √ 2W 0 /2 with d c = (x -x c ) 2 + (y -y c ) 2 + (z -z c ) 2
, W 0 = 2 and R = 80 l.u. (lattice units). The diffuse disk is slotted by imposing φ(x, 0) = 0 if x c -R/6 ≤ x ≤ x c + R/6 and y c -1.1R ≤ y ≤ y c . Components of velocity are imposed by u x (x) = u 0 (2y -1), u y (x) = u 0 (1 -2x) and u z (x) = 0. The value of u 0 is chosen such that the slotted disk performs one complete rotation at T f = 4, i.e. u 0 = 0.7853975 and both parameters of CAC model are set as M φ = 5 × 10 -4 and W = 6δx. The rotation of the slotted disk is presented on Fig. 4 where the interface position φ = 1/2 is superimposed to the initial condition at four times. At the final time of simulation t = T f (Fig. 4d), the contour φ = 0.5 (red) is superimposed to the initial one (black) although the slot corners are slightly rounded. Layered Poiseuille flow. The Navier-Stokes implementation is checked with the analytical solution of a layered Poiseuille flow [START_REF] Zu | Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts[END_REF] for two fluids named A and B:

u x (y) =    Gh 2 2η A -y h 2 -y h η A -η B η A +η B + 2η A η A +η B (-h ≤ y ≤ 0) Gh 2 2η B -y h 2 -y h η A -η B η A +η B + 2η B η A +η B (0 ≤ y ≤ h) (37)
where η A and η B are the dynamic viscosities and 2h is the channel width. The pressure gradient is defined by G = u c (η A + η B )/h 2 with u c = 5 × 10 -5 . For the LB simulation, the mesh is composed of 101 × 101 × 3 nodes and the pressure gradient is replaced by a force term defined by F = (G, 0, 0) T . Periodic boundary conditions are set for all limits except for planes of normal vector directed in y-direction where no-slip conditions are imposed with the half bounce-back method. Two layers of different viscosity are defined as initial condition for φ: φ(x, 0) = 0.5{1 + tanh [2(yy 0 )/W ]} where W = 6δx controls the slope of the hyperbolic tangent function and y 0 = (y max + y min )/2. The mobility coefficient is M φ = 0.1. Comparisons between the LBM code and the analytical solution are presented for two cases. In the first one, the density is identical for both fluids (ρ A = ρ B = 1) and three viscosity ratios are checked on Fig. 6a: η B /η A = 1/3, 1/5, 1/10. For the first ratio ν A = 0.1 and ν B = 0.3; for the second one ν A = 0.07 and ν B = 0.35 and for the third one ν A = 0.01 and ν B = 0.1. For the second test case, the viscosity of each phase is set equal to ν A = ν B = 0.07 and three density ratios are checked on Fig. 6b: ρ A /ρ B = 1/1.658, 1/2, 1/3. The ratio 1/1.658 is used in the simulations of film boiling as well as the viscosity ratio ν A /ν B = 1/6. In Fig. 6c this viscosity ratio is checked for two cases. In the first simulation (red curve), the density ratio is equal to one and the viscosity is interpolated by two methods: the linear (black squares) and the harmonic mean (red circles) defined by

ν(φ) = [1 -φ(x, t)]ν A + φ(x, t)ν B , (38)
and Eq. (1b) respectively. The differences observed with the former method justify the choice of using the latter in the second simulation (blue curve) which combines both ratios of viscosity and density.

Laplace law. The two-dimensional Laplace law is checked by initializing a drop at the center of a square domain of length L x = L y = 2.56 discretized with 256 × 256 nodes. By varying the radius R, the difference between pressure inside the drop (p in ) minus the pressure outside (p out ) must vary proportionally with the surface tension σ: (i In order to check that relationship, an initial drop of radius R and surface tension σ is initialized at the center of the domain (x c = y c = 1.28). The density ratio ρ g /ρ l is set equal to two (ρ g = 2, ρ l = 1) and the viscosities are identical for each phase: ν l = ν g = 0.04. The interface parameters are M φ = 0.04 and W = 0.05 = 5δx. The LBM code is run with a time-step equal to δt = 10 -4 until the stationary solution is obtained. At the end of simulation, the difference between numerical pressures ∆p = p inp out is plotted for three values of surface tension σ = 0.04, 0.08, 0.15. For each value of surface tension, six LBM simulations are run for six values of radius corresponding to each dot on Fig. 6d. On that plot, the slopes of LBM vary linearly and fit quite well to the Laplace law.

p in -p out = σ R . ( 39
) (i) t = 0 (ii) t = T f /2 (iii) t = T f (
) t = T f /2 (ii) t = T f (iii) t = 3T f /2 (iv) t = 2T f (

Verifications with phase change: one-dimensional Stefan problem

In this section, we consider the problem of phase change without flow (u = 0). The objective is to validate the coupling between equations of phase-field and temperature. More precisely, we check the new approximation (Eq. ( 18)) of mass production rate ṁ in the phase-field equation (Eq. (11c)) and the latent heat release in the temperature equation (Eq. (11d)), i.e. the source term -∂φ/∂t. Validation is carried out with the Stefan problem for which several analytical solutions exist [START_REF] Hahn | Heat Conduction. Third Edition[END_REF]Chapter 12]. Here we consider one of the most general one-dimensional problem where the three unknowns are the interface position varying with time x I (t), the liquid temperature T l (x, t) and the gas temperature T g (x, t). Besides, the thermal diffusivities of each phase α l and α g can be different. The one-dimensional domain ]0, ∞[, is initially filled with gas with constant temperature T g (x, t) x>0,t=0 = T ∞ that is greater than the saturation temperature T sat . The left wall x = 0 is maintained at T w for t ≥ 0. As a result, condensation starts at the boundary x = 0 and the liquid-gas interface propagates in the positive direction. At x → ∞, the temperature is kept at T ∞ .

Analytical solutions. The mathematical formulation of this problem writes [66, Section 12-3]

∂T l ∂t = α l ∂ 2 T l ∂x 2 (40a)
for 0 < x < x I (t), with the left boundary condition imposed at T l (x, t)| x=0 = T w . The evolution of the gas phase is formulated as for

∂T g ∂t = α g ∂ 2 T g ∂x 2 (40b) 0 1x10 -5 2x10 -5 3x10 -5 4x10 -5 5x10 -5 6x10 -5 7x10 -5 8x10 -5 9x10
x I (t) < x < ∞ with T g (x → ∞, t) = T ∞ , with the initial condition T g (x, t = 0) = T ∞ and boundary condition T g (x → ∞, t) = T ∞ .
Interfacial conditions are specified by

T l (x, t)| x=x I (t) = T g (x, t)| x=x I (t) = T I , (40c)
K l ∂T l ∂x x=x I (t) -K g ∂T g ∂x x=x I (t) = ρL dx I (t) dt . (40d)
In Eq. (40d), K l and K g are the thermal conductivities of each phase. We consider identical specific heat C l p = C g p = C p and we set C p = 1, L = 1 and ρ = 1. Solutions of interface position and temperature profiles [66, p. 469] are

x I (t) = 2ξ √ α l t, (41a)
θ l (x, t) = θ w + (θ I -θ w ) erf(x/2 √ α l t) erf(ξ) , (41b)
θ g (x, t) = θ ∞ + (θ I -θ ∞ ) erfc(x/2 √ α g t) erfc(ξ α l /α g ) , (41c)
where the temperatures are re-written in dimensionless form with θ = C p (T -T sat )/L. When θ = 0 the temperature of system is at saturation temperature T sat and when θ > 0 (resp. θ < 0), the system is superheated (resp. undercooled). In Eqs.

(41a)-(40d), ξ is solution of the transcendental equation

e -ξ 2 erf(ξ) + α g α l 1/2 θ I -θ ∞ θ I -θ w e -ξ 2 (α l /α g ) erfc(ξ α l /α g ) = - ξ √ π θ w (41d)
where θ w in the right-hand side is the Stefan number defined by St = C p (T w -T sat )/L. Those solutions are compared with LBM saclay, first with identical thermal diffusivities α l = α g and an interface temperature θ I equals to zero. The second validation considers three ratios of diffusivity α j l /α j g (for j = 1, 2, 3) with an interface temperature which is different of the saturation one (θ I = 0). Data entry of LBM simulations. For LBM simulations, the two-dimensional D2Q9 lattice is used for the temperature and phase-field equations. The LBM computational domain is [ x , L x ] × [ y , L y ] = [0, 512] × [0, 32] which is discretized by N x × N y = 512 × 32 nodes i.e. δx = 1. The time-step is also set to δt = 1. Boundary conditions are periodic for s i and g i at y and L y (bottom and top walls respectively) and Dirichlet boundary conditions are applied on left (x = x ) and right (x = L x ) walls by anti-bounceback method on g i and s i . For phase-field, the Dirichlet boundary conditions are φ(x, t) x= x = 0 and φ(x, t) x=L x = 1. For the temperature equation, they are θ(x, t) x= x = θ w and θ(x, t) x=L x = θ ∞ . The temperature is initialized with θ(x, 0) = θ ∞ for 0 < x ≤ L x and the phase-field with φ(x, 0) = 0.5 [1 + tanh(2x/W )]. The mobility parameter is M φ = 0.08, the interface thickness is W = 3δx.

Validations for α l /α g = 1 and θ I = 0. Before considering the more general case α l /α g = 1 and θ I = 0, we assume that thermal diffusivities are the same in liquid and gas (α l = α g = α) and the interface temperature is at saturation (θ I = 0). In that case, whatever the diffusivity value α, the solution of the transcendental equation (Eq. (41d)) depends only on θ w and θ ∞ . With θ w = -0.3 and θ ∞ = 0.3, its solution is ξ = 0.280680. Comparisons between analytical solutions and LBM simulations are presented on Fig. 7a for three values of thermal diffusivity α j g = 0.14, 0.08, 0.03 with j = 1, 2, 3. LBM temperature profiles are superimposed with the analytical solution (Eqs. (41b) and (41c)) at the final time of simulation t f = 2 × 10 5 (Fig. 7a,left). Successive positions of vapor/liquid interface also fit with the analytical solution (Fig. 7a,right) for three values of thermal diffusivity.

Validations for α l /α g = 1 and θ I = 0. Now we consider a more general case for which the diffusivities of liquid and gas can be different. Three ratios are simulated α j l /α j g = 10, 5, 2 for j = 1, 2, 3 with α 1 l = 0.14, α 2 l = 0.125 and α 3 l = 0.08. Same values of θ w = -0.3 and θ ∞ = 0.3 are kept, and the interface temperature is now equal to θ I = 0.05. For those values, the corresponding solutions of the transcendental equation are ξ 1 = 0.349635, ξ 2 = 0.343882 and ξ 3 = 0.331864. For LBM simulations, all numerical values are identical except for interface temperature and diffusivities of each phase. As confirmed by temperature profiles (Fig. 7b, left) and the evolution of interface position (Fig. 7b, right), the model of phase change is well adapted to simulate the phase change problem with different diffusivities in each phase and an interface temperature not equal to zero. Finally this test case validates the approximation of the mass production rate ṁ defined by Eq. ( 18) and implementation of LBM for the phase-field and temperature equations.

Simulations of film boiling

Film boiling is a classical problem of two-phase flows with phase change. It has already been simulated with a lot of different numerical techniques (see [START_REF] Kharangate | Review of computational studies on boiling and condensation[END_REF] for a recent review) for studying the effect of geometries such as an horizontal cylinder [START_REF] Esmaeeli | A front tracking method for computations of boiling in complex geometries[END_REF] or for studying the effect of an electric field [START_REF] Pandey | Effect of superheat and electric field on saturated film boiling[END_REF]. With the lattice Boltzmann method, several simulations use the Cahn-Hilliard model or the pseudo-potential method (respectively in [69, 70, and references therein]). Here we present the capability of the conservative Allen-Cahn equation with a production rate defined by Eq. ( 18) to simulate that problem. In section 5.1, the physical configuration is reminded; in section 5.2 one simulation of bubbles detachment on nodes and anti-nodes is detailed; in section 5.3, indications will be given on computational times for two mesh sizes: 1024 2 for GPU and CPU and 4096 × 3072 for multi-GPUs.

Physical configuration

Inside a two-dimensional domain Ω = Π υ=x,y [ υ , L υ ], a thin film of gas of height y 0 is initialized near the bottom wall y = y which is heated by applying a constant temperature θ| y= y = θ w . The liquid is above the thin film and the gravity acts downward g = (0, -g y ) T . On the top wall y = L y , the temperature is imposed at saturation and the phase-field is equal to φ = +1 (i.e. gas phase). The left and right walls are periodic. If the interface is destabilized by an initial condition defined by

y = y 0 + y 1 sin 2πx λ , (42)
where y 1 and λ are respectively the amplitude and the wavelength of the perturbation, then we can observe bubbles of gas that grow, detach and rise in the domain, provided that the wavelength of perturbation λ is greater than a critical value λ c defined by

where the interface position y is perturbed with several modes λ i which are randomly picked, uniformly distributed between 0.5λ c ≤ λ i ≤ 1.5 √ 3λ c . We simulate two values of wall temperature θ w = 0.025 and θ w = 0.1 corresponding to Jacob numbers respectively equal to Ja = 0.025 and Ja = 0.1. All other values of physical parameters remain identical (Table 2). A comparison on shapes of bubbles is given at t = 95.78 on Fig. 9. When the Jacob number has the value of Section 5.2, discrete bubbles are released periodically from the initial condition (Fig. 9a). When the Jacob number is increased to 0.1, long vapor jets are observed below bubbles (Fig. 9b). That observation is consistent with those simulated with other techniques and even observed on experiments cited in [39, Sec 5.1.2 and Fig. 9]. The simulation took 80 minutes (713 MLUPS) on 8 parallel GPUs to complete 5.33 × 10 5 time iterations.

Conclusion

In this paper, the LBM implementation of two-phase flows was revisited by improving two main points. The first one focuses on the model formulation of phase change and the second one focuses on the portability of the code on various platforms. The interface is tracked by the conservative Allen-Cahn model with a source term involving a mass production rate at the interface. In this work, that source term is simplified compared to approaches of literature, and the approximation avoids to calculate the gradients of temperature numerically. The model is able to simulate two phases of different thermal diffusivities with an interface temperature which is not necessarily at saturation. The phase-field model is coupled with the incompressible Navier-Stokes model where a source term was added in the mass balance equation. The source term is defined as the product of mass production rate times one term inversely proportional to densities. An additional equation on temperature completes the model. The time derivative of phase-field appears in the source term of that equation. It is interpreted as the release or absorption of latent heat at the interface.

The Lattice Boltzmann schemes for all equations are implemented in a new C++ code coupled with the Kokkos library for its performance portability. The new code, called LBM saclay, can be run with good performance on several architectures such as Graphical Process Units (GPUs), Central Process Units (CPUs) and even multi-GPUs and multi-CPUs. Indeed, two levels of parallelism are developed inside the code. The first one uses Kokkos for intra-node parallelism, whereas MPI takes in charge the domain decomposition. Preliminary comparisons between GPUs and CPUs were carried out on a simple diffusive problem. As expected from literature, those tests show clearly that best performance is obtained with GPU compared to CPU (Skylake or KNL) even for best optimization of LBM kernels (CSoA2) which has been developed for Intel Skylake. Here, comparisons were performed with the same C++ source code. No low-level language (CUDA or OpenCL) was used for GPUs.

Numerical implementation was checked with several test cases to validate step-by-step the full model of fluid flows with phase change. The conservative Allen-Cahn equation is validated with two test cases: (i) Zalesak's slotted disk and (ii) interface deformation inside a vortex. The coupling with Navier-Stokes equations is also checked with two test cases: the layered Poiseuille flow and Laplace law. Next, the coupling between equations of phase-field and temperature were compared to the most general one-dimensional analytical solution of the Stefan problem. Comparisons were done first by assuming identical thermal diffusivities, and next by using various ratios of diffusivities with an interface temperature that is different of the saturation one. The full model was simulated on the test case of film boiling on one GPU and one multicore CPU for two mesh sizes. Computational times are clearly in favor of GPUs. Finally, the film boiling problem is simulated with 8 parallel GPUs for mesh size that is twelve times bigger than the previous one.

In this paper, foundations have been laid for improving performance of lattice Boltzmann simulations in a context of quick evolution of HPC platforms. In the future, a three-dimensional extension of the coupling terms is planned. Next LBM saclay could be enriched with other models requiring interface tracking such as crystal growth and demixing of ternary fluids. Besides, the range of physical parameters could be increased and the code stability could be enhanced by using alternative collision operators such as those based on the Two-Relaxation-Times and Multiple-Relaxation-Times. a thermodynamically derived phase-field model, such as the one used for computation of chemical potential (Eq. ( 14)) with bulk phases φ = 0 and φ = 1. The normal derivative of Eq. (A.6) leads to Eq. (A.8) is the Allen-Cahn equation for which the curvature-driven displacement of the interface has been canceled with a counter term. Let us notice that, if ṽ is chosen such as ṽ = α(θ Iθ)/(A W ) then ṽ ∇ ∇ ∇φ ≈ -(4α/A W 2 )(θ Iθ)φ(1φ) can be used in Eq. (A.8) for the problem of phase change. The release or absorption of latent heat at the interface is taken into account in the temperature equation by the time derivative of φ. If the physical problem necessitates a curvature-driven interface motion, the curvature term must be kept in the Allen-Cahn equation and then only the first term in the right-hand side of Eq. (A.3) appears in the derivation. With a = 1/2, the term ∇ ∇ ∇φ •∇ ∇ ∇ ∇ ∇ ∇φ / ∇ ∇ ∇φ is equal to

∇ ∇ ∇φ •∇ ∇ ∇ ∇ ∇ ∇φ ∇ ∇ ∇φ = ∂ 2 φ ∂ζ 2 = 16 W 2 φ(1 -φ)(1 -2φ) (A.9)
The curvatuve-driven term writes Those integrals can be computed analytically and yield a numerical value provided that the interface width W is set. Here, to be consistent with the rescaling of space and the analysis performed in [START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimensions[END_REF], it is enough to set W = 2 √ 2, and the integrals are:

M φ κ ∇ ∇ ∇φ = M φ ∇ ∇ ∇ 2 φ - 16
I = 1 3 √ 2 , J = 1

AppendixC. Discrete lattice Boltzmann equations

In this Appendix, the variable change for the discrete lattice Boltzmann equation is reminded in AppendixC.1. In Ap-pendixC.2, we will show that, for CAC model, the formulation with a source term is equivalent to the formulation with a modification of the equilibrium distribution function.

AppendixC.1. Variable change for discrete lattice Boltzmann equation

The discrete lattice Boltzmann equation with an external force or source term S ϑ i can be written with the BGK collision term:

∂ϑ i ∂t + c i •∇ ∇ ∇ϑ i = - ϑ i -ϑ eq i τ ϑ + S ϑ i . (C.1)
In what follows, the calculations will be performed by setting ϑ ≡ f , S ϑ i = S f i = S i and τ ϑ ≡ τ but the variable change derivation holds also for ϑ ≡ h and ϑ ≡ s. Terms that are evaluated at position x and time t are noted f i ≡ f i (x, t), f eq i ≡ f eq i (x, t) and S i ≡ S i (x, t), whereas terms evaluated at position x + c i δt and time t + δt are noted with a star: f i ≡ f i (x + c i δt, t + δt), f eq i ≡ f eq i (x + c i δt, t + δt) and S i ≡ S i (x + c i δt, t + δt). With those notations, integration of Eq. (C.1) over t and t + δt yields:

f i = f i - δt 2τ f i -f eq i - δt 2τ f i -f eq i + δt 2 S i + δt 2 S i (C.2)
where the trapezoidal rule was applied for the right-hand side of Eq. (C.1). In this expression, the natural variable change for implicit terms is

f i = f i + δt 2τ f i -f eq i - δt 2 S i . (C.3)
The same variable change is used for f i :

f i = f i + δt 2τ f i -f eq i - δt 2
S i .

(C.4)

By inverting the latter relation in order to express f i with respect to f i , we obtain:

f i = 2τ 2τ + δt f i + δt 2τ f
f i = f i - δt τ + δt/2 f i -f eq i + δtS i . (C.8)
Without using the previous variable change for f eq i , Eq. (C.6) is equivalent to

f i = f i - δt τ + δt/2 f i -f eq i + τδt τ + δt/2 S i , (C.9)
where only the factor in front of the source term is modified.

By introducing the dimensionless collision rate which is defined by τ = τ/δt, Eq. (C.8) finally writes

f i = f i - 1 τ + 1/2 f i -f eq i + δtS i , (C.10)
or alternatively, The purpose of this Appendix is to prove the equivalence between the source term and the modification of the equilibrium distribution function. The lattice Boltzmann scheme for the conservative Allen-Cahn equation is (Eq. (29a) with S g i defined by Eq. (30a)):

f i = f i - 1 τ + 1/2 f i -f eq i + τδt τ + 1/2 S i .
g i = g i - 1
τ g + 1/2 g ig eq i + P g i δt + F g i δt (C. [START_REF] Gunstensen | Lattice Boltzmann model of immiscible fluids[END_REF] with the mobility coefficient defined by M φ = τ g c 2 s δt. By using the definition of g eq i for g eq i = φΓ i -δtP g i /2 -F g i δt/2 and gathering the term F g i δt inside the bracket, we obtain g i = g i -1

τ g + 1/2 g i -φΓ i -F g i τ g δt + δt 2

P g i + P g i δt (C.14)
Next, the collision rate is replaced by its mobility τ g = M φ /(c 2 s δt):

g i = g i - 1 τ g + 1/2 g i -φΓ i -F g i M φ c 2 s + δt 2 P g i + P g i δt (C.15)
Finally, if we use the definition of F g i given by Eq. (30b), the Allen-Cahn equilibrium distribution function g eq,CAC i can be defined by [START_REF] Fakhari | A weighted multiple-relaxation-time lattice Boltzmann method for multiphase flows and its application to partial coalescence cascades[END_REF] g eq,CAC i

= φΓ i + M φ 4 W φ(1 -φ)w i c i • n c 2 s (C.16)
and the alternative lattice Boltzmann equation is

g i = g i - 1
τ g + 1/2 g ig eq,CAC i

+ P g i δt (C.17) with g eq,CAC i = g eq,CAC i

-P g i δt/2 with P g i defined by Eq. (30b).

Bibliography

  e 7, 8 = (±1, 1, 0) T , e 9, 10 = (±1, -1, -0) T , e 11, 12 = (±1, 0, 1), e 13, 14 = (±1, 0, -1) T , e 15, 16 = (0, ±1, 1) T , e 17, 18 = (0, ±1, -1) T . For D3Q7 N pop = 6, e 2 = 1/4, w 0 = 1/4 and w 1,...,6 = 1/8. For D3Q15 N pop = 14, e 2 = 1/3, w 0 = 2/9, w 1,...,6 = 1/9 and w 7,...,14 = 1/72. For D3Q19 N pop = 18, e 2 = 1/3, w 0 = 1/3, w 1,...,6 = 1/18 and w 7,...,18 = 1/36. The standard notations will be used: c i = e i c with c = δx/δt where δx and δt are the space-and time-steps respectively and c 2 s = c 2 /3.

Figure 1 :

 1 Figure 1: 3D lattices of LB scheme.

Figure 2 :

 2 Figure 2: Example of using the Kokkos library to compute the zeroth-order moment of distribution function.

  (a) Comparison of computational times for three NVIDIA graphical cards: K80 (oldest), P100 and V100 (newest). (b) Comparisons of computational times for three optimizations of LBM kernel for Intel KNL: fused (left), CSoA (middle) and CSoA2 (right).

Figure 3 :

 3 Figure 3: Computational times (in Million Lattice Updates Per Second -MLUPS) for a diffusive problem with a D3Q19 lattice. (a) GPU for three mesh sizes and (b) on CPU (Intel KNL) for two mesh sizes.

Figure 4 :

 4 Figure 4: Zalesak's slotted disk with the conservative Allen-Cahn model.

  a) Without change of rotation during simulation. Streamlines of u (colored lines) and interface φ = 0.5 (black line) at three times.

  b) With a change of rotation direction during simulation. Streamlines of u (colored lines) and contours φ = 0.5 (black lines) for four times.

Figure 5 :

 5 Figure 5: Deformation of an initial disk standing inside a vortex. (a) Without change of rotation during simulation. (b) With change of rotation.

  ν A /ν B =1/3 LBM Analytical ν A /ν B =1/5 LBM Analytical ν A /ν B =1/10 LBM (a) Comparisons between LBM and the double-Poiseuille analytical solution for three viscosity ratios and ρ A /ρ B = 1. ρ A /ρ B =1/1.658 LBM Analytical ρ A /ρ B =1/2 LBM Analytical ρ A /ρ B =1/3 LBM (b) Comparisons between LBM and the double-Poiseuille analytical solution for three density ratios and ν A /ν B = 1. .u.) Analytical ρ A /ρ B =1 LBM linear interp LBM harmonic interp Analytic ρ A /ρ B =1/1.658 LBM harmonic interp (c) Red line -comparisons between two interpolation methods of viscosity: linear (black squares) and harmonic mean (red circles) defined by Eqs. (38) and (1b) respectively. Blue line -verification for ρ A /ρ B = 1/1.658. sol σ = 0.15 Analytical sol σ = 0.08 Analytical sol σ = 0.04 LBM σ = 0.15 LBM σ = 0.08 LBM σ = 0.04 (d) Laplace's law verification for three values of surface tension.

Figure 6 :

 6 Figure 6: Verification of coupling between the phase-field equation and fluid flow model without phase change. (a) Double-Poiseuille flow with three viscosity ratios. (b) Double-Poiseuille for three density ratios. (c) Effect of linear interpolation and harmonic mean of viscosity. (d) Laplace law.

  Comparisons for three values of thermal diffusivity: α 1 g = 0.14 (red), α 2 g = 0.125 (blue) and α 3 g = 0.08 (black). Left: x-profiles of temperature θ at the end of simulation t f = 2 × 10 5 . Right: evolution of interface position x I (t) tracked by φ = 1/2. The temperature interface is θ I = 0. , 5, 2 for j = 1, 2, 3 with α 1 l = 0.14 (red), α 2 l = 0.125 (blue) and α 3 l = 0.08 (black). Left: x-profiles of temperature θ at the end of simulation t f = 2 × 10 5 . Right: evolution of interface position x I (t) tracked by φ = 1/2. The temperature interface is θ I = 0.05.

Figure 7 :

 7 Figure 7: Comparisons between LBM (dots) and analytical solution of Stefan problem (solid lines). (a) With α l /α g = 1 and θ I = 0. (b) With α l /α g = 1 and θ I = 0.

  Detachment of bubbles occurs on nodes at the early stage of the simulation.

  Coalescence is observed at the top of the domain for bubbles detached from nodes. It is also observed a detachment of bubbles at anti-nodes.

  Later during the simulation, bubbles are detached on nodes, the cycle is pursued periodically. (d) Streamlines (white lines) and interface φ = 1/2 (black lines) superimposed on the velocity magnitude (colored field) at t 158.03.

Figure 8 :

 8 Figure 8: Simulation of film boiling for Ja = 0.025. Interface position φ = 1/2 superimposed on temperature field and for several dimensionless times of simulation. Three successive times from (a) t a 23.94, (b) t b 114.93 and (c) t c 148.45 with δt = 4.79.

Figure 9 :

 9 Figure 9: Velocity magnitude (colored field) and interface position φ = 1/2 (black lines) at t = 95.78 for (a) Ja = 0.025 and (b) Ja = 0.1.

7 )

 7 Finally by setting a = 1/2 the conservative Allen-Cahn equation with a source term is∂φ ∂t +∇ ∇ ∇ • (uφ) = ∇ ∇ ∇ • M φ ∇ ∇ ∇φ -

  W 2 φ(1φ)(1 -2φ) . (A.10)AppendixB. Numerical value of coefficient AWhen the matched asymptotic expansions are carried out on the one-dimensional phase-field model, the coefficient A is defined by four integrals I , J , G and U by (e.g. [49, Eq. (59)]):ζ φ 0 ) 2 , J = ˆ∞ -∞ dζ(∂ ζ φ 0 )p 0 φ , G = ˆ∞ -∞ dζ(∂ ζ φ 0 )p 0 φ ˆζ 0 dζ h 0 , and U = ˆ0 -∞ dζh 0 . (B.2)In Eq. (B.2), the functions φ 0 , p 0 φ and h 0 of our model are defined such as = φ 0 (1φ 0 ), and h 0 = φ 0 (B.3)

6 ,

 6

(C. 11 )AppendixC. 2 .

 112 In Section 3, Eq. (C.10) is the starting point for each lattice Boltzmann equation. The variable change Eq. (C.4) leads to the calculation of the zeroth-order moment: Equivalence of lattice Boltzmann formulations for the Allen-Cahn equation

Table 1 :

 1 List of architectures that are currently compatible with the Kokkos library.

	Intel CPUs	NVidia GPUs ARM	IBM	AMD
	Sandy/Ivy Bridge	Kepler	ThunderX Blue gene Q AMD CPUs
	Haswell	Maxwell	ARMv8.0 Power7
	Skylake	Pascal	ARMv8.1 Power8
	Westmere CPUs	Volta		Power9
	Knights Landing/Corner Xeon Phi Turing
	Broadwell Xeon E-class

Acknowledgments

We would like to thank MATHIS PLAPP for the insightful discussions on theoretical aspects of phase-field models.

Liquid and gas properties

Interface properties

Other parameters

The thermal-hydrodynamics of this problem is controlled by several dimensionless numbers: the Grashof number Gr = ρ g g y (ρ lρ g )λ 3 s /ρ 2 g ν 2 g , the Prandtl number Pr = ν g /α g and the Jacob number Ja = C p (T w -T sat )/L. Moreover the solution is sensitive to parameters that are involved in Eq. [START_REF] Kobayashi | Modeling and numerical simulations of dendritic crystal growth[END_REF]. Several sensitivity simulations on parameters of the initial condition can be found in [START_REF] Singh | Numerical investigation of film boiling on a horizontal wavy wall[END_REF].

Simulations of film boiling with LBM saclay are first carried out inside a two-dimensional domain Ω = [0, 1.28] 2 which is discretized with N x × N y = 1024 × 1024 nodes. The space-and time-steps are respectively equal to δx = 1.25 × 10 -3 and δt = 7.5 × 10 -5 . The D2Q9 lattice is used for all distribution functions f i , g i and s i . For parameters of Table 2, the value of critical wavelength is λ c = 2πλ s = 0.2738, with λ s = 4.358 × 10 -2 . The Jacob number is Ja = 0.025, the Prandtl Pr = 0.2 and the Grashof number is Gr = 871.38.

Simulation of bubble detachment on nodes and antinodes

We present one simulation for which the interface is initialized by Eq. ( 42) with y 0 = 0.03, y 1 = 0.015 and λ = 0.64. The choice λ = 0.64 was done after one first preliminary simulation which was performed with λ = 0.32 (> λ c = 0.2738) to check detachment of bubbles. For λ = 0.64, the maximum value of y is y max = 0.045 for two positions x 2) ymax are called "nodes" and x

(1), [START_REF] Guo | Lattice Boltzmann Method and Its Applications in Engineering[END_REF] ymin are called "anti-nodes". Here, we present one simulation to observe detachment of bubbles alternatively on nodes and anti-nodes. Actually, it is what we observe on Figs. 8a-8c which present the temperature fields and the iso-values φ = 1/2 (black line) at several dimensionless times. The dimensionless time is defined by t = t/t s where t s = λ s /g y = 0.1044. At the early stage of simulation (Fig. 8a), we can observe that the detachment of bubbles occurs on nodes. Later during the simulation (Fig. 8b), the bubbles that are emitted on nodes coalesce on the top on the domain, while two other bubbles grow and are detached from anti-nodes. Finally (Fig. 8c), the cycle is repeated periodically: bubbles emitted at anti-nodes coalesce and new bubbles on nodes detach and rise. Streamlines and velocity magnitude corresponding to the last time t 158.03 are presented on Fig. 8d. In Table 2, the mobility was set to M φ = 1.7 × 10 -3 after a sensitivity analysis. If M φ is too low, the authors have observed the appearance of parasitic bubbles in the liquid phase. The mobility coefficient is directly related to the relaxation time τ g and the algorithm can be unstable if its value is too low. It is expected that a wider range of parameter M φ could be reached with the TRT or MRT collisions operators.

Computational times

A first comparison of computational times between GPUs and CPUs has been indicated in Section 3.5, but only for a diffusive problem. For a single-and double-Poiseuille flow of Section 4.1, the computational times on a 100 × 100 lattice are respectively 56 MLUPS and 38 MLUPS. Those computations have been performed on a computer equipped of one AMD CPU processor (Ryzen 5 2600, 3.4GHz with 12 threads). The MLUPS are higher for the single-phase because the algorithm requires much less floating points computations. There is neither Allen-Cahn equation nor intermediate gradient to update for a single-phase flow. However, let us note that the MLUPS for two distribution functions are higher than half of the value obtained with only one (i.e. 28 MLUPS), which indicates a good code optimization by resolving the phase-field equation.

Simulation of diffusion or single-phase flow requires only one distribution and the double-Poiseuille flow requires two distribution functions. The film boiling simulation requires three Lattice Boltzmann equations with three distribution functions and the computation of additional gradients. In that case, to complete 5.33 × 10 5 time iterations on a computational domain of 1024 2 nodes, the simulation took 1h56m (80.96 MLUPS) on a single GPU NVIDIA K80. The same simulation took 12h57m (11.97 MLUPS) on 16-cores Intel Xeon CPU E5-2630 v3 2.40GHz. The computation on GPU is quicker than on CPU as expected after the preliminary diffusion simulation of section 3.5. The ratio is 6.7 times in favor of GPU compared to CPU. Next, the full grid (1024 2 nodes) is decomposed in four sub-domains composed of 256 × 1024 nodes, each one of them being taken in charge by one GPU. The simulation took 38 minutes (249.99 MLUPS) to perform the same number of time iterations on four parallel GPUs. The computational time is divided by a factor three compared to a single GPU. Finally, the computational domain is increased to Ω = [0, 5.12] × [0, 3.84] and discretized by N x × N y = 4096 × 3072 nodes, i.e. the mesh size is twelve times bigger than the previous one. The initial condition is slightly modified to AppendixA. Removal of the driven-curvature interface motion in Eq. ( 5)

In this Appendix, the derivation of first term in the right-hand side of Eq. ( 5) is reminded. The advection of phase index φ writes

If the total velocity V is defined as the sum of an external advective velocity u of an incompressible fluid plus a normal velocity of the interface

For the second term, we have used the definition of normal vector n = ∇ ∇ ∇φ/ ∇ ∇ ∇φ . If the normal velocity v n is also assumed to be separated into one term, -M φ κ, depending on the curvature κ and another one, ṽ independent on κ then: v n ∇ ∇ ∇φ = -M φ κ ∇ ∇ ∇φ + ṽ ∇ ∇ ∇φ and Eq. (A.1) writes:

For solidification problems, ṽ is the coupling with temperature equation and ensures that the Gibbs-Thomson condition is well recovered. A discussion on ṽ is presented at the end of this appendix. The next stage of the derivation is to cancel the driven-curvature interface motion M φ κ ∇ ∇ ∇φ , without setting M φ = 0, but by adding a supplementary counter term: M φ κ ∇ ∇ ∇φ -M φ κ ∇ ∇ ∇φ = S(φ). The purpose is to transform an hyperbolic-type PDE into a parabolic-type PDE by expanding κ in the first term with its definition κ = ∇ ∇ ∇ • n = ∇ ∇ ∇ • (∇ ∇ ∇φ/ ∇ ∇ ∇φ ) in order to obtain an expression involving the laplacian of φ:

The main advantage of this formulation (Eq. (A.3)) is that, for a plane interface, i.e. κ = 0, the equilibrium solution of S(φ) = 0 is an hyperbolic tangent. By using the definition of n, Eq. (A.3) becomes S(φ where ζ is the normal coordinate of the interface, a controls the slope of the hyperbolic tangent and W is the interface width. The above kernel function ensures an hyperbolic tangent profile at equilibrium. It is consistent with the profile obtained in