ConvectiveFoam1.0: development and benchmarking of a infinite-Pr number solver

We developed a new fluid-dynamical numerical model, which we call convectiveFoam, designed to simulate fluids with very large Prandtl number. First we implemented the high-Pr case, in which advection still acts explicitly, and then the Pr→∞ version, where the momentum equation becomes diagnostic (that is, without time derivatives) and it is formalized as an elliptic problem. The new solver, based on a finite volume integration method, is developed on the OpenFOAM platform 5 and it exhibits a good performance in terms of computational costs and accuracy of the results. Scaling properties show a maximum performance around 16000 cells/core, in agreement with other works developed on the same platform. A systematic validation of the solver was performed for both 2D and 3D geometries, showing that convectiveFoam is able to reproduce the main results of several iso-viscous cases. This new solver can thus simulate idealized configurations of natural geophysical convection, such as in the Earth Mantle where Pr = 10. This solver represents a starting point for general exploration of the 10 behaviour and parameter dependence of several fluid systems of geological interest.


Introduction
Convection in the Earth mantle is driven by internal heating associated with radioactive decay, cooling at the surface, crustal subduction and interaction between rock and liquid water. Mantle convection in turn affects plate tectonics and the global carbon cycle. Since our planet is the only known one which shows an active tectonics, several researchers also suggested a link 15 between this phenomenon and Earth habitability. Despite such crucial importance in determining planetary equilibria, however, mantle convection is still only partially understood, owing to the extreme difficulties in examining planetary interior dynamics (Bercovici and Ricard, 2014;O'Neill et al., 2007). To cite a few, the Earth mantle exhibits a viscoelastic nature, behaving like a solid on short time scales but like a fluid on the long ones. This behaviour is synthesized by the Maxwell time: 20 defined as the ratio between dynamic viscosity ([η] = Pa · s) and elasticity (or shear modulus [µ R ] = Pa). A typical value for the Earth mantle is τ M ≈ 10 10 s which relates the evolution to the geological time scales. Also the characteristic spatial scales (about 3000 km) are a non trivial obstacles to direct geological measurements (Bercovici and Ricard, 2014). For these reasons, indirect methods have usually been employed, such as seismic data analysis or numerical simulations. In particular, for the long time-scale behaviour (t τ M ), numerical simulations produce diagnostic and prognostic descriptions of geological processes 25 (Khaleque et al., 2015). Furthermore, because of the high viscosity values (ν = 10 18 − 10 21 m 2 s −1 ), the infinite-Pr number assumption, where Pr = ν κ , ν is the kinematic viscosity and κ the thermal diffusivity, is usually adopted for the Earth mantle (Bercovici, 2007;Turcotte and Schubert, 2014). Simulating high-Pr regimes poses peculiar numerical difficulties due both to time steps constraints and changes of the mathematical structure of governing equations (Busse, 1979).
This work is a first step towards the exploration of Earth Mantle processes through a model suitable to simulate both high-30 and infinite-Pr fluid flows. In its first implementation, we will present the model for an isoviscous, single-phase fluid. However, the numerical model is conceived to include future applications with non-Newtonian viscosities (dependent on temperature and pressure), multi-phase and multi-component flows, and for non Cartesian problems, potentially suitable to simulate complex geometries (e.g., slabs or spherical shells) and flows with strong interfaces, needing locally refined numerical meshes.
To develop the new tool for the numerical simulations we adopt the OpenFOAM® (Open Field Operation and Manipulation) 35 toolkit, a free, open source, and widely diffused CFD software. OpenFOAM offers a variety of C++ libraries for Finite Volume (FV) discretisation of partial differential equation systems on a three-dimensional unstructured, colocated mesh. It includes libraries for data manipulation and linear algebra, and a number of numerical solvers for the solution of CFD problems. Thanks to the high level of abstraction, numerical solvers can be developed following pre-built templates. The segregated, iterative solution method at the base of most pre-built solvers allows the solution of CFD problems with complex, non-linear rheolo-40 gies, compressible/incompressible regimes, turbulent flows, and the solution of multi-phase, multi-component fluids, withouth increasing the size of the linear systems, thus keeping the computational complexity affordable. OpenFOAM is parallelized with a domain decomposition approach using MPI libraries. It displays satisfactory parallel efficiency up to thousands of cores and offers a variety of pre-and post-processing integrated tools. OpenFOAM is released open-source and is supported by a broad users and developers community worldwide, making it suitable for the development of community models.

45
All built-in OpenFOAM solvers are implemented for low-Pr viscous fluids and several rheological laws are already implemented for fluids like water, air, oil, but also for multiphase mixtures of gas, liquids and granular fluids (OpenCFD, 2007). Some extensions of this provided solvers have been recently developed and applied to Earth Science (Cerminara et al., 2016;Dietterich et al., 2017;Rosi et al., 2018). Since an infinite-Pr solver is not provided in the current version of OpenFOAM, we built a new solver to simulate idealized problems in geological convection. We started exploring high-Pr fluid behaviour, up to 50 that of infinite-Pr number fluids (under the assumption of isoviscosity). For infinite-Pr number fluids, we implemented a brand new solver, which is then validated for 2D and 3D configurations.
The structure of this paper is as follows: a theoretical description of the problem is presented in Sec.2; Sec.3 introduces the numerical setting with an overview on code development. The validation of the new solver is presented in Sec.4. Conclusions and perspectives are given in Sec.5.
As in the milestone work on convection done by Rayleigh (Rayleigh, 1916) and Bénard (Bénard and Avsec, 1938) at the beginning of last century, the convection model adopted here is based on classical (incompressible) fluid dynamics equations (Busse, 1979;Chandrasekhar, 1981). The equations for momentum, energy and mass conservation in a rotating frame read: where u = (u, v, w) is the three-dimensional velocity, Ω = (0, 0, Ω) is the angular velocity, T is the temperature field and other parameters are listed in 1 Pr where primed variables are non-dimensional. In Eqs.
(3), Pr = ν κ and Ra = gβ∆T D 3 νκ are the Prandtl and Rayleigh numbers,70 representing respectively the ratio between convective and conductive heat transfer and the ratio between kinematic and thermal heat transport.  (Ricard, 2007). parameters are strongly dependent on pressure and temperature, here we assume constant parameters, fixed as in Ricard (2007).
This leads to constant values for Pr and Ra. Moreover, since Ek ≈ 10 12 , the Pr → ∞ assumption holds. In this limit, Eqs. (3) become (removing primes for simplicity): where an elliptical PDE governs the dynamics of the velocity field. In this new problem, the velocity field is strictly dependent on the temperature and pressure fields, since acceleration has been neglected together with the nonlinear inertial term. Never-80 theless, the advective term in the temperature equation still acts as a source of nonlinearity. It follows that analytical solutions are not available and a numerical approach is warranted.

Numerical Methods
The problem described by Eqs.(4) is solved with OpenFOAM. OpenFOAM uses the FV method discretising the integral form of governing equations over each control volume (CV). Each CV is indentified by a polyhedral cell and all field values are 85 calculated in the centre of such cell (i.e. colocated) (Jasak, 1996). Once the mesh is defined, the surface fluxes are reconstructed using the Gauss theorem. This requires an interpolation step to infer the surface value from the volume one. To this aim OpenFOAM provides a set of pre-defined interpolation rules (Ferziger and Peric, 2012). In the following we briefly introduce two rules adopted in our work.  Jasak (1996). Each field is calculated in the centre of a polyhedral cell ( P). N is the centre of the nearby cell and f the point on the surface on which the field is reconstructed by interpolating between P and N.
The first is the Linear method, also known as Central Differencing (CD), which reconstructs the field on the surface S, 90 assuming its linear variation between the two points P and N, located in the centre of two nearby cells. The field value in point f is calculated as: where f x is defined as the ratio between f N and P N . CD is second order accurate on all kind of meshes (Ferziger and Peric, 2012) even if affected by some non physical oscillations in the solution of convection-dominated problems (Jasak, 1996). The second rule, the LimitedLinear interpolation, removes the oscillations which are tipically associated with second-order schemes (Sweby, 1984). The Total Variation Diminishing (TVD) criterion was introduced to characterise such oscillation-free flux-limited schemes (Harten, 1983).
OpenFOAM uses a general approach, the segregated method, to solve systems of coupled equations as Eqs.(4) (Patankar and Spalding, 1972;Issa, 1986). At each integration step, equations are solved separately for each variable with an iterative 100 procedure, in order to enforce the coupling, until a prescribed global residual level is achieved. Non-linear differential equations are linearized before discretisation and the non-linear terms are lagged (Jasak, 1996). The segregated solution approach consists of two main phases: the PISO (Pressure Implicit with Splitting of Operator) algorithm (Issa, 1986)  To solve an unsteady, incompressible Navier Stokes problem as in Eqs.
(3), OpenFOAM solvers merge the PISO and the SIMPLE procedures into the PIMPLE algorithm. Accordingly, an elliptic equation for the pressure field is derived by a combination of the continuity and momentum equations (Patankar, 1981). The integration procedure for each time step is described by the following (Jasak, 1996): 7. PIMPLE loop: repeat from step 2 for a prescribed number of iterations or until residual constraints are satisfied; As part of this work we extended the built in buoyantBoussinesqPimpleFoam solver, designed to solve heat transfer problems, to approach the infinite-Pr number problem. First, the overall structure of the solver has been modified by moving 120 the temperature equation, originally solved outside the PISO loop, within it. A flowchart comparing the structure of the two solvers is illustrated in Fig.3 . We compared the performances of the two solvers by monitoring the number of PIMPLE iterations required to converge, that is when residuals, for a given variable  Computation of temperature, originally done outside the PISO loop, has been moved into the loop.

Simulated Time Execution Time TOT PIMPLE Aver.PIMPLE
Teq. IN 10 −2 τκ 75 s 3.1 · 10 3 2.28 ±0.26 Teq. OUT 10 −2 τκ 239 s 1.3 · 10 4 9.56 ±19.53 Table 2. Execution time, number and average number of PIMPLE iterations for the two solvers for a simulated time of 10 −2 τκ (transient regime). the T eq.IN solver is about 3 time faster than T eq.OU T . Teq. OUT τκ 3.8 · 10 4 s 2.3 · 10 6 17.99±56.43 Table 3. Execution time, number and average number of PIMPLE iterations for the two solvers for a simulated time of τκ. The performance is comparable with that of the transient regime: the T eq.IN solver is about 5 time faster than T eq.OU T . Required loops for TeqIN (blue) are lower than those for TeqOU T (orange) in both cases. a graphical representation). We concluded that the T eq IN solver performance is clearly improved respect to that of T eq OU T for both the transient and the statistically stationary regime.
A more detailed representation of pressure and temperature residuals behaviour is shown in Fig.5 for a fixed time step. Each PIMPLE loop groups a fixed number of four points which are the initial residuals at every PISO loop. This is valid for the pressure field, which is always calculated into the PISO loop, but only for T eq IN temperature field. Looking at the pressure 135 panel, the first two PISO loops of each PIMPLE loop result the more important for convergence (cf. 2nd and 3rd point of each group). The same is valid for PIMPLE loops: only two cycles are required by T eq IN solver to converge against the six required by T eq OU T . This behaviour, in agreement with results of Tabs.(2,3), confirms that an average value of 2-3 PIMPLE loops is sufficient, for the T eq IN solver, to reach the prescribed convergence.
3.1 Implementation of the Pr → ∞ case.

