Georges Sadaka
email: georges.sadaka@univ-rouen.fr

Aina Rakotondrandisa
email: aina.rakotondrandisa@etu.univ-rouen.fr

Pierre-Henri Tournier
email: tournier@ljll.math.upmc.fr

Francky Luddens
email: francky.luddens@univ-rouen.fr

Corentin Lothodé
email: corentin.lothode@univ-rouen.fr

Ionut Danaila
email: ionut.danaila@univ-rouen.fr

Parallel finite-element codes for the simulation of two-dimensional and three-dimensional solid-liquid phase-change systems with natural convection

Keywords: phase change, PCM, parallel computing, finite element, 3D melting, 3D solidification, mesh adaptivity, Navier-Stokes-Boussinesq, FreeFem Mac, Super-computer, Linux

come

Introduction

Accurate and efficient numerical simulations of solid-liquid phase-change problems are needed in many practical applications. Metal casting, Earth's mantle formation and food freezing are well explored topics in this area. Recently, a great deal of attention was paid to the simulation of latent thermal energy storage (LTES) systems based on phase-change materials (PCM). Such devices are used for thermal energy storage (e. g. for solar power generation) or passive temperature control (e. g. for modern portable electronics) devices. For a review of various applications of PCMs, see recent reviews by [START_REF] Agyenim | A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (lhtess)[END_REF] and [START_REF] Kalnaes | Phase change materials and products for building applications: a state-of-theart review and future research opportunities[END_REF].

In all these problems, melting and solidification are fundamental processes that are difficult to simulate if accurate physical models are used. A sketch of the melting problem is shown in Fig. 1a. Buoyancy forces in the liquid (melted) phase generate a significant convective flow and thus deform the liquid-solid interface. Solidification is a similar process, but with a slower evolution and the possibility to generate several solid-liquid interfaces propagating simultaneously. Resolving all the scales in the liquid region using the Navier-Stokes-Boussinesq equations and accurately capturing the solid-liquid interfaces are the main challenges for a numerical system addressing these problems. In our recent contribution [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF] we presented an in-depth review of physical and numerical models dealing with solid-liquid phase-change problems with convection. The approach retained in [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF] was based on the widely used enthalpy-porosity single-domain model, also called the fixed-domain model [START_REF] Brent | Enthalpy-porosity technique for modeling convectiondiffusion phase change: Application to the melting of a pure metal[END_REF]. In the present contribution, we extend this model to parallel computations using domain-decomposition (DD) methods, as illustrated in Fig. 1b.

The main advantage of the enthalpy-porosity single-domain model is that the solidliquid interface is not explicitly tracked. The position of the interface is computed a posteriori by post-processing the obtained temperature field T (in Fig. 1b, the interface is represented by the iso-contour T = 0). This makes the model appropriate for the use of domain decomposition models and parallel computing. Front tracking or front fixing deforming-grid methods (e. g. [START_REF] Stella | Modelling methodologies for convection-diffusion phase-change problems[END_REF]; [START_REF] Tenchev | Finite element moving mesh analysis of phase change problems with natural convection[END_REF]) would obviously introduce an algorithmic complexity penalizing parallel computing performances. Note that a different category of models was recently suggested in the literature, based on the Lattice Boltzmann Method [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF][START_REF] Gong | Numerical simulation of melting with natural convection based on lattice Boltzmann method and performed with cuda enabled gpu[END_REF]. Such methods based on non-deterministic models are also well adapted for parallel computations. For a comprehensive review of models for phase-change problems with convection, see [START_REF] Kowalewski | Phase change with convection: modelling and validation[END_REF].