140
The central goal of this work is to reproduce the infinite-Pr structure. As introduced, Eqs.(3) reduce to the elliptical problem: Hyperbolic and elliptic equations require different numerical strategies to be solved. However, to improve the convergence of the linear system solution, we want to maintain the original hyperbolic structure of the solver. This requires the presence of a time derivative, to preserve the diagonally dominant structure of the matrix resulting from the numerical discretization of the 145 momentum equation (Moukalled et al., 2016). With this aim, we introduced an explicit artificial time derivative on the left hand side of Eq. (6), designed so that it does not affect the result when convergence is reached. This derivative is computed within Is also clear that, in both cases, the first two PISO loops are the main steps to reach convergence.
the PISO loop as: where n indicates the n-th PIMPLE iteration and m the m-th PISO iteration. Once convergence is reached, |u m n − u m−1 n | < , 150 with the prescribed threshold, the auxiliary time derivative of Eq.(7) goes to zero recovering a steady problem.
It is worth noticing that in the Boussinesq approximation the density is held constant in the unsteady and in the advective terms but not in the gravity one, where it is assumed to depend on temperature (Oberbeck, 1879). The Boussinesq OpenFOAM solver defines two auxiliary variables in order to avoid errors in calculating the buoyant term on non-orthogonal and distorted meshes (OpenCFD, 2007): and it solves the momentum equation with respect to the effective pressure p rgh . Since in this work a Cartesian mesh is used, we explicitly rewrote the buoyancy term defining ρ as the sum of a mean value and a perturbation: where a linear relation between temperature and density is assumed with T 0 , ρ 0 are reference values.
Similarly to the spatial scheme, OpenFOAM provides several temporal discretization strategies. Among them we selected the following (using OpenFOAM definitions (Ferziger and Peric, 2012)): -Euler scheme (first order, implicit), defined as : -Backward scheme (second order, implicit): indicating with n the time step and with ∆t the time step width.
In the following section, all the previous features are tested, in order to validate the new solver.