Another advantage of the enthalpy-porosity single-domain model is that the same Navier-Stokes-Boussinesq (NSB) system of equations is solved in both liquid and solid phases. Keeping in mind that the NSB equations are pertinent only in the liquid phase, their extension into the solid phase has to preserve the single-domain formulation. This is generally obtained by introducing in the momentum equations a penalization source term that brings the velocity to zero in the solid and do not affect the liquid. One of the most used expression of the penalization source term follows the Carman-Kozeny model for the permeability of a porous medium [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF][START_REF] Hannoun | A reference solution for phase change with convection[END_REF]Belhamadia et al., 2012), but other mathematically equivalent expressions were suggested [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Favier | Rayleigh-Bénard convection with a melting boundary[END_REF]. The energy equation is also modified to be valid in both phases using an enthalpy (temperature-transform) model introducing a regularized latent-heat term. Different formulations and implementations of the enthalpy-porosity model are presented in [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; Giangi et al. (2000); [START_REF] Stella | Modelling methodologies for convection-diffusion phase-change problems[END_REF].

The enthalpy-porosity single-domain model was implemented in [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF] using an adaptive finite-element method. The corresponding Toolbox for FreeFem++(a free software under LGPL license)1 was distributed with that paper. Extensive validation tests proved the ability of the numerical system to deal with melting and/or solidification problems of increasing difficulty. Problems with complex shapes of the computational domain or with multiple solid-liquid interfaces were successfully computed. The mesh adaptivity capabilities of the method ensured reasonable computational costs, since the mesh was refined only in the zones of large gradients of variables (artificial mushy region, boundary layers) and coarsened in the solid zones with low gradients. However, some cases, such as the water freezing or the melting-solidification cycle of a PCM, demanded hours or days on a personal computer to simulate physically pertinent time evolution.

Consequently, the main purpose of the present contribution is to extend the singledomain model of [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF] to parallel computing using domaindecomposition methods. For two-dimensional (2D) configurations, the new Toolbox can be used to reduce the computational time on personal computers with multi-core processors. Using high-performance computing (HPC) facilities, the Toolbox is well adapted to simulate 3D configuration of phase-change problems. It is important to note that very few 3D simulations with accurate capturing (mesh adaptivity) of the solid-liquid interface were reported in this research area. Adaptive FE methods were suggested for classical two-phase Stefan problem in 2D and 3D (Belhamadia et al., 2004a,b). For phasechange systems with convection, adaptivity strategies were suggested and tested only for 2D problems [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF][START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF][START_REF] Belhamadia | A two-dimensional adaptive remeshing method for solving melting and solidification problems with convection[END_REF]. An attempt to adapt the mesh in 3D simulations for melting phenomena was undertaken by [START_REF] Zimmerman | Monolithic simulation of convection-coupled phase-change: Verification and reproducibility[END_REF] using an AMR (Adaptive Mesh Refinement) technique based on a dual-weighted residual method. The method displayed a major drawback, since the mesh was refined during the advancement of the liquid-solid front, but never coarsened behind. As a consequence, only a very preliminary coarse grid simulation of the 3D melting was possible using the AMR method. We show in this contribution that our new method using domain decomposition and mesh adaptivity is perfectly suited to the accurately simulate 3D melting or solidification, in simple or complex-shape geometries.

The new parallel Toolbox is based on a single-domain enthalpy-porosity model for solid-liquid phase change problems with convection. For the energy conservation equation, a temperature-based formulation takes into account the latent heat by introducing a discontinuous source term. For the mass and momentum conservation equations, we solve in the entire domain the incompressible Navier-Stokes equations with Boussinesq approximation for buoyancy effects. To bring the velocity to zero in the solid phase, we introduce in the momentum equation a penalty term following the Carman-Kozeny model. The coupled system of momentum and energy equations is integrated in time using a second-order implicit Gear scheme. The resulting discretized equations are solved using a Newton method [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. For the space discretization we use Taylor-Hood triangular finite elements, i. e. P 2 for the velocity and temperature and P 1 for the pressure. Temperature is discretized using P 2 or P 1 finite elements. Discontinuous variables (latent heat, thermal diffusivity, etc) at the solid-liquid interface are regularized through an intermediate artificial mushy region.

To enable parallel computing, we use the ffddm framework [START_REF] Tournier | ffddm online documentation[END_REF], which is a set of FreeFem++ scripts implementing Schwarz domain-decomposition methods for the efficient solution of linear systems.The mesh is first split into subdomains using an automatic graph partitioning library, such as Scotch [START_REF] Pellegrini | SCOTCH: A Software Package for Static Mapping by Dual Recursive Bipartitioning of Process and Architecture Graphs[END_REF] or Metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF]. Each subdomain is assigned to a MPI process. The linear systems are assembled and solved in parallel using a GMRES Krylov method with a Restricted Additive Schwarz (RAS) preconditioner [START_REF] Dolean | An Introduction to Domain Decomposition Methods: algorithms, theory and parallel implementation[END_REF]. Mesh adaptivity using metrics control makes possible the optimization of the distribution of mesh elements. For 3D cases, the mesh adaptivity is more involved: the metric is first computed using mshmet, which is a module inside FreeFem++, and then used within the mmg2 remeshing tool [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] to generate the adapted new mesh.

The paper is organized as follows. Section 2 introduces the enthalpy-porosity single domain model based on the Navier-Stokes-Boussinesq equations. Section 3 presents the adaptive finite-element numerical method using a Newton algorithm for the non-linear discretized equations. We also discuss in this section the Schwarz domain decomposition method used for parallel computations and the algorithm for mesh adaptivity (with emphasis on the 3D adaptivity technique). A description of the programs contained in the provided 2D and 3D toolboxes is given in Section 4. The next two sections are devoted to extensive numerical validations of the method for 2D benchmarks ( §5) and corresponding 3D configurations ( §6). The robustness of the algorithm is demonstrated by comparing our results with reference data available in the literature. The capabilities of the Toolbox to deal with complex geometries are also illustrated. The main features of the software and possible extensions are summarized in Section 7.

Navier-Stokes-Boussinesq equations and enthalpy-porosity model

We consider a solid-liquid system placed in a three-dimensional domain Ω. The dimensionless system of equations to be solved in both liquid and solid regions is based on the incompressible Navier-Stokes equations, with Boussinesq approximation for buoyancy effects, and a temperature transforming model for the energy equation [START_REF] Voller | An enthalpy method for convection/diffusion phase change[END_REF][START_REF] Cao | A numerical analysis of Stefan problems for generalized multidimensional phase-change structures using the enthalpy transforming model[END_REF]:

∇ • u = 0, (1) ∂u ∂t + (u • ∇)u + ∇p - 1 Re ∇ 2 u -f B (θ) e z -A mushy (θ)u = 0, (2)
∂ (Cθ) ∂t + ∇ • (Cθu) -∇ • K ReP r ∇θ + ∂ (CS) ∂t = 0. (3)
Non-dimensional space, velocity, temperature and time variables in (1)-(3) were obtained from physical ones after applying the following scaling:

x → x H , u → u V ref , θ = T -T ref δT , t → V ref H t, ( 4
)
where H is the reference length (usually the height of the cavity when a rectangular domain is considered) and V ref a reference velocity that will be defined differently for melting and solidification problems. T ref is the reference temperature and in most cases T ref = T f (the temperature of fusion), unless otherwise specified. Consequently, the non-dimensional temperature of fusion is set to θ f = 0. Temperature difference δT defines a temperature scale, that will be set differently for melting and solidification cases.

The linearized Boussinesq buoyancy force (f B ), the Reynolds (Re) and Prandtl (P r) numbers are defined as (subscripts s and l refer to the solid and the liquid phases, respectively):

f B (θ) = Ra P rRe 2 θ, Re = ρV ref H µ l = V ref H ν l , P r = ν l α l , (5)
with ν the kinematic viscosity and α = k/(ρc) the thermal diffusivity. In the expression of f B , the Rayleigh number of the flow is defined as:

Ra = gβH 3 δT ν l α l , (6)
with β the thermal expansion coefficient and g the gravitational acceleration.

If previous non-dimensional numbers are pertinent only in the liquid phase, the non-dimensional conductivity K and specific heat C are defined in both phases:

K(θ) = k k l = 1, θ ≥ θ f , k s /k l , θ < θ f . , C(θ) = c c l = 1, θ ≥ θ f , c s /c l , θ < θ f . ( 7
)
The non-dimensional function S = s/s l in the energy equation ( 3) takes a similar nondimensional form:

S(θ) = s s l =    h sl /c l δT = 1 Ste , θ ≥ θ f , 0, θ < θ f , (8)
with h sl the latent heat of fusion and Ste the Stefan number. Discontinuous step-functions defined in ( 7) and ( 8) are replaced by continuous and differentiable hyperbolic-tangent functions, generically defined for all θ by the formula [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]:

F (θ; a s , θ s , R s ) = f l + f s -f l 2 1 + tanh a s θ s -θ R s , ( 9
)
where f l , f s are the imposed values in the liquid and solid phases, a s a smoothing parameter, θ s the central value (around which we regularize) and R s the smoothing radius. For example, we use for the non-dimensional source term in (3) the following regularization over the artificial mushy region θ ∈ [-ε, ε]:

S(θ) = 1 Ste - 1 2 Ste 1 + tanh θ r -θ R s , ( 10
)
where θ r is the central value around which we regularize (typically θ r = θ f = 0) and R s the smoothing radius (typically R s = ε).

Finally, the penalty term A mushy (θ)u in momentum equation (2) follows from the Carman-Kozeny model [START_REF] Voller | An enthalpy method for convection/diffusion phase change[END_REF][START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF][START_REF] Kheirabadi | The effect of the mushy-zone constant on simulated phase change heat transfer[END_REF]:

A mushy (θ) = - C CK (1 -L f (θ)) 2 L f (θ) 3 + b , ( 11
)
where L f (θ) is the local liquid fraction, which is 1 in the fluid region and 0 in the solid. L f is regularized inside the artificial mushy-region using a hyperbolic-tangent similar to (10). The Carman-Kozeny constant C CK is set to a large value (10 6 ) and the constant b = 10 -6 is introduced to avoid divisions by zero.

Numerical method

Finite-element formulation

Finite-element methods for solving Navier-Stokes type systems of equations like (1)-(3) are generally based on a separate discretization of the temporal derivative (using finite differences, splitting or characteristics methods) and the generalization of the Stokes problem for the resulting system [START_REF] Temam | Navier-Stokes equations and nonlinear functional analysis[END_REF][START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]. To simplify the presentation, we consider in the following that C = 1. For the phase-change problems considered in this paper, this is a physically valid assumption.

For the time integration, we use a second-order Gear (BDF2) finite-difference scheme (see also [START_REF] Belhamadia | An enhanced mathematical model for phase change problems with natural convection[END_REF]):

dφ dt 3φ n+1 -4φ n + φ n-1 2δt , ( 12
)
computing the solution φ n+1 at time t n+1 = (n + 1)δt by using two previous states (φ n , φ n-1 ). We use this scheme to advance in time both velocity (φ = u) and temperature fields (φ = θ). The other terms in equations ( 1)-( 3) are treated implicitly. We obtain the following implicit semi-discretization in time of the single-domain model ( 1)-( 3):

∇ • u n+1 = 0, (13) 3 2 u n+1 δt + (u n+1 • ∇)u n+1 + ∇p n+1 - 1 Re ∇ 2 u n+1 -A mushy (θ n+1 )u n+1 -f B (θ n+1 ) e z = 2 u n δt - u n-1 2δt , ( 14
) 3 2 θ n+1 + S(θ n+1 ) δt + ∇ • u n+1 θ n+1 -∇ • K(θ n+1 ) ReP r ∇θ n+1 = 2 θ n + S(θ n ) δt - θ n-1 + S(θ n-1 ) 2δt . ( 15
)
To solve the system of equations ( 13)-( 15) we use a classical Galerkin finite-element method. We consider homogeneous Dirichlet boundary conditions for the velocity, i. e. u = 0 on ∂Ω, and set the classical Hilbert spaces for the velocity and pressure:

V = V × V × V , V = H 1 0 (Ω), Q = q ∈ L 2 (Ω) Ω q = 0 (16)
Following the generalization of the Stokes problem [START_REF] Temam | Navier-Stokes equations and nonlinear functional analysis[END_REF][START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF], the weak formulation of the system (13)-( 15) can be written as: find

(u n+1 , p n+1 , θ n+1 ) ∈ V × Q × V , such that: b u n+1 , q -γ(p n+1 , q) = 0, ∀ q ∈ Q (17) 3 2δt u n+1 , v + c u n+1 ; u n+1 , v + 1 Re a u n+1 , v -(A mushy (θ n+1 ) u n+1 , v) + b v, p n+1 -f B (θ n+1 ) e z , v = 2 δt (u n , v) - 1 2δt u n-1 , v , ∀ v ∈ V (18) 3 2δt θ n+1 + S(θ n+1 ), φ + u n+1 • ∇θ n+1 , φ + K(θ n+1 ) ReP r ∇θ n+1 , ∇φ = 2 δt (θ n + S(θ n ), φ) - 1 2δt θ n-1 + S(θ n-1 ), φ , ∀ φ ∈ V , ( 19
)
where (u, v) = Ω u • v denotes the scalar product in L 2 (Ω) or L 2 (Ω) 2 ; the bilinear forms a, b and trilinear form c are defined as [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Quarteroni | Numerical Approximation of Partial Differential Equations[END_REF]:

a : V × V → R, a(u, v) = Ω ∇ t u : ∇v = 3 i,j=1 Ω ∂jui • ∂jvi, b : V × Q → R, b(u, q) = - Ω ∇ • u q = - 3 i=1 Ω ∂iui • q, c : V × V × V → R, c(w; z, v) = Ω [(w • ∇) z] • v = 3 i,j=1 Ω wj(∂jzi)vi.
Note that we introduced in Eq. ( 17) a penalty term on the pressure. An extensive discussion of the role of this parameter in reinforcing the incompressibility constraint and in the stabilization of the method is provided in [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF]. We recall the main lines of this discussion. Following the incompressibility constraint (13), the pressure is defined up to an additive constant and imposing that p n+1 ∈ Q removes this uncertainty. Equation ( 17) numerically ensures that the pressure field is of zero average in Ω. The penalty term in (17) acts at an algebraic level, by modifying the matrix of the final system (a zero lower diagonal block is avoided). Since in the present toolbox we use only iterative solvers (GMRES), this modification is not essential. It is in exchange important when LU-type direct solvers without pivoting are used, essentially for 2D calculations (UMFPACK solver in [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF]; [START_REF] Aldbaissy | A full discretisation of the time-dependent Boussinesq (buoyancy) model with nonlinear viscosity[END_REF], SuperLU in [START_REF] Woodfield | Stability and finite element approximation of phase change models for natural convection in porous media[END_REF]). Finally, the penalty constant γ has no role in stabilizing the method. Since the Taylor-Hood finite elements used in our simulations satisfy the inf-sup condition, this technique is assimilated to a stable penalty, or, in other words, just a computational trick to obtain a good solution [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF][START_REF] Boffi | Mixed finite element methods and applications[END_REF]. Note that, compared to classical penalty methods, we use in our calculations very low values of the penalty parameter (γ = 10 -7 ). This ensures very low values of the average on Ω for the pressure and also for the divergence of the velocity field.

The system of non-linear equations ( 17)-( 19) is solved using a Newton method. To advance the solution from time t n to t n+1 , we start from an initial guess w 0 = (u n , p n , θ n ) (which is the solution at t n ), and construct the Newton sequence w k = (u k , p k , θ k ) by solving for each inner iteration k:

b (u k+1 , q) -γ(p k+1 , q) = 0, (20) 3 2δt (u k+1 , v) + c (u k+1 ; u k , v) + c (u k ; u k+1 , v) + 1 Re a (u k+1 , v) - dA mushy dθ (θ k ) θ k+1 u k , v -(A mushy (θ k ) u k+1 , v) + b (v, p k+1 ) - dfB dθ (θ k ) θ k+1 ez, v = 1 δt 2u n - 1 2 u n-1 , v +c (u k ; u k , v) - dA mushy dθ (θ k ) θ k u k , v - dfB dθ (θ k ) θ k -fB(θ k ) ez, v , (21) 3 2δt θ k+1 + dS dθ (θ k ) θ k+1 , φ + (u k • ∇θ k+1 , φ) + (u k+1 • ∇θ k , φ) + K(θ k ) ReP r ∇θ k+1 , ∇φ + dK dθ (θ k ) θ k+1 ReP r ∇θ k , ∇φ = 2 δt (θ n + S(θ n ), φ) + 3 2δt dS dθ (θ k ) θ k -S(θ k ), φ (u k • ∇θ k , φ) - 1 2δt θ n-1 + S(θ n-1 ), φ + dK dθ (θ k ) θ k ReP r ∇θ k , ∇φ . ( 22
)
Note that the last term of Eq. ( 21) cancels in the case of a linear Boussinesq force f B (see Eq. ( 5)); this is not the case when non-linear variations of the density of the liquid are considered (convection or solidification of water). Note also that the previous system of equations ( 20)-( 22) depends only on u n , u n-1 , θ n and θ n-1 and is independent of p n , the pressure being in this approach a Lagrange multiplier for the divergence free constraint. The Newton loop (following k) has to be iterated until convergence for each time step δt following the algorithm:

Navier-Stokes time loop following n set w 0 = (u n , p n , θ n ) Newton iterations following k solve (20) to get w k+1 stop when w k+1 -w k < ξ N actualize (u n+1 , p n+1 , θ n+1 ) = w k+1 . (23)
For the space discretization of the system (20)-( 22) we use standard Taylor-Hood finite elements [START_REF] Taylor | A numerical solution of the Navier-Stokes equations using the finite element technique[END_REF], approximating the velocity with P 2 (piecewise quadratic) finite elements (V h space), and the pressure with the P 1 (piecewise linear) finite elements (Q h space):

V h = v ∈ C 0 ( Ω) 2 ∀K ∈ T h , v |K ∈ P 2 , ( 24
)
Q h = v ∈ C 0 ( Ω) ∀K ∈ T h , v |K ∈ P 1 , ( 25
)
where K is an element of the triangulation T h , with characteristic mesh size h. Temperature and enthalpy variables are discretized using either P 1 or P 2 finite elements. This algorithm was implemented using the open-source software FreeFem++ [START_REF] Hecht | FreeFem++ (manual)[END_REF][START_REF] Hecht | New developments in Freefem++[END_REF]. The FreeFem++programming framework offers the advantage to hide all technical issues related to the finite element method. The high level programming language with syntax close to mathematical formulations, makes the implementation of the present numerical algorithm very easy. Similar algorithms based on FreeFem++were successfully used for solving different systems of equations with locally sharp variation of the solution, such as Gross-Pitaevskii equation [START_REF] Danaila | A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose-Einstein condensates[END_REF][START_REF] Vergez | A finite-element toolbox for the stationary Gross-Pitaevskii equation with rotation[END_REF] or Laplace equations with nonlinear source terms [START_REF] Zhang | Existence and numerical modelling of vortex rings with elliptic boundaries[END_REF].

The FreeFem++syntax to implement the Newton algorithm is very close to the mathematical formulation given above. We present below the main elements of syntax for the 3D formulation. We start by creating a new type Wh that will be used to define variables gathering in a vector all the unknowns of the problem: the three components of the velocity, the pressure and the temperature. Wh thus corresponds to V × V × V × Q × V (see Eq. ( 16)). Each component of a vector of type Wh is independently assigned with a finite-element discretization (P 2 , P 1 or other available in FreeFem++). In our case, we use Taylor-Hood finite elements for the fluid part and P 2 for the temperature and thus we define the vector finite-element space as fespace Wh(Th,[P2,P2,P2,P1,P2]). The unknowns of the problem are then defined by declaring: Wh [u1,u2,u3,p,T];. Note that Wh is associated to the mesh Th. If the mesh Th changes (following mesh adaptation), the definition of Wh is automatically associated to the new mesh. Corresponding test functions are defined similarly: Wh [v1,v2,v3,q,TT];. The next step in building the FreeFem++program is to define the weak (variational) formulation. The built-in function varf facilitates the implementation of the algorithm (20)-( 22) with a syntax close to the mathematical formulation. The use of the macro environment (which is a pre-processor command for a simple syntax replacement in the script) makes the reading of the programs very intuitive, when comparing each term to its mathematical expression. For example, the following macros define the necessary operators for the variational formulation (20)-( 22):

macro grad ( u ) [ dx ( u ) , dy ( u ) , dz ( u ) ] // macro Grad ( u ) [ grad ( u #1) , grad ( u #2) , grad ( u #3) ] // macro div ( u ) ( dx ( u #1) + dy ( u #2) + dz ( u #3) ) // macro ugrad (u , v ) ([ u #1 , u #2 , u #3] '* grad ( v ) ) // macro UgradV (u , v ) [ ugrad (u , v #1) , ugrad (u , v #2) , ugrad (u , v #3) ] // macro a ( Mu ,u , v ) ( Mu *( Grad ( u ) : Grad ( v ) ) ) // macro b (u , q ) ( -div ( u ) * q ) // macro c (w ,Z , v ) ( UgradV (w , Z ) '*[ v #1 , v #2 , v #3]) //
Note that in some macros we used the concatenation operator #. This means that when we use in a script, for example, Grad(uw), the preprocessor will automatically replace this part with [grad(uw1),grad(uw2),grad(uw3)] in the final script (that will be executed). Vector operators (transposition ' and contraction :) are naturally implemented in FreeFem++. The varf syntax is used to build the matrix and the rhs-vector of the final linear system. Consequently, it follows closely the mathematical formulation ( 20)-( 22). To simplify the presentation, we start by giving below the varf formulation for the simple case of the stationary convection of air (the terms coming from the discretization of the time derivative are not present): macro VarfStatNATCONV ( varfName , meshName , VhName ) varf varfName ([ uw1 ,uw2 ,uw3 ,pw ,Tw ] ,[ v1 ,v2 ,v3 ,q , TT ]) = int3d ( meshName , qforder = ord ) ( b ( uw , q ) -gamma * pw * q + c ( uw ,u , v ) + c (u , uw , v ) + IRe * a ( uw , v ) + b (v , pw ) -dfB ( T ) * Tw * v3 + ugrad (u , Tw ) * TT + ugrad ( uw , T ) * TT + grad ( Tw ) '* grad ( TT ) * IPr ) + bcu1 + bcu2 + bcu3 + bcT ; // EOM macro V arf rh sS ta tN AT CO NV ( varfName , meshName , VhName ) varf varfName ([ uw1 ,uw2 ,uw3 ,pw ,Tw ] ,[ v1 ,v2 ,v3 ,q

, TT ]) = int3d ( meshName , qforder = ord ) ( c (u ,u , v ) + ugrad (u , T ) * TT -dfB ( T ) * T * v3 + fB ( T ) * v3 ) + bcu1 + bcu2 + bcu3 + bcT ; // EOM
The syntax +bcu1+ ... is used to include the macros that implement boundary conditions. The correspondence with variables in Eqs. ( 20)-( 22) is the following: [v1,v2,v3,q,TT] to test functions and [u1,u2,u3,p,T] to previous Newton step

[uw1,uw2,uw3,pw,Tw] corresponds to unknowns w k+1 = (u1 k+1 , u2 k+1 , u3 k+1 , p k+1 , θ k+1 ),
w k = (u1 k , u2 k , u3 k , p k , θ k ).
For the time evolution formulations, we use [u1p,u2p,u3p,pp,Tp] and [u1pp,u2pp,u3pp,ppp,Tpp] for variables at time t n and t n-1 , respectively. Consequently, the full formulation for the Newton algorithm (20)-( 22) for the case of the melting of a PCM becomes: macro VarfPCM ( varfName , meshName , VhName ) varf varfName ([ uw1 ,uw2 ,uw3 ,pw ,Tw ] ,[ v1 ,v2 ,v3 ,q [ v1 , v2 , v3 ] + ugrad (u , Tw ) * TT + ugrad ( uw , T ) * TT + grad ( Tw ) '* grad ( TT ) * IPr + c1 * dS ( T ) * Tw * TT ) + bcu1 + bcu2 + bcu3 + bcT ; // EOM macro VarfrhsPCM ( varfName , meshName , VhName ) varf varfName ([ uw1 ,uw2 ,uw3 ,pw ,Tw ] ,[ v1 ,v2 ,v3 ,q

, TT ]) = int3d ( meshName , qforder = ord ) ( b ( uw , q ) -gamma * pw * q + c1 *[ uw1 , uw2 , uw3 , Tw ] '*[ v1 , v2 , v3 , TT ] + c ( uw ,u , v ) + c (u , uw , v ) + IRe * a ( uw , v ) + b (v , pw ) -dfB ( T ) * Tw * v3 -Amushy ( T ) *[ uw1 , uw2 , uw3 ] '*[ v1 , v2 , v3 ] -dAmushy ( T ) * Tw *[ u1 , u2 , u3 ] '*
, TT ]) = int3d ( meshName , qforder = ord ) ( c (u ,u , v ) + ugrad (u , T ) * TT -[ u1 , u2 , u3 ] '*[ v1 , v2 , v3 ]* dAmushy ( T ) *T -dfB ( T ) * T * v3 + fB ( T ) * v3 -c2 *[ u1p , u2p , u3p , Tp ] '*[ v1 , v2 , v3 , TT ] -c3 *[ u1pp , u2pp , u3pp , Tpp ] '*[ v1 , v2 , v3 , TT ] + c1 *( dS ( T ) * T * TT -S ( T ) * TT ) -c2 * S ( Tp ) * TT -c3 * S ( Tpp ) * TT ) + bcu1 + bcu2 + bcu3 + bcT ; // EOM
New variables and coefficients can be easily identified. For example, c1, c2, c3 are the coefficients (depending on δt) corresponding to the Gear scheme, Amushy is the Carman-Kozeny penalty term, S the enthalpy source term, etc.

New terms can be added to the variational formulation expressed in the varf structure, without affecting other parts of the program. The implementation of new models or numerical methods for this problem is greatly facilitated by this modular structure of programs.

Domain decomposition method

The main time consuming part in the algorithm (23) lies in the solution of the sequence of linear systems of the form Ax = b. To reduce the computational time, we implement an overlapping Schwarz domain decomposition method to solve these linear systems in parallel. This task becomes an easy job by using the parallel framework ffddm, which was recently made available in FreeFem++. ffddm is a set of high-level FreeFem++ macros that the user can call in the script to perform different steps needed for the parallel solution of a linear system using a domain decomposition preconditioner. We recall these steps below and give the corresponding ffddm macro calls that are used in the implementation of our toolbox. The detailed description of these macros can be found in the online ffddm documentation [START_REF] Tournier | ffddm online documentation[END_REF].

The first step in using ffddm is the partition of the global mesh T := T h into N S non-overlapping meshes {T i } 1≤i≤N S (see Fig. 1b). Standard graph partitioners available in FreeFem++, such as Scotch [START_REF] Pellegrini | SCOTCH: A Software Package for Static Mapping by Dual Recursive Bipartitioning of Process and Architecture Graphs[END_REF] or Metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF], can be used for this task. If δ is a positive integer, the overlapping decomposition {T δ i } 1≤i≤N S is defined recursively as follows: T δ i is obtained by including all mesh elements of T δ-1 i and adding recursively one layer of elements. For δ = 0, T δ i = T i and a non-overlapping decomposition is obtained. In practice we use δ = 1 which corresponds to an overlap of width 2. The mesh partitioning step is performed in macro ffddmbuildDmeshAug. Metis is the default partitioner. Scotch can be selected by simply declaring ffddmpartitioner = 2. Each MPI process i then holds the local mesh T δ i corresponding to subdomain i. The second important step is to build an appropriate preconditioner for the linear system. Several types of preconditioners are available in ffddm. For the mathematical background of domain decomposition methods and associated preconditioners, we refer to [START_REF] Dolean | An Introduction to Domain Decomposition Methods: algorithms, theory and parallel implementation[END_REF]. In the present toolbox, we use the Restricted Additive Schwarz (RAS) preconditioner that proved very efficient. We describe below the main ideas of the corresponding algorithm. Let W h = V h × Q h × V h be the global finite element space for velocity, pressure and temperature variables. The mesh decomposition induces a natural decomposition of the global space W h on T into N S local finite element spaces {W δ i } 1≤i≤N S , each of them defined on the corresponding local mesh