Solver Benchmarking
Our solver has been tested considering both 2D and 3D cases. For the first, we simulated numerically the 2D results reported 170 by Whitehead et al. (2013) for increasing Pr. For the second, instead, we referred to Sotin (1999) in which the infinite-Pr configuration is studied with varying Ra number.

2D Benchmark
We replicated the same conditions of Whitehead et al. (2013): the xz-squared domain has fixed BCs for bottom and top temperature (T = 0 at z = 1, T = 1 at z = 0∧ 1 4 ≤ x ≤ 3 4 ) and zero normal gradient boundary condition at x = 0 and x = 1. Velocities

175
BCs are free slip on every side of the domain. A 128 × 128 grid has been used with the exception of the resolution study where different choices are shown. In Whitehead et al. (2013), the dynamics of a convective process is simulated in the case of a fixed Ra = 10 6 number and varying Pr ∈ [1, ∞). Fig.6 shows the evolution of temperature field for Pr = 1, 10, 10 2 , 10 3 , ∞. The process evolves through four different stages: plume formation, the plume head reaching the top of the domain, the head spreading and turning back to the bottom. These phases correspond, for Pr = 10 3 , at times t = 12.4 · 10 −4 , t = 27 · 10 −4 , t = 50 · 10 −4 180 and t = 1000 · 10 −4 (Fig.6, panel d) and all times are in unit of the non-dimensional diffusive time τ κ . The dynamics varies at different values of Pr and the new solver is able to reproduce all of them (Fig.6, panels a-e) even for longer times when the plume starts to oscillate (Fig.6, panel f).
All simulations have been run on the Laki cluster at INGV, Sezione di Pisa (2 × 8-cores Intel Xeon 2.40GHz).