T δ i . Consider the restrictions {R i } 1≤i≤N S from W h to {W δ i } 1≤i≤N S ,

and a local partition of unity

{D i } 1≤i≤N S such that N S i=1 R T i D i R i = I n×n . ( 26
)
At the algebraic level, if n is the global number of unknowns and {n i } 1≤i≤N S are the numbers of unknowns for each local finite element space, then R i is a Boolean matrix of size n i × n, and D i is a diagonal matrix of size n i × n i , for all 1 ≤ i ≤ N S . Note that R T i , the transpose of R i , is a n × n i matrix that gives the extension by 0 from W δ i to W h . The construction of the local finite element spaces {W δ i } 1≤i≤N S and partition of unity matrices {D i } 1≤i≤N S are implemented in the macro ffddmbuildDfespaceAug.

Using these matrices, we define the RAS preconditioner as:

M -1 RAS = N S i=1 R T i D i A -1 i R i , (27)
with local subdomain matrices

{A i } 1≤i≤N S = {R i AR T i } 1≤i≤N S .
The preconditioner ( 27) is naturally parallel since its assembly requires the concurrent factorization of each {A i } 1≤i≤N S . In practice, this operation is performed locally on different processes in a distributed computing context, as one subdomain is assigned to each MPI process. Likewise, applying (27) to a distributed vector only requires peer-to-peer communications between neighboring subdomains, and a local forward elimination and backward substitution (see Chapter 8 of [START_REF] Dolean | An Introduction to Domain Decomposition Methods: algorithms, theory and parallel implementation[END_REF], for a more detailed description). Local matrices {A i } 1≤i≤N S are also used to perform the parallel matrix-vector product A * v that are main operations in a parallel GMRES algorithm with preconditioner M -1 RAS . The parallel assembly and factorization of local matrices {A i } 1≤i≤N S are performed in the macro ffddmsetup, that also defines the parallel matrix-vector product and preconditioner operators. Additionally, the macro ffddmbuildrhs computes the distributed right-hand side {b i } 1≤i≤N S = {R i b} 1≤i≤N S . Macros ffddmsetup and ffddmbuildrhs use the corresponding variational varf formulations of the problem defining the bilinear and linear parts, respectively. Finally, the linear system Ax = b is solved in parallel with a GMRES algorithm, called in the function fGMRES. We apply in our algorithms a left preconditioning by M -1 RAS . The error tolerance ε G in GMRES is variable and adapted to the convergence of the Newton algorithm by monitoring ε N = w k+1 -w k . We set

ε G = 10 -10 if ε N < 10 -4 , ε G = 10 -9 if 10 -4 < ε N < 10 -2 and ε G = 10 -2 if ε N > 10 -2 . The output of the fGMRES function is the distributed solution {x i } 1≤i≤N S = {R i x} 1≤i≤N S .
The macro fromVhi is then used to recover the global solution x from the distributed solution. The global solution is finally used for visualization and also in the mesh adaptation procedure.

A detailed description of the ffddm calls used for the different steps of the domain decomposition method can be found in the User's manual of ffddm [START_REF] Tournier | ffddm online documentation[END_REF]. In the entire simulation process, the parallel algorithm described above to assemble and solve the linear system is performed at each time step and at each Newton iteration. However, the first two steps concerning the construction of the overlapping mesh and the finite element space decomposition are performed only once at the beginning of each time step. These steps are reiterated only if a mesh change occurs, after calling the mesh adaptivity procedure (see next section).

Mesh adaptivity

For 2D simulations, we use the standard mesh adaptivity function (adaptmesh) offered by FreeFem++ [START_REF] Hecht | New developments in Freefem++[END_REF]. The key idea implemented in this function (see also [START_REF] Borouchaki | Anisotropic adaptive mesh generation in two dimensions for cfd[END_REF]; [START_REF] Castro-Diaz | Anisotropic grid adaptation for inviscid and viscous flows simulations[END_REF]; [START_REF] Hecht | Mesh adaptation by metric control for multi-scale phenomena and turbulence[END_REF]; [START_REF] George | Delaunay triangulation and meshing[END_REF]; [START_REF] Frey | Maillages. Hermès[END_REF]; [START_REF] Mohammadi | Applied Shape Design for Fluids[END_REF]) is to use the Delaunay algorithm to generate a new triangular mesh with edges close to the unit length in the metric M = |H| E , where |H(x)| is the Hessian of the variable χ at point x (after being made positive definite) and E the interpolation error for χ. This implies to modify the scalar product used in the automatic mesh generator to evaluate distance and volume by defining a scalar product based on the evaluation of the Hessian H of the variables of the problem. Equilateral elements are thus constructed, with an equally distributed interpolation error E over the edges of the mesh. The previous approach could be generalized for a vector variable χ = [χ 1 , χ 2 ]. After computing the metrics M 1 and M 2 for each variable, the retained metric is the intersection M = M 1 ∩ M 2 , defined such that the unit ball of M is included in the intersection of the two unit balls of metrics M 2 and M 1 (for details, see the procedure defined in [START_REF] Frey | Maillages. Hermès[END_REF]).

The use of the mesh adaptivity algorithm for 2D phase-change systems is described in detail in [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF]. In 2D problems considered below, we took into account several metrics computed from different variables monitoring the evolution of the phase-change system. For the natural convection system, the mesh was adapted using the values of the two velocity components and the temperature. For phase-change systems, to accurately track the solid-liquid interface we added the variation of the enthalpy source term in the adaptivity criterion. For water systems (convection or freezing), we also added an extra function tracking the anomalous change of density around 4 o C. To reduce the impact of the interpolation on the global accuracy for time-depending problems, we considered, for each variable used for adaptivity, the metrics computed at actual (t n+1 ) and previous (t n ) time instants (see also Belhamadia et al. (2004a)). The capabilities of the mesh adaptivity algorithm in 2D are illustrated in §5.

For 3D simulations, the metric was first computed using mshmet, which is a FreeFem++ module, following the code: load " mshmet " real [ int ] hmetric (6* ThBackup . nv ) ; hmetric = mshmet ( ThBackup ,[ u1Backup ,u2Backup ,u3Backup ] , normalization =1 , aniso =1 , nbregul =1 , hmax = hmax , hmin = hmin , err = errh ) ;

Several parameters of the mesh can be thus controlled (h min , h max , anisotropy, etc). Then, the metric was regularized in order to avoid flat tetrahedra and saved in m11[] :

load " aniso " boundaniso (6 , hmetric ,40) ; fespace Wh6Backup ( ThBackup ,[ P1 , P1 , P1 , P1 , P1 , P1 ]) ; Wh6Backup [ m11 ,m21 ,m22 ,m31 ,m32 ,m33 ]; m11 []= hmetric ;

The external software mmg [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] was finally loaded to read the saved metric and the old mesh to finally generate the adapted new mesh:

load " mmg " ThBackup = mmg3d ( ThBackup , metric = m11 [] , verbose =0 , hmin = hmin , hmax = hmax , hgrad = adaptratio , mem =10000) ;

Description of the programs

In this section, we first describe the architecture of the programs and the organization of the provided files. Then we focus on the list of input parameters and the structure of output files.

Program architecture

Figure 2 gives a schematic overview of the content of the Toolbox. To facilitate the reading and the assimilation of the programs, we separated the 2D and 3D Toolboxes in directories PCM Toolbox DDM 2D and PCM Toolbox DDM 3D. Many detailed comments were included in programs, with direct reference to the mathematical expressions used in this paper. The used FreeFem syntax was intentionally kept at a low level of technicality and supplemented with detailed comments when specific more technical syntax was used.

This directory is organized as follows:

1. The directory Common Macros contains five files:

• Macro operator.idp includes macros and functions defining mathematical operators,

• Macro problem.idp: macros defining the variational formulation of the problem,

• Macro restart.idp: macros used to start a new simulation from a saved field,

• Macro output.idp: macros used to save the solution with different formats.

2. The directory Test Cases contains four subdirectories, each of them defining one of the following applications:

• natural convection of air or water in a differentially heated square/cubic cavity,

• melting of a PCM stored in containers of different shapes,

• melting followed by solidification of a rectangular/cubic PCM,

• freezing of pure water in a square/cubic cavity. Each subdirectory contains three files: NEWTON $case.edp is the main FreeFem++script file, param phys.inc defines the physical parameters and param num.inc the numerical parameters. The 3D Toolbox is supplemented with the script Read 3D data.edp, which reads data from the RST folder and generates the files for visualization with Medit, Paraview or Tecplot. For example, to run the natural convection case of air in a square cavity, the user can use the following command in a terminal window: mpirun -np 6 FreeFem++-mpi NEWTON stat natconv ddm.edp -wg -v 0 -ns.

The folder structure of each test case is illustrated in Figure 3. The obtained solutions are saved in the folder OUTPUT/Data. Depending on the output format selected by the user, data files are generated in specific folders for being visualized with Tecplot, Paraview, Gnuplot or Medit. We also provide in the folder Figures ready-made layouts for each visualization software. The user can thus obtain the figures from this paper using newly generated data. More details about the output structure are given below.

Input parameters

Physical parameters and parameters related to the run are separated into two files.

(1) The file param phys.inc contains the physical descriptions of the problem:

• typeT: is the finite-element type for the temperature, with possible values P2 or P1,

• Torder: is the order of accuracy of the time integration scheme, with possible values 1 (Euler scheme) or 2 (Gear scheme), • scalAdim: defines the characteristic scales of the problem, see (4). Possible values 1, 2 or 3 correspond to the following choice of the characteristic scales [START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]:

(1) :

V (1) ref = ν l H =⇒ t (1) ref = H 2 ν l =⇒ Re = 1, (28)
(2) :

V (2) ref = α H =⇒ t (2) ref = t (1) ref Pr =⇒ Re = 1/Pr, (29)
(3) :

V (3) ref = ν l H Ra Pr =⇒ t (3) ref = t (1) ref Pr Ra =⇒ Re = Ra Pr , (30)
• x l , x r , y l , y r : are the values defining the dimensions of the cavity [x l , x r ] × [y l , y r ],

• Pr, Ra, Ste: are the Prandtl, Rayleigh and Stefan numbers, see ( 6) and ( 5),

• T hot , T cold : are dimensionless temperatures according to (4),

• bcu 1 , bcu 2 , bcT: are macros defining the velocity (u) and the temperature (T) boundary conditions,

• epsi: is the half width ε of the mushy region. Default value = 0.01,

• dt: is the dimensionless time step,

• t max : is the dimensionless final time,

• Parameters for regularization functions:

The parameters of the hyperbolic-tangent function (9) used to regularize discontinuous functions are set by default as follows:

f s f l a s θ s R s C CK b Enthalpy 0 1/Ste 1 0.01 0.01 - -
Carman -Kozeny 0 1 1 0.01 0.01 10 6 10 -7

Conductivity (water) 1 2.26/0.578 1 θ f 0.015 --

• rho(T) and Drho(T): (water cases only) define the density and its derivative as functions of the temperature, following the model [START_REF] Gebhart | A new density relation for pure and saline water[END_REF]:

ρ(T ) = ρ m (1 -ω|T -T m | q ), ρ m [kg/m 3 ] ω [ o C -q ] q T m [ o C] 999.972
9.2793 • 10 -6 1.894816 4.0293

• f B (T), df B (T): define the buoyancy force and its derivative.

(2) The file param num.inc contains the parameters controlling the run.

Restart parameters:

• Nsave: the solution is saved every Nsave time steps in the Data folder (see Figure 3). The temperature and velocity fields are saved in Tecplot and Medit folders, while the liquid fraction, the Nusselt number, and the accumulated heat input are saved in the Gnuplot folder.

• Nrestart: restart files (mesh and solution) are saved every Nrestart time steps. Solutions at current and previous iterations, the CPU time, the accumulated heat input Q 0 , and the time step dt are saved in the folder RST.

• Ncondt: allows the user to stop the run and save the solution properly. The file OUTPUT/zz.condt is read every Ncondt time steps: if the user replaces the value "0" in this file by "1" the run is stopped. This is a simple solution for a clean stop of the job by the user. Default value = 20.

• Nremesh: the mesh is adapted every Nremesh iterations. If this parameter is set to "1" the mesh is adapted every time step.

• IFrestart: is a Boolean controlling the set up of the initial field.

IF restart = 0, the initial condition is built in the code for each test case. For the PCM melting cases, the PCM is initially motionless at isothermal temperature. To set-up a smooth initial field, a few time steps (with very small δt) are computed by increasing progressively the boundary temperature at the hot wall and the Rayleigh number (by continuation). Outputs are saved in OUTPUT/Data-RST-0.

IF restart > 0, (positive integer values) the solution field previously computed at iteration IF restart is loaded from the folder OUTPUT/Data-RST-filenameRST/RST, with filenameRST a variable selecting the restart folder. Note that for 3D cases, the user has to copy manually the solution used as restart field to the folder INIT. IF restart < 0, (negative integer values), the same principle for loading a solution is used, but from the folder INIT (see Figure 3). The solution fields stored in this folder could come from different previous calculations (e. g. a steady state solution or, for the water, the natural convection field before freezing).

Newton parameters:

• epsconv: is the value of the stopping criterion for steady cases,

• gamma: is the penalty parameter in (13). Default value = 10 -7 ,

• tolNewton: is the Newton tolerance ξ N (see ( 23)). Default value = 10 -6 ,

• newtonMax: limits the maximum number of iterations in the Newton algorithm (23). Default value = 50.

Mesh parameters:

• nbseg: is the number of segments for the discretization along the x and y directions,

• errh: is the interpolation error level. Default value = 0.02,