185
As in Whitehead et al. (2013), we fixed Ra = 10 6 , but we focused only on Pr = 10 3 and Pr → ∞ cases, setting the residual threshold to 3·10 −7 . We performed our analysis comparing temperature and velocity profiles and streamfunction time evolution against the originals. The study has been done for three different solvers: the original solver buoyantBoussinesqPimpleFoam, the new solver for finite-Pr numbers convectiveFoam and for infinite-Pr numbers convectiveFoamInf. For each solver, a systematic analysis has been performed on both short and long term behaviour of the process. In both cases several choices of 190 resolution, algorithm or numerical schemes have been tested, in order to test the sensitivity of the results.
The main scope of this study is to explore how the resolution affects results, since higher resolution allow to resolve more detailed fluid structures. We led this analysis following the same approach adopted in Whitehead et al. (2013) where increasing 195 resolution provides different results.

Numerical Scheme
FV method requires interpolation rules to reconstruct face fields values from the volume ones. We tested the Linear and Lim-itedLinear spatial schemes implemented in OpenFOAM (OpenCFD, 2007). Concerning temporal domain, instead, we tested OpenFOAM's Euler and Backward schemes. We performed several attempts combining spatial and temporal interpolation 200 rules as follows: -Linear + Euler (LE); -LimitedLinear + Euler (LLE); -Linear + Backward (LB); -LimitedLinear + Backward (LLB); 205 We used the LB/LLB schemes in conjunction with convectiveFoam solver (CLB/CLLB for brevity) and with the original buoyantBoussinesqPimpleFoam solver (BLB/BLLB). The same naming convention has consistently extended to the LE and LLE schemes.