• hmin, hmax: are the minimum and the maximum edge size, respectively,

• adaptratio: is the ratio for a prescribed smoothing of the metric. For a value less than 1.1 no smoothing is done. Default value = 1.5,

• nbvx: is the maximum number of vertices allowed in the mesh generator. Default value = 50000.

Output parameters:

• dircase: is the name of the output folder,

• fcase: is the prefix-name for output files.

• Paraview, Tecplot, Medit, Gnu: correspond to the name of the visualization software to be used; the format of the outputs written in OUTPUT/Data (see Figure 3) is set accordingly.

Outputs

When a computation starts, the OUTPUT directory is created (see Figure 3). It contains two folders storing the output data and the echo of the run parameters. The folder Data contains four subdirectories with different output format files of the computed solution. File names are created using the prefix defined by the parameter fcase, the current iteration and the current dimensionless time t. Solution files can be visualized using either Tecplot or any other CFD Visualization tools (Paraview, Visit, etc.). Moreover, .gmsh (mesh) and .rst (fields) files are generated in the folder RST to enable restarts of the computation. Note that the folder FFglut contains FreeFem++scripts that read and visualize the RST-files to facilitate the selection of a restart field. An .echo file with a summary of the main parameters, information on the run and the names of the output files is saved in the folder RUNPARAM. This directory additionally contains a copy of the .inc parameter files, allowing an easy identification of each case and preparing an eventual rerun of the same case.

2D parallel simulations

In this section we consider well defined 2D benchmarks used to validate numerical codes for natural convection and phase-change problems. All these cases were described in detail in [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF] and computed using sequential codes. We present below the results obtained for the same cases, but using parallel computations with the new Toolbox: (i) the natural convection of air ( §5.1), (ii) the melting or solidification of a phase-change material (PCM) ( §5.2) and (iii) the convection and the freezing of pure water ( §5.3). This approach allows us to test the programs by adding progressively non-linearities in the Newton algorithm. In the following, we validate each case with respect to the physical reference data available in the literature. To have a quick overview of the benefits of using parallel computing, we compare in Table 1 computational (CPU) times necessary to run sequential and parallel computations of each case (with default parameters) on a personal computer. We also display the maximum number of triangles of adapted mesh and the ratio between the CPU time required to adapt the mesh to the CPU time for a complete time step. Mesh adaptation is performed every time step. To make the simulations accessible and testable by using a personal computer with a multi-thread processor, we use in this section only 6 MPI processes. 20

Case

Natural convection of air

We start by testing the Newton algorithm ( 20)-( 22) for the case of natural convection, i. e. C = K = 1, A mushy (θ) = S(θ) = 0. We consider the classical problem of the thermally driven square cavity [0, 1] × [0, 1], filled with air. The Boussinesq term f B (θ) is then linear and takes the form (5). Top and bottom walls are adiabatic, while the temperature is fixed on the left (hot) wall and the right (cold) wall. Natural convection flows are computed for three values of the Rayleigh number: Ra = 10 4 , 10 5 , 10 6 . The Prandtl number is set to P r = 0.71. For these values of the Ra number, the flow is known to be steady (Le [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF] and therefore, we also solve the time-independent version of equations ( 20)-( 22).

We provide programs for both steady (time-independent) and time-dependent cases. The steady case is performed using a continuation following the Rayleigh number: a smaller value for Ra is set initially and is then smoothly increased until reaching the wanted value. The time-dependent case is computed until a steady state with a single convection cell is reached.

The temperature is imposed at the right (cold) wall as θ c = -0.5 and at the left (hot) wall as θ h = 0.5. Top and bottom walls are adiabatic. The initial condition models a cavity filled with motionless air (u = 0), with a linear distribution of the temperature. Both steady and time-dependent codes converge to the same flow state, with a single convection cell. For this final state, horizontal u(y) and vertical v(x) velocity profiles were extracted at mid-domain (y = 0.5 and x = 0.5, respectively) and plotted in Figure 4. Our results are in very good agreement with reference numerical results obtained by Le [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF] with a spectral code. Table 2 offers a quantitative assessment of the accuracy of the present Newton method. The values of u max and the location Y of this maximum are compared to reference values from Le [START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF]. The Newton method gives results identical to reference values, with a difference less than 0.1%.

Run

Melting or solidification of a phase-change material (PCM)

We continue our validation tests by considering the full system (20)-( 22) for the case of the melting or solidification of a phase-change material. Two new non-linearities are now present in the system: the Carman-Kozeny penalty term A mushy (θ) and the enthalpy source term S(θ). The function S is regularized using (10). We also consider that the material properties in the liquid and solid are the same, i. e. C = K = 1. This is a frequent assumption [START_REF] Wang | A comprehensive numerical model for melting with natural convection[END_REF][START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF]. Five cases were computed (the exact values of the defining parameters are summarized in Table 3):

• PCM-Case #1 simulates the experimental study of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF]. It consists of a differentially heated square cavity, filled with octadecane paraffin.

• PCM-Case #2 is extracted from the collective publication by [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF], compiling the results of different numerical approaches used for the simulation of the melting of a PCM at high Rayleigh numbers.

• PCM-Case #3 simulates the melting of a cylindrical PCM with heated inner tubes, as in [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF].

• PCM-Case #4 simulates the melting of Gallium in a rectangular cavity heated by the side-wall, as in [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF].

• PCM-Case #5 simulates the solid crust formation in a highly distorted PCM domain, as in [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF].

We summarize in Table 3 the values of physical parameters and the scales used to simulate these cases. The values of the CPU time necessary to run each case (with default parameters) are given in Table 1

δt 1 10 -2 10 -2 10 -5 1 V ref ν l H ν l H ν l H ν l H ν l H Ra P r
Table 3: Parameters for the cases simulating the melting of a phase-change material.

PCM-Case #1:

Melting of an octadecane PCM in a square cavity [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] studied experimentally the melting of an octadecane PCM in a square cavity of height H = 1.5 cm. His results were often used to validate numerical methods [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF][START_REF] Wang | A comprehensive numerical model for melting with natural convection[END_REF][START_REF] Ma | Solid velocity correction schemes for a temperature transforming model for convection phase change[END_REF][START_REF] Danaila | A Newton method with adaptive finite elements for solving phase-change problems with natural convection[END_REF]. The material is initially solid (θ 0 = -0.01) and melts progressively starting from the left boundary, maintained at a hot temperature θ h = 1. The right boundary is also isothermal, with cold temperature θ c = -0.01. Horizontal boundaries are adiabatic. The other parameters of this case are reported in Table 3. The computation starts from a refined mesh near the hot boundary. Mesh adaptivity is applied every time step. Figure 5 shows in the left panel the adapted mesh at t = 80.1 and the 6 subdomains of the domain decomposition used for the parallel computation. The velocity vectors in the liquid phase, represented in the middle panel, show the strong convection cell affecting the shape of the liquid-solid interface. To offer a quantitative validation of this computation, we compare in Fig. 5 (right panel) the position of the solid-liquid interface at two time instants. Since δt = 1 in the simulation, we consider t = 40.1 and t = 79.1 and compare with experimental data of [START_REF] Okada | Analysis of heat transfer during melting from a vertical wall[END_REF] taken at slightly different time instants, t = 39.9 and t = 78.7. The obtained shape and position of the liquid-solid interface are very close to experimental results. We show in the same figure that the parallel computation (6 MPI processes) gives the same results as the sequential simulation (1 MPI process).

PCM-Case #2: Melting of an octadecane PCM with high Rayleigh number

This case considers the same problem of the melting of a PCM, but with a very high value of the Rayleigh number, Ra = 10 8 (see Table 3). This case is very challenging since the natural convection becomes important in the fluid flow, and enhances considerably the heat transfer. [START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF] compiled results provided by five different authors (Lacroix, Le Quéré, Gobin-Vieira, Delannoy and Binet-Lacroix) who used different numerical methods to compute the basic configuration presented in this section. This benchmark is a very demanding numerical test. The high velocity, inducing a very narrow thermal boundary layer, can lead to unrealistic results if under-resolved. This explains why two investigators among the five failed to correctly predict the process.

We show in Fig. 6 an illustration of the flow configuration at dimensionless time t = 5. The adapted mesh and the domain decomposition for the parallel computation with 6 MPI processes is presented in the left panel. The convection flow is depicted in the middle panel. Note the very narrow boundary layers, which impose a very refined mesh in these regions. This flow particularity explains why this case is rather costly in terms of computational time (see Table 1), when compared to the previous one. The right panel of Fig. 6 offers a quantitative validation of this computation. The position of melting front at t = 5 is in good agreement with that reported by Gobin and Le Quéré. Details of their numerical method are presented in [START_REF] Gobin | Melting from an isothermal vertical wall. synthesis of numerical comparison exercise[END_REF]. Gobin used a front-tracking method based on a coordinate transformation with a finite volume method and a 62 × 42 grid, while Le Quéré solved a single domain model using a second-order finite volume method with a 192 × 192 grid. The interest of the mesh adaptation is clearly demonstrated for this case, since we initially used a coarse 40 × 40 grid.

PCM-Case #3: Melting of cylindrical PCM with inner heated tubes

A more complex geometry, suggested by [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF], is simulated in this section. It consists of a cylindrical PCM of radius R = 1 with heated tube inclusions of different arrangements. This configuration is also interesting from a practical point of view. [START_REF] Agyenim | A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (lhtess)[END_REF] pointed out that more than 70% of the PCM containers used for heat storage are using shell-tube systems. We simulate three configurations, with one, four and nine heated tubes (see Fig. 7). The size of the tubes is adjusted to have the same total tube area for all configurations. The radius R i of the inner tube is R i = R/4 for the case with one tube, R i = R/8 for the four heated tubes case and R i = R/12 for the case with nine tubes. The inner tubes are heated at constant θ = θ h (Dirichlet boundary condition). A Neumann boundary condition (∂θ/∂n = 0) is used for the outer boundary. For the velocity, all boundaries are considered as non-slip walls (u = 0).

For the three configurations, Fig. 7 shows the adapted mesh and the domain decomposition with 6 subdomains (left column), the velocity vectors (central column) and the temperature field (right column). The time instants are chosen to have the same liquid fraction in the system, L f = 80%. Only half of the domain is simulated, since all configurations are symmetric with respect to the vertical axis. The mesh is refined initially around the inner tubes, and is dynamically adapted at each time step around the melting front and in the thermal boundary layer area.

To estimate the heat-storage efficiency of each configuration, we plot in Figure 8 the time evolution of the liquid fraction L f . Including more heated tubes results in an enhanced heat transfer and a faster melting process. The nine-tube configuration melts 5 times faster than the reference configuration with one tube. Note also from Figure 8 the good agreement between our results and those reported by [START_REF] Luo | Lattice Boltzmann simulation of convection melting in complex heat storage systems filled with phase change materials[END_REF], obtained using a completely different model based on the Lattice Boltzmann Method. We checked again that the results obtained with 6 MPI processes are identical to those computed with the sequential code.

PCM-Case #4: Melting of Gallium in a rectangular cavity

The melting of the Gallium in a rectangular cavity is a challenging case that generated an animated debate in the literature, since [START_REF] Dantzig | Modelling liquid-solid phase changes with melt convection[END_REF]. Le [START_REF] Quéré | A note on possible flow instabilities in melting from the side[END_REF] showed that the flow has to display a multi-cellular structure, resulting from the hydrodynamic instabilities during the conduction regime before the onset of convection. This multi-cellular flow configuration was also found numerically by [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF]. These authors pointed out that coarse grids or inconsistencies in the mathematical model could generate unphysical flows, with a mono-cellular convection cell. Therefore, this case is a relevant exercise to test the accuracy of our method. The parameters of this case are reported in Table 3. To capture the very small convection cells during the first step of the melting, [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF] used a 800 × 1, 120 fixed grid. With our adaptive method, a maximum of 4820 triangles are necessary to reproduce the numerical result of [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF]. The grid size is thus reduced with our method by a factor of 100. We illustrate in Fig. 9 the configuration of the phase-change system at dimensionless time t = 0.02. Note in the left panel the refined mesh at the liquid-solid interface and also at the borders of each convection cell. The multi-cellular structure of the flow is illustrated by the velocity field presented in the middle panel. The structure of the temperature field in the right panel shows that the rolls of the convection cells are well resolved. The number of cells decreases later through a process of roll merging, as it was also reported by [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF]. Our numerical results are in good agreement with the observations of [START_REF] Hannoun | Resolving the controversy over tin and gallium melting in a rectangular cavity heated from the side[END_REF], [START_REF] Cerimele | Numerical modelling of liquid/solid phase transitions: Analysis of a gallium melting test[END_REF] and Giangi and Stella (2000).

PCM-Case #5: Solid crust formation in a highly distorted domain

This final PCM case considers the solidification of a PCM in a domain with a more complicated shape (Fig. 10). This case was simulated by [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF], using a discontinuous Galerkin finite element method. The fluid is initially motionless, with an initial dimensionless temperature θ 0 = 2. The temperature of fusion is set to θ f = 1.4, according to [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF] (see Table 3 for the values of all parameters). The left side boundary is maintained at a cold temperature θ c = 1.39 in the initial stage. The right wall is isothermal, with hot temperature θ h = 2. A nearly steady-state natural circulation is induced in the early time evolution of the flow. Then, the cold temperature at the left wall is decreased to θ c = 1, below the temperature of solidification. At this point, the formation of a solid crust layer starts at the left boundary. Figure 10 depicts the configuration of the solid-liquid system at dimensionless time t = 30. Mesh adaptivity used metrics computed from the two components of the velocity and from the phase-change variables (enthalpy and temperature). Consequently, it is shown in the left panel that the mesh is well adapted near the solid-liquid front and following the distorted shape of the convection cell in the liquid phase. In the middle panel, we notice the high velocity values defining the convection cell, waving in the liquid phase. The final panel of Fig. 10 shows the temperature field, with the solid phase well identified in blue. This structure of the flow is in a very good agreement with that reported by [START_REF] Nourgaliev | Fully-implicit orthogonal reconstructed discontinuous galerkin method for fluid dynamics with phase change[END_REF].

Natural convection of water and water freezing

In this section we consider phase-change systems using water. Since pure water exhibits non-linear density variation for T < 10.2 • C, with a maximum at T m = 4.0293 • C, the Boussinesq force becomes non-linear. Convection and freezing of water are therefore interesting cases to test the flexibility of our Toolbox to deal with additional non-linear terms. We used the following density-temperature formula suggested by [START_REF] Gebhart | A new density relation for pure and saline water[END_REF]:

ρ(T ) = ρ m (1 -w |T -T m | q ) , (31)
with ρ m = 999.972 [kg/m 3 ], w = 9.2793 • 10 -6 [( • C) -q ], and q = 1.894816. Hence, the buoyancy term 1) is not any more linear and becomes after scaling:

f B = g(ρ ref -ρ)/ρ ref in (
f B (θ) = Ra P r Re 2 1 βδT ρ(θ f ) -ρ(θ) ρ(θ f ) , (32)
where β = (1/ρ m ) (dρ/dT ) is the thermal expansion coefficient taking the value β = 6.91 [START_REF] Scanlon | A numerical analysis of buoyancy-driven melting and freezing[END_REF]. We simulate a differentially heated square cavity filled with liquid pure distilled water. This problem was investigated experimentally and numerically by Giangi et al. (2000); [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; [START_REF] Michalek | Simulations of the water freezing process -numerical benchmarks[END_REF]. The non-dimensional parameters describing the problem are (see [START_REF] Michalek | Simulations of the water freezing process -numerical benchmarks[END_REF] for physical details): Ra = 2.518084 • 10 6 , Pr = 6.99 and Ste = 0.125.

• 10 -5 [(K) -1 ] (Scanlon

Natural convection of water

The initial temperature is linearly distributed in the square cavity, with a hot temperature T h = 10 • C at the left wall and a cold temperature T c = T f = 0 • C at the right wall. The temperature field and the streamlines of the steady state are presented in Figure 11a. The isoline θ = θ m , corresponding to the line of maximum density is represented by a dashed line. Due to the anomalous thermal variation of water density, two recirculating zones are formed in the flow: a lower (abnormal) recirculation in the vicinity of the cold wall where θ < θ m and an upper (normal) one where the density decreases with temperature (θ > θ m ). A more precise comparison with previously published results is shown in Figure 11b. The obtained temperature profile θ(x) along the horizontal central line of the cavity (y = 0.5) is in good agreement with the numerical results of [START_REF] Michalek | Simulations of the water freezing process -numerical benchmarks[END_REF]. Their results were obtained with finite-volume and finite-difference codes (FLUENT and FRECONV3V). 29

Water freezing

We finally consider the difficult case of water freezing in a square cavity. The initial state for this computation is the convection steady pattern in the cavity presented in Figure 12 shows the flow configuration at dimensionless time t = 1.61. The left panel depicts the domain composition and the adapted mesh following two lines of interest: the solid-liquid interface T = 0 and the line T = T m separating the two recirculating zones (see the middle panel of the same figure). The metrics used for adaptivity were computed from the two components of the velocity, the temperature and a P 1 tanh-function φ(T ) "tracking" the value T m . To reduce the impact of the interpolation on the global accuracy we used both φ(T n ) and φ(T n+1 ) in the adaptivity procedure (see also Belhamadia et al. (2004a)). This adaptivity strategy allowed us to accurately capture the structure and the extent of the two recirculating zones, features that are difficult to obtain with fixed meshes (discrepancies in numerical results are described in Giangi et al. (2000); [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]; [START_REF] Michalek | Simulations of the water freezing process -numerical benchmarks[END_REF]). In the right panel of Fig. 12 we superimpose the experimental image from [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF] with our numerical results for the same physical time t ϕ = 2340 [s]. The position of the solid-liquid interface and the flow pattern in the liquid phase correspond very well qualitatively to the experimental image. The simulation was performed with a small time step (δt = 10 -2 ≈ 15[s]). Note from Table 1 that the computational time for this case is more than reasonable (23 minutes using 1 thread and only 13 minutes with 6 MPI processes).

3D parallel numerical results

All the 2D benchmarks presented in the previous section could be simulated in 3D using the second provided Toolbox PCM Toolbox DDM 3D. To avoid an excessively long paper, we shortly present in this section only the following selection of 3D cases: (i) the steady natural convection of air in a cubic cavity with Ra = 10 6 ( §6.1), (ii) the melting of a phase-change material ( §6.2), within a cubic ( §6.2.1) or cylindrical domain ( §6.2.2), and (iv) the convection and the freezing of pure water in a cube ( §6.3). Running parameters for these cases are summarized in Table 4 andTable 5, for computations using P 1 and P 2 finite elements, respectively. Again, to make these simulations accessible with fast-running jobs in a batch-queuing system, we gave priority to jobs with a reasonable number of MPI processes (less than 400). We used 4GB memory for each CPU-Core node for all simulations, excepting for the case PCM-Case #3 with 9 heated tubes, when 9GB of memory were necessary. In the subsequent sections, we present the results obtained with P 2 finite elements for the temperature, since these simulations are more challenging. The results using P 1 elements for the temperature are provided in the Supplemental Material http://lmrs-num.math.cnrs.fr/2020CPCP2.html containing images and animations of the cases presented in this section. 4, but with all computations performed using P 2 finite elements for the temperature.

Case

Natural convection of air in a cubic cavity

We go back to the benchmark of natural convection of air described in §5.1 and add the third dimension in space. We thus simulate the thermally driven cubic cavity [0, 1] 3 , filled with air. The temperature is now fixed on the left (cold) wall surface and the right (hot) wall surface. All the other lateral surfaces are adiabatic (∂T /∂n = 0). No-slip walls are applied for the velocity on all boundary surfaces. We show in Fig. 13 the final (steady) flow configuration for Ra = 10 6 obtained using the time-dependent solver with P 2 finite-elements for the temperature. The mesh is mainly adapted following the boundary layers of the ascending and descending flows (Fig. 13a), which are clearly visible when representing the 3D velocity vectors (Fig. 13b). A quantitative assessment of the quality of the simulation is offered in Fig. 13c and Fig. 13d, comparing our results with those reported by [START_REF] Wakashima | Benchmark solutions for natural convection in a cubic cavity using the high-order time-space method[END_REF]. These authors solved the vorticity-stream function formulation of the Navier-Stokes equations using a fourth-order finite difference scheme and a uniform mesh of 120 3 grid nodes.

Melting of 3D phase-change materials

We simulate in this section 3D configurations of PCMs, in cubic or cylindrical geometries.

Melting of a cubic octadecane PCM

We extend the 2D simulations presented in §5.2.1 to a 3D cubic domain Ω = [0, 1] 3 . We impose cold dimensionless temperature θ c = -0.01 at x = 1 (right wall), hot temperature θ h = 1 at x = 0 (left wall), and a homogeneous Neumann boundary condition at the remaining walls. A Dirichlet boundary condition u = 0 is prescribed for all surfaces defining ∂Ω. The temperature distribution and the corresponding adapted mesh at t = 80.1 is shown in Fig. 14a. The blue region denotes the solid phase. A zoom of the mesh in the mid-plane, at t = 80.1, refined along the iso-surface θ = 0, is shown in Fig. 14b. As the heating progresses, the natural convection intensifies enough to have a pronounced influence on the shape of the interface. The temperature difference between the hot wall and the solid PCM induces a clockwise recirculation of the flow in the melted PCM (Fig. 14c). The shape of the liquid-solid interface at t = 80.1 is also displayed in Fig. 14c, showing a non-uniform melting front receding from the top to the bottom of the domain. The position of the solid-liquid interface at the mid-plane for 2D and 3D configurations is plotted in Fig. 14d. Differences between 2D and 3D results are not visible in this mid-plane. Three-dimensional effects are more important near the lateral walls, as explained in [START_REF] Nikrityuk | Three-dimensionality of fluid flow in the benchmark experiment for a pure metal melting on a vertical wall[END_REF]. comparison between 2D and 3D simulations (both using P 2 finite elements for the temperature).

Melting of a cylindrical PCM with inner heated tubes

We simulate the melting of a cylindrical PCM with inner heated tubes, as in §5.2.3. The 3D domain is now a horizontal cylinder of radius R = 1 and length L = 0.25. In addition to the 2D configuration, we impose the following boundary conditions. For the temperature, a Dirichlet boundary condition (θ = θ h = 1) is applied at the inner tubes and a Neumann boundary condition ( ∂θ ∂n = 0) at the outer sections of the cylinder. For the velocity, a homogeneous Dirichlet boundary condition (u = 0) is prescribed on the lateral surfaces of external and inner cylinders and a symmetry boundary condition is imposed on the ending sections of the cylinder. The PCM is initially solid and when the melting process starts, a liquid layer grows around the heated tubes and expands toward the lateral boundary. A slice in the adapted mesh, the velocity vector field, and temperature iso-surfaces are illustrated in Fig. 15 for configurations with one, four, and nine heated tubes. The mesh is adapted along the melting front and in the boundary layers around the heated tubes. Adaptivity is efficient even when the solid-liquid interface touches the outer wall (see left panels). Counter rotating flows are observed at the left and right parts of the cylindrical domain (middle panels). Iso-surfaces of temperature shown in the right panels also confirm that the evolution of the liquid flows keeps the symmetry of the configuration.

Natural convection of water and water freezing

We simulate in this final section the challenging case of the natural convection of water and the ice formation inside a differentially heated cubic cavity.The dimensionless parameters are the same that those used in §5.3.2. For the natural convection of water, the vertical walls at x = 0 and x = 1 are isothermal and have different temperatures θ h = 1 and θ c = 0, respectively. The remaining walls are adiabatic. A homogeneous Dirichlet boundary condition u = 0 is applied on all walls. The fluid is initially at rest and the temperature is linearly distributed from the cold to the hot wall (see Fig. 16). Given the anomalous thermal variation of water density, we also adapted the 3D mesh along θ = 0.4 to capture correctly the flow structure. The temperature distribution and the corresponding adapted mesh for the steady (time-independent) state computation are shown in Fig. 16a. The blue region identifies the cold water trapped by the abnormal fluid recirculation and the red region the hot fluid driven by the upper clockwise circulation. A cut in the mesh following the mid-plane (Fig. 16b) shows the internal structure of the mesh. Smaller tetrahedra are clearly observed in the vicinity of the walls and between the two counter-rotating circulation patterns. Figure 16c shows that the three-dimensional evolution of the flow (spiral movement of the fluid along the walls) affects the topology of the iso-surface θ = 0.4. Giangi et al. (2000) analyzed three-dimensional effects in both convection and freezing of water and they noted that only the flow in the symmetry plane is similar to the two-dimensional flow. They pointed out that no-slip velocity and adiabatic thermal boundary conditions at the side walls enhance three-dimensional effects near the walls. To validate this observation, we plot in Fig. 16d the profile of the vertical velocity along the z direction, at x = 0.93 and passing through the velocity saddle point, where normal and abnormal convection streams collide in the vicinity of the cold wall. We compare both 2D and 3D simulations with 3D results of [START_REF] Michalek | Simulations of the water freezing process -numerical benchmarks[END_REF] (symbols), obtained with a finite-difference code. 3D (red solid line) profile agrees well with the benchmark solution with a maximum difference of 3%.

For the water freezing case, the simulation starts from the steady solution shown in Fig. 16. The cold temperature at the right wall is suddenly dropped from θ c = 0 to θ c = -1 (corresponding to physical temperature T c = -10 o C). Ice starts to form at the cold wall and expands toward the left wall. The temperature distribution and the adapted mesh at the final time t ϕ = 2340[s] (t = 1.61) are displayed in Figs. 17a and17b.

When compared to the natural convection case (Fig. 16), the mesh is also adapted along a second moving front represented by the solid-liquid interface (iso-surface θ = 0 shown in Fig. 17c). For the identification of flow structures during the ice formation, Fig. 17c shows the 3D velocity vectors, the iso-surface θ = 0.4 along the anomalous density variation, and the phase-change front at time instant t = 1.61. The effect of the secondary flow (see [START_REF] Nikrityuk | Three-dimensionality of fluid flow in the benchmark experiment for a pure metal melting on a vertical wall[END_REF]) is visible from the curved shape of the iso-surface θ = θ m in the transverse y direction. We note, however, that the shape of the solidification front is almost 2D. The buoyancy-induced fluid motion in the abnormal recirculation region is too weak to influence the solid front. Finally, the superimposition of the experimental image of [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF] and the current simulation give good agreement for the location of the solid-liquid interface. Differences come mainly from the fact that the undercooling phenomenon during the solidification stage is not taken into account in our physical model. [START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]. Simulations using 224 MPI processes with P2 for the temperature and a max of degrees of freedom 5, 107 229.

Scalability test for 3D parallel computations

We present the results for both strong and weak scalability tests. We simulate the 3D unsteady natural convection of air with Ra = 10 6 presented in §6.1. Mesh adaptation is used only for the strong scalability test.

The CPU is Intel Broadwell E5-2680 v4 @ 2.40GHz (14 cores per socket). Each node has two sockets and 128 GB of DDR4 2400 MHz RAM. They communicate through Intel Omnipath 100Gb/s low latency network. First, we report the simulation time, the efficiency and speed-up for the weak (Fig. 18) and strong scalability tests (Fig. 19). Assuming a perfect speed-up for 8 cores, we obtained for 256 cores a speed-up of approximately 122.56 (weak) and 70.6 (strong), resulting in an efficiency of 47.9%(weak) and 27.6%(strong). In Figs. 20 (weak) and 21 (strong), we report the computational time for the main steps of the computational algorithm, namely the mesh adaptation, the matrix construction and the resolution of the linear system. For 8 cores, the mesh adaptation takes 221.9s (strong), the construction of the matrix 83.5s (strong) and 0.96s (weak) and the resolution of the linear system 4866.9s (strong) and 552.1s (weak). For 256 cores, the mesh adaptation takes 223.0s (strong), the construction of the matrix 76.4s (strong) and 6.91s (weak) and the resolution of the linear system 422.2s (strong) and 893.5s (weak). It results from this analysis that the only step that benefits from the parallelization is the resolution of the use of an adaptive finite element method with regularized functions representing the variation of thermodynamic properties at the solid-liquid interface, (iii) a fully implicit discretization with a Newton algorithm for solving the non-linear system of equations, and (iv) a parallel final algorithm using automatic partitioners (Scotch or Metis) and the recent library ffddm that makes available in FreeFem++state-of-the-art scalable Schwarz domain decomposition methods. Mesh adaptivity is an essential ingredient of the toolbox and a new algorithm was implemented using external remeshing tools (mmg) for 3D configurations.

We presented careful validations of the toolbox by considering parallel 2D and 3D computations of well-known benchmark cases of increasing difficulty: natural convection of air, natural convection of water, melting or solidification of a phase-change material, and, finally, a water freezing case. For 2D configurations, parallel computations with the present toolbox using modest resources (6 MPI processes) can bring a CPU gain factor between 2 and 3, when compared to previous sequential algorithms [START_REF] Rakotondrandisa | A finite-element toolbox for the simulation of solid-liquid phase-change systems with natural convection[END_REF]. But the most significant progress brought by the new toolbox is to make affordable well-resolved 3D computations with mesh adaptivity, which are hardly accessible with sequential algorithms. The efficiency of the toolbox in 3D was demonstrated by simulating difficult cases (melting of a PCM, water freezing) in simple (cubic) or complex geometries (cylindrical with inner tubes). To facilitate its use, the toolbox is provided with separate folders containing all the necessary files (parameters, restart files) necessary to run all the cases described in the paper. Ready-made scripts and layouts allow the user to generate the figures presented in this paper with newly generated data after running the programs. Validation data sets from experiments or previous publications are included in these layouts.

The present parallel approach can be easily tested and adapted to address other computational challenges related to different physical or mathematical models in this field. For instance, other possible choices for the finite-element discretization could be tested for the implementation of the algorithm presented in this paper. To prove the versatility of the toolbox, we illustrate in Fig. 22 the results obtained using the mini-element introduced by Arnold, Brezzi and Fortin (see [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF][START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. This is the simplest element for Stokes-type problems, offering inf-sup stability and global linear convergence. Compared to the Taylor-Hood element, in the mini-element discretization velocities are piecewise P b 1 (P 1 -bubble), while pressure is still piecewise P 1 . Since the mini-element already exists in FreeFem++ the only necessary change in the distributed scripts is to declare the unknowns of the problem accordingly to the syntax discussed in Sect. 3.1: fespace Wh(Th,[P1b,P1b,P1b,P1,P1]).

Figure 22 shows that, for the test case of the melting of a PCM presented in Sects. 5.2.1 (2D case) and 6.2.1 (3D case), the results are remarkably accurate when compared to those obtained using the Taylor-Hood discretization. The CPU time is reduced by 20% for the 2D computation and 42% for the 3D simulation with the same number of MPI processes (6 for 2D and 224 for 3D). However, knowing that the mini-element offers a minimal accuracy for the velocity field, a careful assessment of this type of discretization in computing more complicated cases (with several convection rolls and complex dynamics in the fluid part) is necessary. Also, the efficiency of the mesh adaptivity in this case has to be carefully investigated to obtain a robust algorithm (as in the present toolbox using the Taylor-Hood finite element). This topic will be addressed in a future contribution, together with the possibility to enrich the physical model for the solidification process.

Figure 1 :

 1 Figure 1: Sketch of the liquid-solid phase change problem considered in this paper. (a) Illustration of the melting, with convection in the liquid phase (streamlines of the velocity field) and a bent interface separating the two phases. (b) Illustration of the domain-decomposition used in present numerical simulations. Coloured patches represent subdomains created by automatic graph partitioning libraries. Note the adapted mesh, especially around the liquid-solid interface.

Figure 2 :

 2 Figure 2: Folder tree structure of the PCM Toolbox DDM 2D and PCM Toolbox DDM 3D to solve solid-liquid phase-change problems.

Figure 3 :

 3 Figure 3: Structure of each Test-case folder. Solid lines in the diagram correspond to folders which come with the Toolbox, the dashed lines to folders generated after running the script NEWTON case.edp, the hashed lines to folders generated after running the script NEWTON case.edp in the 2D case. Finally, the dotted lines indicate folders generated after running the script Read 3D data.edp, which reads data from the RST folder and generates the Medit, Paraview and Tecplot files for visualization.

Figure 4 :

 4 Figure 4: Natural convection of air in a differentially heated cavity for values of Ra ranging from 10 4 to 10 6 . (a) Vertical velocity profile v(x) along the horizontal symmetry line (y = 0.5). (b) Longitudinal velocity profile u(y) along the vertical symmetry line (x = 0.5). Results obtained using the present Toolbox, with a mesh resolution nbseg = 80. Comparison with the spectral simulations by Le Quéré (1991) (solid lines).

Figure 5 :

 5 Figure 5: PCM-Case #1. Melted PCM at time instant t = 80.1. Adapted finite-element mesh and subdomains used for the domain decomposition and parallel computing with 6 MPI processes (left panel). Velocity vectors in the liquid phase (central panel). Position of the solid-liquid interface and comparison with experimental data of Okada (1984) for two time instants (right panel).

Figure 6 :

 6 Figure6: PCM-Case #2. Melting of the PCM with high Rayleigh number (Ra = 10 8 ). Simulated configuration at time instant t = 5. Adapted finite-element mesh and subdomains used for the domain decomposition and parallel computing with 6 MPI processes (left panel). Velocity vectors in the liquid phase (central panel). Position of the solid-liquid interface and comparison with five sets of results presented in[START_REF] Bertrand | Melting driven by natural convection a comparison exercise: first results[END_REF].

Figure 7 :

 7 Figure 7: PCM-Case #3. Melting of a cylindrical PCM with heated inner tubes. Flow configuration for one tube at t = 2.491 (first raw), four tubes at t = 0.981 (middle row) and nine tubes at t = 0.411 (last row). Time instants correspond to the same liquid fraction in the system, L f = 80%. Adapted mesh and domain decomposition with 6 subdomains (left column), velocity vectors with background color representing the velocity intensity (central column) and temperature field (right column).

Figure 8 :

 8 Figure 8: PCM-Case #3. Melting of a cylindrical PCM with heated inner tubes. Time evolution of the liquid fraction for configurations with one, four, and nine heated tubes. Comparison with numerical results of Luo et al. (2015).

Figure 9 :

 9 Figure 9: PCM-Case #4. Melting of Gallium. Simulated configuration at t = 0.02. Adapted finite-element mesh and subdomains used for the domain decomposition and parallel computing with 6 MPI processes (left panel). Velocity vectors in the liquid phase (central panel). Temperature field (right panel).

Figure 10 :

 10 Figure 10: PCM-Case #5. Solid crust formation in a distorted domain. Simulated configuration at t = 29.33. Adapted finite-element mesh and subdomains used for the domain decomposition and parallel computing with 6 MPI processes (left panel). Velocity vectors in the liquid phase (central panel). Temperature field (right panel).

Figure 11 :

 11 Figure 11: Natural convection of water in a differentially heated cavity. Non-dimensional temperature θ at steady state. (a) Two-dimensional temperature field and streamlines showing the two recirculating zones. (b) Temperature profile along the horizontal central line. Comparison with the numerical results of Michalek and Kowalewski (2003).

Figure 11 .

 11 The freezing starts by dropping smoothly the temperature of the cold right wall from T c = 0 • C to T c = -10 • C. The new boundary condition on the right cold wall is imposed by setting a very thin layer of δx = 0.01, with constant temperature T = T c and zero velocity.

Figure 12 :

 12 Figure 12: Freezing of pure water. Simulated configuration at t = 1.61. Adapted finite-element mesh and subdomains used for the domain decomposition and parallel computing with 6 MPI processes (left panel). Velocity vectors in the liquid phase (central panel). Comparison with the experimental image from Kowalewski and Rebow (1999); the thick red line represents the solid-liquid interface computed with the present method (right panel).

Figure 13 :

 13 Figure 13: Natural convection of air in a cubic differentially heated cavity (Ra = 10 6 ) computed with P 2 finite elements for the temperature. (a) Adapted mesh with a maximum number of degrees of freedom of 913 057. (b) 3D velocity vector field in the final (steady) state. Temperature contour-lines in the central section of the cube (y = 0.5): results of Wakashima and Saitoh (2004) (c) and present results (d). Parallel computation with 56 MPI processes.

Figure 14 :

 14 Figure 14: Melting of a cubic PCM (Okada, 1984) computed with P 2 finite elements for the temperature. Configuration at t = 80.1, computed using 224 MPI processes with P 2 finite elements for the temperature. (a) Adapted finite-element mesh with a maximum number of degrees of freedom ndof = 4, 871 640. (b) Illustration of the adapted mesh in the centre of the domain. (c) 3D velocity vectors and isosurface θ = 0 representing the liquid-solid interface. (d) Position of the interface in the central plane (y = 0.5): comparison between 2D and 3D simulations (both using P 2 finite elements for the temperature).

Figure 15 :

 15 Figure 15: Melting of cylindrical PCM computed with P 2 finite elements for the temperature. Adapted finite-element mesh (left panels), 3D velocity vector field (middle) and temperature iso-surfaces (right). Configurations with one inner heated tube (time instant t = 1.591, simulation with 224 MPI processes and a maximum number of degrees of freedom ndof = 3, 919 350), four tubes (t = 0.621, 280 processes and max ndof = 9, 224 708) and nine tubes (t = 0.331, 364 processes and max ndof = 14, 981 519).

Figure 16 :

 16 Figure 16: Natural convection of water in a differentially heated cubic cavity computed with P 2 finite elements for the temperature. (a) Temperature distribution and 3D adapted finite element mesh at the steady state. (b) Illustration of the adapted mesh in the center of the domain. (c) Velocity vector field and temperature iso-surface θ = 0.4. (d) Profile of the vertical velocity along the z-direction at the mid-plane (y = 0.5) and x = 0.93. Comparison with the numerical benchmark suggested by Kowalewski and Rebow (1999) (symbols). Simulations using 112 MPI processes with P2 for the temperature and a maximum number of degrees of freedom ndof=4, 700 233

Figure 17 :

 17 Figure 17: Freezing of pure water in a 3D cubical cavity computed with P 2 finite elements for the temperature. (a) Temperature distribution and adapted mesh at time instant t = 1.61. (b) Cut through the adapted mesh at mid-plane. (c) 3D velocity vectors in the liquid phase and temperature iso-surfaces θ = 0.4 and θ = 0. (d) Superimposition of the interface obtained in the present simulation (red line) on the experimental image of[START_REF] Kowalewski | Freezing of water in differentially heated cubic cavity[END_REF]. Simulations using 224 MPI processes with P2 for the temperature and a max of degrees of freedom 5, 107 229.

Figure 18 :

 18 Figure 18: Weak scalability test for 3D simulations: simulation time (left), efficiency (center) and speed-up (right) for a number of processes ranging from 8 to 256.

Figure 19 :

 19 Figure 19: Strong scalability test for 3D simulations: simulation time (left), efficiency (center) and speed-up (right) for a number of processes ranging from 8 to 256.

Figure 22 :

 22 Figure 22: Melting of a PCM. Results obtained using the mini-element (P b 1 for the velocity and P 1 for the pressure). 2D simulation: (a) Temperature field and position of the interface at t = 80.1; the red dashed line is the experimental result of Okada (1984). (b) Comparison with the Taylor-Hood discretization (red dashed lines) for the position of the interface at t = 40.1 and t = 79.1. 3D simulation: (c) Temperature field and adapted mesh at t = 80.1. (d) Comparison with the Taylor-Hood discretization (red dashed line) for the position of the interface at t = 79.1 in the central plane (y = 0.5). To be compared with Figs. 5 and 14.

Table 1 :

 1 Comparison of CPU times for 2D simulations performed with one or 6 MPI processes. Speed-up of the parallel computation, maximum number of triangles of the adapted mesh, total number of time steps and ratio between the CPU time for adapting the mesh and the CPU time for a complete time step.

		CPU time 1 MPI proc.	CPU time 6 MPI procs.	speed-up	max number of trian-gles	number of time steps	ratio CPU adapt / CPU time step
	Unsteady NC of air,Ra = 10 6	00:01:39	00:00:50	1.98	2,902	53	18%
	PCM-Case #1	00:06:39	00:04:35	1.45	2,590	81	5%
	PCM-Case #2	01:31:11	00:40:44	2.23	5,928	501	8%
	PCM-Case #3: 1 tube	00:25:09	00:19:18	1.3	5,052	401	6%
	PCM-Case #3: 4 tubes	00:28:57	00:16:13	1.78	13,507	151	6%
	PCM-Case #3: 9 tubes	00:33:28	00:17:44	1.88	24,866	100	5%
	PCM-Case #4	06:31:31	01:38:35	3.97	5,371	2,001	15%
	PCM-Case #5	00:07:52	00:03:17	2.39	2,898	117	12%
	Cycle of a PCM	00:13:58	00:08:49	1.58	3,160	383	17%
	NC of water	00:00:36	00:00:21	1.71	2,686	20	14%
	Unsteady NC of water	00:01:49	00:00:54	2.02	3,368	33	10%
	Water freezing	00:22:49	00:13:38	1.67	3,743	296	10%
	All computations were performed on a Macbook pro 2.2GHz Intel Core i7, 16GB of DDR4 2400 MHz
	RAM.

Table 2 :

 2 Natural convection of air in a differentially heated cavity. Maximum value umax of the horizontal velocity profile at mid-domain (x = 0.5) and location Y of this maximum. Comparison to reference values reported by Le[START_REF] Quéré | Accurate solutions to the square thermally driven cavity at high Rayleigh number[END_REF].

  .

		case #1 case #2 case #3 case #4 case #5
	Ra 3.27 • 10 5	10 8	5 • 10 4	7 • 10 5	10 6
	P r	56.2	50	0.2	0.0216	0.1
	Ste	0.045	0.1	0.02	0.046	4.854

Table 4 :

 4 Running parameters for 3D simulations: total CPU time, maximum number of degrees of freedom (dof) after mesh adaptivity, number of MPI processes, total number of time steps, maximum CPU time necessary for mmg library to adapt the 3D mesh and ratio between the CPU time for adapting the mesh and the CPU time for a complete time step. All computations were performed using P 1 finite elements for the temperature. All computations were performed using a parallel computer (CRIANN Computing Center and MATRICS platform) based on Intel Broadwell E5-2680 v4 @ 2.40GHz (14 cores per socket) architecture with two sockets per node and 128 GB of DDR4 2400 MHz RAM. An Intel Omnipath 100Gb/s low latency network was used for communications.

		CPU time	max of dof	MPI procs.	number of time steps	max mmg CPU time	ratio CPU adapt / CPU time step
	Unsteady NC of air,Ra = 10 4	00:55:21	249 625	6	23	9 sec.	6%
	Unsteady NC of air,Ra = 10 6	01:14:49	725 097	56	50	20 sec.	19%
	PCM-Case #1	07:42:32	3,299 661	224	80	53 sec.	16%
	PCM-Case #3: 1 tube	16:02:31	3,809 760	224	401	65 sec.	17%
	PCM-Case #3: 4 tubes	13:29:23	8,813 510	280	151	150 sec.	20%
	PCM-Case #3: 9 tubes	10:37:44	14,860 117	364	71	262 sec.	24%
	Cycle of a PCM	13:06:26	2,150 052	224	382	36 sec.	32%
	NC of water	00:59:26	3,701 916	112	20	62 sec.	20%
	Unsteady NC of water	06:39:40	5,492 302	112	37	76 sec.	14%
	Water freezing	1-10:11:55	4,448 031	224	296	53 sec.	21%

Table 5 :

 5 Same caption as for Table

FreeFem++for different OS can be downloaded from http://www.freefem.org/.

http://www.mmgtools.org/.

Acknowledgements

This project was co-financed by the European Union with the European regional development funds and by the Normandy Regional Council via the M2NUM (ERDF, HN0002081) and M2SiNUM (ERDF, 18P03390/18E01750/18P02733) projects. Part of this work was performed using computing resources of CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, France) and MATRICS plateform (Amiens, France).

linear system. The speed-up of the resolution of the linear system is 78.9 (strong) and 158.2 (weak), resulting in an efficiency of 36.0% (strong) and 61,8% (weak).

Summary and conclusions

The main advantage of the new toolbox distributed with this paper is to enable parallel computing of liquid-solid phase-change systems in 2D or 3D complex geometries. The physical model implemented in the software is one of the most accurate existing for such systems: the enthalpy-porosity model, based on the incompressible Navier-Stokes equations with Boussinesq approximation for thermal effects. This model captures all the characteristics of the flow developing in the liquid phase, dominated by convection. Using a Carman-Kozeny-type penalty term, added to the momentum equations to bring progressively the velocity to zero into the solid, results in a single-domain formulation of equations for both liquid and solid phases. We proved that this single-domain enthalpyporosity model is well adapted to the use of Schwarz domain decomposition methods (DDM), and, consequently, to parallel computing.

The numerical method was implemented using FreeFem++, a free software offering multiple advantages to easily implement a finite-element algorithm. The key ingredients of the implemented method are: (i) a second order accuracy in space (P 2 finite elements for velocity and temperature, P 1 for the pressure) and time (Gear scheme); (ii) the Supplementary images and movies depicting the dynamics of some cases simulated in this paper are provided as Supplemental Material at http://lmrs-num.math.cnrs.fr/2020CPCP2.html.