Algorithm
Algorithm analysis was performed by varying the number of inner PISO loops (which enforces incompressibility) and outer 210 PIMPLE loops (which enforces the coupling with the temperature equation). We named each configuration adopting the pXpY convention, where pX stands for X PISO loops and pY for Y PIMPLE loops. Since we observed that a p2p2 configuration is sufficient to reach the satisfactory residual, we focused on a limited number of PISO/PIMPLE iterations. In the following, only a selection of significant results is discussed, since some PISO/PIMPLE combinations, (p1p3 and p3p1), were computationally unstable.

Analysis of the Transient Regime
In this section we report, at different times temperature (Fig. 7) and velocity (Fig. 8) profiles related to the transient scenario comparing them with the originals. To compare the time-dependent behaviour, we also compare the value of the streamfuntion as a function of time (Fig. 9). Figs.7-9 show results from the resolution study performed with the momentum predictor, the LB scheme and the p2p2 algorithm. Results are in agreement with Whitehead et al. (2013), confirming that the 128 × 128 220 resolution is enough to reproduce all cases, while a 64 × 64 grid is insufficient, especially in reproducing the first time steps.
The results for the different tests on numerical schemes in Figs.10-12 indicate an overall scheme independence, with some discrepancies around t = 15·10 −4 , when the velocity and the temperature fields exhibit a highly dynamic behaviour. In this time steps, LB and LLB provide the closest match. The same is valid for the horizontal profile of vertical velocity at t = 1000 · 10 −4 ( Fig.11). For all simulations we used the momentum predictor, a 128 × 128 resolution and the p2p2 algorithm. Concerning algorithm, for all simulations, we used the momentum predictor, the LB scheme and a 128 × 128 resolution . The analysis of the results suggests p2p2 algorithm as a good compromise both in terms of computational time and accuracy. From the analysis reported above varying resolutions, algorithms and schemes, we selected a reference configuration aimed to develop further studies. In particular, the final configuration has: Linear spatial interpolation scheme, Backward temporal scheme, 128 × 128 spatial resolution and p2p2 algorithm. We remark this is not the more accurate setting but the best 230 compromise between simulation time and accuracy. A better performance of the p2p2 algorithm is evident.

Analysis of the long term regime
With the same configuration of Sec. 4.1.2, we then analyze the statistically stationary regime until two thermal times. Referring to the long term results of Whitehead et al. (2013), we reproduced the temporal evolution of the streamfunction and compared it with the original outcomes. As in the original work, we observed the emergence of two subsequent oscillating regimes with 235 different frequencies, whose, unlike those of Whitehead et al. (2013), didn't lock in a steady oscillation of fixed frequency.
To identify these oscillations we adopted the same approach of the original work, computing the heat flux (HF) and the streamfunction. The heat flux is defined as function of time: where T (x, z, t) is the temperature field, N the number of used level, the streamfunction is defined as: 240 where Ψ is the proper streamfunction, but, for simplicity we will refer to Sm(t) as streamfunction. Results for both the Sm(t) and HF are reported in Fig.16, in which a shifting related to the increasing resolution (right panel) is observable as in the original work. A complete characterization of the oscillating phenomena has been performed evaluating the amplitudes of the oscillations in the long term regime. The main difference of our work is the absence of a long term regime with locked   Table 4. Comparison of several parameters for the Pr = ∞, Ra = 10 6 case. With the same convention adopted in Whitehead et al. (2013)(1 refers to first oscillation, 2 to the second): Ams1 = time average of Sm, Af1 = amplitude peak to peak. T11 = period of the first oscillation, T12 = the doubled period after the previous oscillation (not observed in our simulations), Ams2 = the steady value between the two oscillations (or the time average when a non steady behaviour is present), Af2 = amplitude of oscillations peak to peak, T2 = period of the second oscillation.

250
Three-dimensional simulations are run following the work of Sotin (1999 Tmean = 0.500 Tmean = 0.561 Table 5. Columns 1-4: configuration for each 3D simulation. Columns 5-6: comparison of heat flux and mean temperature evaluated at the stationary state for the original data and convectiveFoamInf. For each value of the heat flux is reported the mean value and the relative fluctuation, for the temperature field, the value averaged over the volume. We reproduced some cases selected from the original work and validated our solver comparing the statistically stationary 255 behaviour of the flux at the top boundary Q t and the mean temperature in the volume of domain. As shown in columns 5-6 of Tab.5, the convectiveFoam results are close to those of Sotin (1999) and catch the increase of the heat flux and fluctua-Parallel efficiency has been estimated by evaluating the strong scaling, that is how the computation time varies with increasing number of cores, keeping the problem size fixed. We defined the speedup as the ratio: where T (1) and T (n) are the execution times on one and n cores, respectively. Results of the simulations are reported in

Further tests and portability 270
For completeness, additional test considering a relaxation parameter have been done. Results with relaxation exhibit a smoothing effect for both the temperature field and streamfunction with respect to those without relaxation, thus introducing a non acceptable deformation of both finite and infinite-Pr number results.
ConvectiveFoamInf, originally developed on OpenFOAM5, is available also on OpenFOAM7.
The relaxation and the portability tests, have been performed monitoring the same observables of the previous benchmark: 275 the vertical temperature profile and isolines of temperature and streamfunction. For all the described studies results are available under request, as for the complete results for the benchmark.
We built and tested convectiveFoam, a new solver dedicated to the simulation of infinite-Pr number fluids, using the Open-FOAM toolkit. The development of the new solver moved from very high, finite-Pr number to the infinite-Pr case. We called 280 this different variants convectiveFoam and convectiveFoamInf respectively. We tested the solver reliability by quantitative comparisons with the 2D simulations described by Whitehead et al. (2013) and the 3D simulations reported by Sotin (1999). A satisfactory agreement between our results and such cases was found. In particular, the analysis of the 2D cases indicates that: i) Moving the temperature equation into the PISO loop improved the performance, by reducing significantly the number of required PIMPLE loops to achieve convergence. The T eq IN solver is about 3-5 times faster than T eq OU T .

285
ii) Resolution higher than 128×128 did not considerably improve the results, with the exception of highly dynamic times for which a finer grid can give a more accurate description in all the domain. As in Whitehead et al. (2013), the bigger discrepancies are present for the 64 × 64 resolution.
iii) Concerning schemes, even if the LLB and LLE schemes produced good results, the LB choice, combining the spatial Linear and temporal Backward schemes, gave the best compromise between time simulation and matching with original results.
290 iv) Although we discarded some algorithms because of their inconsistency with original results, several combinations of PISO and PIMPLE loops gave good results in terms of reproducibility. The p2p2 algorithm was chosen as the optimal compromise between simulation cost and agreement with the previous results.
v) Further studies performed with and without relaxation parameter revealed that the use of relaxation gave results that were not consistent with the reference study. Further investigation is needed to further explore the origin of this behaviour.

295
For the 3D case, comparing our results with Sotin (1999) we found that: vi) The convectiveFoamInf solver was able to reproduce main results of Pr → ∞, isoviscous fluid simulations even if a more accurate study on resolution is necessary to avoid non-physical results.
vii) Strong scaling analysis showed a maximum speedup around 16000 cells/cores. viii) Portability was evaluated for the latest release OpenFOAM7. 300 We remark that in this work we analyzed idealized models. This step is just the first of a longer project: further developments should consider rheological aspect of the Earth mantle such as the temperature/pressure viscosity dependence, which still has unsolved numerical issues in 3D domains (Khaleque et al., 2015). Also, analysis of both Newtonian and non Newtonian viscosity has to be considered.