Melvin Creff

Jean-Luc Guermond

PRECONDITIONING OF THE PRESSURE SCHUR COMPLEMENT OF THE AUGMENTED LAGRANGIAN FORMULATION OF THE GENERALIZED STOKES PROBLEM *

Keywords: Augmented Lagrangian, Stokes problem, preconditioning, Saddle point problem AMS subject classifications. 65M60, 65M12, 65M15, 35L45, 35L65

Cette thèse, ces trois années de travail, furent une belle aventure grâce à la présence, la patience et la bienveillance de nombreuses personnes.

En premier lieu, je souhaite chaleureusement exprimer ma pleine gratitude envers mes deux directrices de thèse, Caroline Nore et Bérengère Dubrulle, pour leurs conseils, leurs aides et leurs encouragements, m'ayant permis de mener à bien les travaux de cette thèse. Je leur suis sincèrement reconnaissant de m'avoir suivi durant ces trois années et d'avoir contribué à ma formation de chercheur.

Je suis très reconnaissant envers Yannick Ponty et Nathanaël Schaeffer pour leur lecture et commentaires pertinents sur mon manuscrit ainsi qu'à Sébastien Galtier et Francky Luddens pour avoir participé au jury de ma thèse et pour les discussions intéressantes qui ont suivi.

Je remercie grandement Jean-Luc Guermond pour l'opportunité qu'il m'a offerte de venir travailler à TAMU et, par la même occasion, le reste de l'équipe du département de mathématiques à TAMU, qui m'a chaleureusement accueilli en terre étrangère.

Mes remerciements reviennent aussi aux nombreuses équipes techniques et administratives des différentes structures (LISN, CEA, IDRIS, SMEMAG, IUT Orsay, Paris-Saclay, TAMU) avec qui j'ai collaboré m'ayant permis de faire cette thèse dans la sérénité.

Je remercie par ailleurs l'équipe enseignant de l'IUT d'Orsay pour ses conseils sur l'enseignement ainsi que les élèves rencontrés durant ces trois années de monitorat.

Au final, je suis heureux d'avoir passé ces trois années, entouré d'un grand nombre de collègues (chercheurs, doctorants, stagiaires, étudiants) qui ont tous apporté à mon vécu par leurs discussions.

Il me reste à remercier mes parents, mes frères et ma soeur pour leur soutien mental sans qui je n'aurai pas fait les études m'ayant permis de faire cette thèse.

Résumé:

Dans cette thèse, l'objectif est de mieux comprendre l'effet dynamo pouvant apparaître dans les fluides électriquement conducteurs en étudiant les transferts d'énergie de la magnétohydrodynamique (MHD). Motivés par l'observation que l'effet dynamo est un phénomène de conversion d'énergie entre champs de vitesse et magnétique, nous avons développé une nouvelle méthode se basant sur l'étude locale (en temps, espace et échelle) des bilans énergétiques décrivant ainsi les différentes dissipations ainsi que les transferts inter-échelles. Notre approche se base sur un nouveau type de filtrage des champs pouvant être utilisé sur n'importe quel maillage de simulation, incluant les maillages non-structurés.

Pour simuler la MHD, nous utilisons le code SFEMaNS qui est une méthode hybride mêlant éléments finis et décomposition spectrale. Tout d'abord, nous présentons quelques améliorations de ce code : une meilleure répartition du maillage sur les processeurs et une étude portant sur l'implémentation d'une méthode de résolution des équations de Navier-Stokes à base de méthodes matricielles plutôt que celle de prédiction-correction actuellement utilisée. Ensuite, nous expliquons les nombreux outils informatiques implémentés pour pouvoir calculer, visualiser et étudier les transferts d'énergie.

Finalement, ces outils sont utilisés pour étudier l'expérience de von Kármán Sodium, dans les régimes de croissance et de saturation, qui possèdent deux types de dynamo en fonction du matériau constituant les turbines (acier ou fer doux). Bien que ces deux dynamos diffèrent peu du point de vue de la théorie de champ moyen (les champs de vitesses étant les mêmes), la localité de notre formalisme permet d'élucider l'origine de la différence entre les deux dynamos : pour les turbines en acier, l'effet dynamo provient d'un transfert venant du fluide dans le coeur du volume, alors que, pour les turbines en fer doux, l'effet dynamo provient exclusivement des turbines. Cette approche nous permet aussi de discuter des signes précurseurs à des phénomènes de dissipations anormales et de dynamo rapide, qui peuvent devenir importants dans la limite de fluide non visqueux. Title: How can magnetic fields be generated in fluids, stars and planets? Keywords: Magnetohydrodynamics, numerical analysis, fluid mechanics

Abstract:

In this thesis, the aim is to better understand the dynamo effect that takes place in conducting fluids by studying the various energy transfers in magnetohydrodynamics (MHD). Motivated by the observation that dynamo is a conversion mechanism between magnetic and kinetic energy, we develop a new approach to unravel the dynamo mechanism based on local (in space, scale, and time) energy budgets describing dissipation and scale-by-scale energy transfers. Our approach is based upon a new filtering approach that can be used effectively for any type of meshes,including unstructured ones.

In order to simulate the MHD, we use the SFEMaNS code which is a hybrid method using both a finite element and spectral decomposition. We first present some developments in SFEMaNS: a better repartition of the unstructured mesh on the processors and a study on solving the Navier-Stokes equations using preconditioning matrix methods instead of the current prediction-correction method. Then, we explain how various post-processing tools have been implemented in order to compute, visualize, and study the energy transfers.

Finally, these tools are used to study the von Kármán Sodium setup, showing dynamo action for two types of impellers (steel or soft iron) in the magnetic field growth and saturation phases. Although the two types of dynamo hardly differ from the mean-field theory point of view (the velocity fields are the same in both cases), the locality of our formalism allows us to trace the origin of the differences between these two types of dynamo: for steel impellers, the dynamo is due to the transfer of kinetic energy both in the bulk and in the vicinity of the impellers, whereas for soft iron impellers, the dynamo effect comes from the impellers. We also discuss possible signatures of precursors to anomalous dissipation and fast dynamo, that could become relevant in the inviscid limit.

-Introduction

. Context

The dynamo effect is a well-known process that generates a magnetic field from the movement of an electrical conductive matter. In the case of planets, it is a fluid (liquid iron in the outer core of the Earth) put in movement by different mechanisms (rotation, thermal convection, solutal convection, etc) while, in stars, it is ionized gas. In all cases, what is required is an electrical conductive matter, and an energy source maintaining the motion of the matter. The generation of a magnetic field has also been reproduced in the laboratory using liquid sodium set in motion by various mechanical devices [1,2,3]. The only experiment to have generated a magnetic field with an unconstrained flow is the von Kármán Sodium (VKS) experiment [4]. Previous experiments [1,2] have generated a magnetic field using mathematically predetermined flows.

The von Kármán Sodium (VKS) experiment is an impressive device resulting from a collaboration between the ENS Paris, ENS Lyon, CEA Saclay, and CEA Cadarache. The aim of this experiment is to create and maintain a magnetic field by transferring energy from a conductive fluid, thus experimentally reproducing the dynamo effect at play inside planets and stars. The fluid used is liquid sodium, kept at 120°C, and put in turbulent motion thanks to two counter-rotating impellers which are disks fitted with curved blades.

A sustained magnetic field was first measured in 2006 in this experiment. The magnetic field was generated only when the movement was driven by soft-iron impellers [4] (and not by copper or steel impellers)] and with a flow containing limited fluctuations [3]. This flow was created by the two impellers rotating in opposite directions with an unscooping motion for the curved blades. No dynamo was observed with the scooping direction of rotation. The last point proves that this experiment cannot be interpreted as a simple magnet rotating in a conductive fluid.

However, it does indicate that the fluctuations impede the dynamo mechanisms as already seen in numerical simulations for simpler geometry [5] or asymptotics [6]. The resulting field was on average dipolar, aligned on the symmetry axis, with an azimuthal component concentrated near the impellers. However, due to experimental constraints, it was not possible to measure fluid velocity or the non-axisymmetric variation of the magnetic field. For some parameters, magnetic field inversions were also observed, as is the case for the Earth's magnetic field.

The main difficulty in finding the conditions for which the dynamo effect occurs comes from the turbulent movement of the fluids in stars and planets. A turbulent flow corresponds to a disordered and rapidly varying veloc-ity field containing numerous vortices. This kind of flow appears when the viscous force (associated with the capacity of the flow to dissipate motion) is small compared to the inertial force of the system. This makes the flow span multiple spatial scales from the global one generated by external forces to the smallest vortices dissipating the kinetic energy. As the dynamo effect comes from this flow, it is difficult to pinpoint which events in the flow generate the magnetic energy. This can lead to a magnetic field generated by very rough velocity fields [7] through a mechanism called fast-dynamo [8,9] (which occurs when the exponential growth of the magnetic field is as fast as the advective time scale of the field). This can occur for very turbulent flows for which the stretching of the magnetic lines grows exponentially through vortices. The turbulent fluctuations may also be detrimental to the dynamo action by desorientation of the magnetic field lines [?].

To deal with this difficulty, one popular approach is to build simpler dynamo models using mean-field or stochastic theories [10,11]. These approaches revealed the importance of two aspects of the fluid's flow: the helicity (through the α-effect) and the differential rotation (through the Ω-effect). For the VKS experiment, various combinations of these effects were considered for the emergence of the magnetic field [12,13]. However, as these models only analyze the impact of the velocity field on the magnetic field, it is not possible to understand the differences experimentally observed when using different materials for the impellers. This requires other models for the interaction between the Ω-effect and the temporal or spatial variation of electric conductivity [14,15,16] or magnetic permeability [17,18,19,20]. It also leads to a few analyzes of increase in the α-effect by eddy collimation near the impellers [21,22]. Mean-field theory also makes it difficult to understand the initial phase of magnetic field growth, as it necessitates distinguishing between a mean and a fluctuating part of the magnetic and velocity fields. Lastly, as all the small-scale events are parametrized in coefficients for the α and Ω effects, observing events linked to possible fast-dynamo is not possible.

In this thesis, we propose a new method for understanding the dynamo mechanisms responsible for magnetic field generation in turbulent conducting fluids. The main idea of our approach is to analyze the local (in time and space) energy budgets by describing dissipations and energy transfers between scales and fields. In the case of turbulent fluids, this approach was proposed by Duchon and Robert [23]. This method is based on a filtering process that permits the separation of the different scales of the flow. Moreover, compared to other filtering techniques [24], it is possible to detect potential singularities thanks to the introduction of anomalous dissipation corresponding to energy transfers at smaller and smaller scales. This filtering process can be extended to magnetohydrodynamics [25,26]. However, in these references, variations in magnetic permeability or electrical conductivity are not taken into account, even though they are essential for the VKS experiment and our simulations. In addition, as these articles are specialized in the study of astrophysics (particularly solar winds), magnetic and kinetic energies are combined. Since we want to study energy transfer between the fluid and the magnetic field, we need to revise this approach by separating the two energies. Therefore, in this thesis, we extend this approach to magnetohydrodynamics with varying magnetic permeability and conductivity, providing a complete description of where all energy transfers (from fluid to magnetic field and from one scale to the other) take place. By combining some numerical simulations carried out previously with this new description of local energy transfers, it is now possible to determine where and when the dynamo effects occur.

. Magnetohydrodynamics equations

Firstly, the equations of magnetohydrodynamics (MHD) are used to model the interactions of a magnetic field with the velocity field of a conducting fluid. They originate from Maxwell's equations taken in the non-relativistic limit and coupled to the Navier-Stokes equations through the Lorentz force and Ohm's law. In this section, we will present how they are derived and some key quantities in relation to turbulence and the dynamo effect. In particular, the introduction of Reynolds numbers and selected anti-dynamo theorems are important for understanding the choice of configuration and its limits.

. Maxwell's equations

Maxwell's equations are a combination of four equations that predict the behaviour of the electric and magnetic fields in the presence of moving charged particles. In the laboratory frame, they are given by:

∂ t (µH) = -∇ × E, (1.1) ∇ × H = j + χu + ∂ t ϵE, (1.2) ∇ • ϵE = χ, (1.3) ∇ • µH = 0, (1.4)
where H is the magnetic field, E the electric field, j the electric current density, χ the volumic density of charges, u the fluid velocity, ϵ the electric permittivity and µ the magnetic permeability. Usually, these equations are presented using the constitutive relation between the electric displacement field D and the electric field and the relation between the magnetic flux density (or magnetic induction) b and the magnetic field: D = ϵE and b = µH.

When the fluid is neutral in charges (χ = 0) and goes at a non-relativistic speed we can use the quasi-static approximation. This means that the time needed for the charges to follow the variation of the electric field is negligible compared to the advective time at which the fluid flows. Hence, these equations can be simplified as:

∂ t (µH) = -∇ × E, (1.5)
∇ × H = j,

(1.6)

∇ • ϵE = 0, (1.7)
∇ • µH = 0.

(1.8)

In order to link the electric current density and the electric field, Ohm's law is used: j ′ = σE ′ in the frame of the fluid. In the laboratory frame, this relation reads, using the non-relativistic transformation for an electromagnetic field:

j = σ(E + u × b).
(1.9)

Using this relation to eliminate the electric field gives the equation for the evolution of the magnetic flux density:

∂ t b = ∇ × (u × b) -∇ × j σ .
(1.10) This is the equation for the time evolution of the magnetic flux density. It consists of two terms on the right-hand side: the first is the source of b and the only way to maintain it; the second is the encoding of the Joule effect, which is a loss of energy by the magnetic field through the electrical resistivity of the fluid.

. From the Navier-Stokes equations and the Lorentz force

The Navier-Stokes equations for an incompressible fluid, of kinematic viscosity ν and density ρ, under an external volumic force f is given by:

∂ t u + ∂ i (u i u) -ν∆u + ∇p ρ = f ρ , ∇ • u = 0 (1.11)
In MHD, the external force is the Lorentz force which, for a particle of charge q at velocity v under an electromagnetic field (E,H), is q(E + v × b).

Hence, here the total volumic force is, respectively on the positive and negative charges:

f + = χ + (E + u + × b), (1.12
)

f -= χ -(E + u -× b), (1.13)
where χ ± are the charge densities and u ± their velocity.

Using that, for a neutral fluid χ + = -χ -and that the electric current density is j = χ + u + + χ -u -, we get:

f = f + + f -= j × b.
(1. 14) Combining these equations with the equations obtained for the magnetic field gives the magnetohydrodynamics equations for a neutral conductive incompressible fluid in the quasi-static approximation:

∂ t u + ∂ i (u i u) -ν∆u + ∇p ρ = j × b ρ , (1.15)
∂ t b = ∇ × (u × b) -∇ × j σ , (1.16
)

∇ • u = 0, (1.17)
∇ • b = 0.

(1.18)

. Dynamo effect and Reynolds numbers

The dynamo effect is the amplification and maintenance of a magnetic field by the movement of an electrically neutral material. In our case, this material is a fluid, and the dynamo effect in a conductive fluid results from an instability. Indeed, b = 0 is always a solution for equation (1.16). Hence, for a non-trivial magnetic field to appear and be maintained, the non-linear term (∇ × (u × b)) needs to overcome the Joule dissipation (-∇ × j σ ), meaning that, for a system of characteristic size L and fluid velocity v:

∇ × (u × µH) ∇ × ((∇ × H)/σ) ∼ σµvL = R m > 1.
(

R m is called the magnetic Reynolds number and gives a minimum bound to when a magnetic field can be generated. More precisely, for a fixed velocity field, equation (1.16) is linear in b, and thus R m is the control parameter of this instability. In this framework with a prescribed-velocity dynamo (called kinematic dynamo), many analytical fields were derived (Ponomarenko dynamo [27], Roberts dynamo [28]) and experimentally reproduced [1,2]. However, the most interesting dynamo effect is for the non-linear coupling between the Navier-Stokes and the induction equations. For this coupled system, there is no theorem that states in which circumstances a flow will produce a magnetic fluid and, worst, there are many anti-dynamo theorems (conditions in which the dynamo cannot appear whatever the R m value) [10].

For example, a purely axisymmetric magnetic field can never be produced whatever the velocity field [29,30].

The second control parameter of interest is the (kinetic) Reynolds number for a turbulent fluid. This characterizes the degree of turbulence of a flow by comparing the nonlinear inertial term (∂ i (u i u)) with the viscous dissipation (ν∆u) in Eq. (1.15):

∂ i (u i u) ν∆u ∼ vL ν = R e .
(1.20)

The higher the kinetic Reynolds number, the smaller the structures can be in relation to the size of the setup.

Both numbers are linked by the magnetic Prandtl number which only depends on the physical properties of the fluid:

P m = R m R e = µσν.
(1.21)

Since the fluid used in the experiment is prescribed (i.e. liquid sodium), the only way to increase the chances of a magnetic field appearing is to increase the fluid's velocity. For liquid sodium at 120°C: ν = 7 10 -7 m 2 s -1 , σ = 9.6 10 6 S m -1 and µ = µ 0 = 4π 10 -7 H m -1 . This yields P m = 0.8 10 -5 and imposes that R e ≳ 10 5 at least for a magnetic field to appear. In the experiment, experimentalists needed to go as high as R e ≳ 6 10 6 in order to generate a magnetic field. This high Reynolds number implies that something may happen on a very small scale that may contribute to the emergence of the magnetic field.

. Numerical simulations

To study the magnetohydrodynamics equations and the dynamo mechanism, various numerical codes exist. Most of the time, each geometry requires its own specific numerical framework. In addition, depending on the objectives, a kinematic dynamo equation or complete MHD equations are calculated. When the focus is on the instability threshold, a given velocity field (from experimental results, analytical models, or Navier-Stokes solvers) is used to integrate the magnetic field equation (1.16) with a variable magnetic Reynolds number. The critical threshold R c m is defined when the solution b = 0 becomes unstable. When linear and saturation phases of magnetic field generation are the objective, the complete MHD equations (1.15-1. 16) must be computed.

For the Ponomarenko dynamo studied in [31] a 1D eigenvalue code for kinematic dynamo is developed. In [32] a parallelized pseudospectral code for the full MHD equations in a periodic box is used to better understand the Roberts dynamo.

Dealing with more complex flows, the authors in [33] study the dynamo driven by rotation and shear in the accretion disk of galaxies, while [34] investigate the dynamo action in rotating convective layers. The two codes used in these references employ periodic boundary conditions in at least two directions. For infinite cylinders (periodic boundary conditions in θ and z), there exist a few other codes [35,36,37], the last one being used to study Taylor-Couette dynamo. There are also a few codes dedicated to spheres: pseudo-spectral codes with the usual poloidal-toroidal decomposition of u and B [38,39,40,41], or finite volume algorithms [42], or finite elements [43] and even commercial codes like COMSOL used in [44]. However, to simulate finite domains, only a few exist: for a spheroid [45], for a torus [START_REF] Wu | On a dynamo driven by topographic precession[END_REF], and for cylinders [START_REF] Morales | Magnetohydrodynamically generated velocities in confined plasma[END_REF][START_REF] Giesecke | Generation of axisymmetric modes in cylindrical kinematic mean-field dynamos of vks type[END_REF].

In order to better understand the VKS experimental results and the dynamo mechanism in place, several numerical simulations were carried out. Simple periodic codes of the complete MHD equations with Cartesian geometry were used to estimate the critical magnetic Reynolds and led to other results such as Lorentz force hysteresis cycles and the non-linear dynamo condition [START_REF] Iskakov | An integro-differential formulation for magnetic induction in bounded domains: boundary element-finite volume method[END_REF][START_REF] Ponty | Dynamo action at low magnetic prandtl numbers: mean flow versus fully turbulent motions[END_REF][START_REF] Dubrulle | Bifurcations and dynamo action in a taylor-green flow[END_REF][START_REF] Ponty | Subcritical dynamo bifurcation in the taylor-green flow[END_REF]. A kinematic code in infinite cylinders was used to optimize the geometry of the impellers for the VKS experiment, and initial results were obtained on the structure of the magnetic field generated [START_REF] Krstulovic | Axial dipolar dynamo action in the taylor-green vortex[END_REF]. With a finite cylinder, the detrimental role of the lid layers behind the impeller disks was highlighted [START_REF] Ravelet | Bifurcations globales hydrodynamiques et magnetohydrodynamiques dans un ecoulement de von karman turbulent[END_REF][START_REF] Stefani | Ambivalent effects of added layers on steady kinematic dynamos in cylindrical geometry: application to the vks experiment[END_REF]. However, the magnetic field structure obtained was completely different from that observed in the experiment later on: the numerical magnetic field was dominated by a component perpendicular to the cylinder axis and another parallel to this axis (i.e. an azimuthal m = 1 mode compatible with Cowling's anti-dynamo theorem). Then, by introducing the α and Ω effects from the mean-field theory into the kinematic dynamo simulation, the averaged structure of the magnetic fields observed in the experiment was reproduced [START_REF] Laguerre | Effects of conductivity jumps in the envelope of a kinematic dynamo flow[END_REF][START_REF] Laguerre | Impact of impellers on the axisymmetric magnetic mode in the VKS2 dynamo experiment[END_REF]17]. However, it either necessitated very high values of α or a homogenous one which are both unrealistic. Finally, direct numerical simulations (DNS) and large eddy simulations (LES) that simulate the complete MHD equations and impeller counter-rotation have been developed [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF].

In this thesis, we will use new DNS computations performed with some parameters set in [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF]. These simulations use the SFEMaNS code [START_REF] Nore | Numerical simulation of the von Kármán sodium dynamo experiment[END_REF] which we have chosen for several reasons: it uses cylindrical coordinates to easily represent the cylinder, the impellers are well represented using a penalty method (see [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] p.5), data from various runs for the VKS experiment [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] are available and some tools to study turbulence have already been implemented [START_REF] Guermond | SFEMaNS documentation[END_REF]. To simulate a complete turbine rotation, it takes around ten hours of simulation on 2040 processors of the Jean-Zay supercomputer supplied by IDRIS. The geometry (see Fig. 1.1) is an approximated version of the experimental setup.

The inner volume for r ∈ [0, 1] corresponds to the liquid sodium while the outer volume for r ∈ [1, 1.6] is composed of two conductive solid layers: one with the same conductivity as the liquid for r ∈ [1, 1.4] (modeling the sodium layer) and one with 4.5 times the conductivity r ∈ [1.4, 1.6] (modeling the copper wall). The turbines have different magnetic permeability depending on the material used. In this setup, it is possible to numerically observe the generation of the magnetic field in very good agreement with experimental results [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] using DNS for small Reynolds numbers R e or LES [START_REF] Faller | Dissipation, cascades et singularités en turbulence[END_REF] for large R e (the maximum R e is 10 5 , corresponding to a magnetic Prandtl number P c m ≃ 7 10 -4 at threshold). In this thesis, we will use only the results obtained with R e = 1500 using DNS to study energy transfer down to the smallest scale without modeling. This Reynolds number is low compared to the real one, however, it still gives a good idea of the main ingredients for which the dynamo effect occurs. In particular, it has been shown that using impellers with high magnetic permeability does lower the critical magnetic Reynolds number at which a magnetic field is produced. As observed in the experiment, the threshold at which a magnetic field was generated depends on the material used for the impellers. However, the origin of the dynamo effect is still elusive.

. Outline

In order to analyze these numerical simulations and better understand the dynamo effect in the VKS experiment, in this thesis we will implement the local energy transfer approach derived from the MHD equations in the SFEMaNS code.

First, in chapter 2, we present the program used for the simulations as well as its framework. We also present the improvements that have been made to it and a study of a potential new way of simulating our equations.

Then, in chapter 3, we describe various tools developed in order to efficiently compute, analyze, and visualize the local energy transfers. These tools make it possible to have an efficient workflow.

Finally, in chapter 4, we derive the equations for the energy transfers and present the numerical setup. This leads to the analysis of all the different transfers for two illustrative cases using steel and soft iron impellers. 15

-SFEMaNS improvements

In order to simulate the magnetohydrodynamics equations and calculate the various energy transfers, we use the SFEMaNS code. This chapter is devoted to a brief description of this code and the various improvements that have been made to it. We will first describe the numerical framework, and then present the improvements made to mesh distribution on the processors. We will next examine the possibility of modifying the simulation method used. All the work in this chapter was carried out in collaboration with J.-L. Guermond during my stay at Texas A&M University from the 5th of January to the 3rd of February 2022 and from the 15th of July to the 3rd of December 2022.

The SFEMaNS program [START_REF] Nore | Numerical simulation of the von Kármán sodium dynamo experiment[END_REF] is a code written in Fortran-90 dedicated to axisymmetric systems. It uses the cylindrical coordinates and decomposes the fields on a finite element basis in the plane (r, z) and on Fourier modes along θ. This is done in order to highly parallelize the computation using the modules MPI for the Fourier modes and PETSC for a domain parallelization on the finite elements. In addition, for the sake of efficiency, Fast Fourier Transforms [START_REF] Cappanera | Turbulence in realistic geometries with moving boundaries: when simulations meet experiments[END_REF] are used to perform calculations involving any product of the fields.

For our setup, the simulated volume is presented in Fig. 1.1 where the turbines are given motion using a pseudo-penalty method (see [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] p.5). Most of the functions needed to do the basic calculations and parallelization were already implemented. Hence, the main focus is to understand how the different decompositions work in order to implement the new terms to calculate.

. Numerical approximations used in the SFEMaNS Code

In order to implement the energy transfers in SFEMaNS we need to transform the equations and other operations into the SFEMaNS framework. We will thus present the various numerical approximations done in SFEMaNS and how to use them. The goal is to be able to do various operations (integrals, derivatives, and solving equations) that are needed to calculate the energy transfers.

. Finite elements and Fourier decomposition

The objective of this section is to explain how to solve an equation using the finite element theory in the SFEMaNS code. The study will focus on the temperature field equation (which is already implemented) in order to present how different types of derivatives and other operations are approximated.

The finite element method approximates the space of functions by a subset called test functions. These test functions are defined on a meshgrid (see Fig. 2.1) of triangles in 2D such that, when the characteristic size of a mesh cell goes to zero, the test functions can generate all space of functions. For example, for what is called P 1 approximation (see the left panel of Fig. 2.3), the space function is approximated by the space of continuous piece-wise linear functions and the test functions are thus equal to 1 at one summit (or node) of the triangle and 0 at the two others. There exist higher polynomial approximations using as test functions Lagrange polynomials of higher order, meaning that there are more test functions per mesh cells.

In SFEMaNS, the mesh used is a combination of other meshes: the meshgrid used to define the flow field u is, most of the time, a subdomain of the one that defines the magnetic field. This is because the magnetic field also exists outside the area where the fluid is moving. It is thus possible to create more complex meshgrids by sticking together many smaller grids. The meshgrids for magnetic and velocity fields use P 2 polynomials while (see the central panel of Fig. 2.3) the meshgrid for pressure uses P 1 elements. This is because only some pairs of finite elements for pressure and velocity fields ensure convergence to proper solutions of the Navier-Stokes equations. P 3 finite elements (see the right panel of Fig. 2.3) are also implemented in the code but have not yet been used in simulations.

. Weak formulation of an equation

The time evolution of the temperature field T in a fluid is given by:

∂ t T + ∇ • (uT ) -K∆T = f, (2.1)
where u is the velocity field, f a thermal source and K the thermal diffusivity of the fluid.

The first step is to discretize the time for time-dependent functions. Using the notation f n = f (nτ ), the discretization uses the BDF2 method [START_REF] Loan | Computational Frameworks for the Fast Fourier Transform[END_REF]:

∂ t T -→ 1 2τ (3T n+1 -4T n + T n-1 ), (2.2)
∇ • (uT ) -→ ∇ • ((2u n -u n-1 )(2T n -T n-1 )), (2.3)
∆T -→ ∆T n+1 .

(2.4)

The next step is to multiply the equation by a test function ϕ and to integrate by parts to get rid of the highest derivatives (with Ω the domain where T is defined):

Ω 3 2τ ϕ + K ⃗ ∇ϕ • ⃗ ∇ T n+1 = (2.5) Ω 4T n -T n-1 2τ + f n+1 -∇ • (2u n -u n-1 )(2T n -T n-1 ) ϕ.
This is called the weak formulation of an equation (here, also discretized in time). This is done in order to project the equation on the subspace generated by the test functions. Moreover, getting rid of the highest derivative makes the calculation more precise.

For clarity purposes, let us rename the right-hand side as follows:

Ω 4T n -T n-1 2τ + f n+1 -∇ • ((2u n -u n-1 )(2T n -T n-1 )) ϕ = Ω F ϕ .

. Test functions and Fourier transform

Now that the weak formulation of the equation is obtained, the aim is to transform this equation into a linear problem using both Fourier and finite element decompositions of T .

First, the Fourier decomposition is done as such:

T (r, θ, z, t) = T 0 h (r, z, t)+ M m=1 T m,c h (r, z, t) cos(mθ)+T m,s h (r, z, t) sin(mθ), (2.6)
with the coefficients T m,c/s h living in the space generated by the test functions on the 2D meshgrid containing n p total points. The test functions used are also in Fourier space and thus defined for every point j ∈ [[1, n p ]] as:

ϕ m,c j (r, θ, z) = ψ j (r, z) cos(mθ) with 0 ⩽ m ⩽ M, (2.7
)

ϕ m,s j (r, θ, z) = ψ j (r, z) sin(mθ) with 1 ⩽ m ⩽ M, (2.8)
where ψ i (r, z) are the test functions defined on the triangles of the 2D meshgrid.

T can then be written as:

T (r, θ, z, t) = np j=1 T 0 j (t)ψ j (r, z)+ (2.9) M m=1 np j=1 T m,c j (t)ψ j (r, z) cos(mθ) + T m,s j (t)ψ j (r, z) sin(mθ).
Hence, the left-hand side of equation (2.5) reads for a cosinus test function

ϕ m ′ ,c i (after integration over θ), m ′ ∈ [[0, M ]] and i ∈ [[1, n p ]]: np j=1 T n+1,0 j r,z r 3 2τ ψ j ψ i 2πδ 0,m ′ + (2.10) np j=1 M m=1 r,z r T n+1,m,c j δ m,m ′ π 3 2τ ψ j ψ i + K(∂ r ψ j ∂ r ψ i + ∂ z ψ j ∂ z ψ i + m 2 r 2 ψ j ψ i ) ,
and for the right-hand side:

np j=1 r,z rF 0 j ψ j ψ i 2πδ 0,m ′ + np j=1 M m=1 r,z rF m,c j δ m,m ′ πψ j ψ i .
(2.11)

Here r,z f (r, z) means f (r, z) drdz.

Hence solving (2.5) for T n+1 becomes solving for each Fourier mode m and each test function i (and for the cosine part):

A m,c ij T n+1,m,c j = B m,c i , (2.12)
with:

A m,c ij = r,z r 3 2τ ψ j ψ i + K(∂ r ψ j ∂ r ψ i + ∂ z ψ j ∂ z ψ i + m 2 r 2 ψ j ψ i ) , (2.13)
B m,c i = np j=1 r,z rF m,c j ψ j ψ i . (2.14)
It is the same linear system for the sinus part in the case of the temperature field. Now the Fourier modes m are uncoupled and thus can be solved simultaneously in parallel. Moreover, as this is a linear system requiring only a matrix inversion, powerful tools of linear algebra can be used, such as preconditioning and LU factorization. However, in order to compute the coefficients of A and B, some integrals need to be calculated.

. Evaluation of integrals

A 2D mesgrid Ω r,z is composed of n m cells such that Ω r,z = nm n=1 K n with a total number of points n p . For a function G defined in the (r, z) plane for the {ψ i } i∈[ [1,np]] test functions, its integral is first cut over all n m mesh cells :

Ωr,z G(r, z) drdz = nm n=1 Kn G(r, z) drdz. (2.15)
Then each cell is transformed into a reference cell K with J n the associated Jacobian and ψi the reference local test function:

Ωr,z G(r, z) drdz = nm n=1 K G(r, z)J n (r, z) drdz.
(2.16) Thanks to this transformation, all different test functions ψ become generic functions. Hence, it is not needed to keep in memory all test functions but only the ones on the reference cell K and the Jacobian associated with each cell.

Then, using Gauss quadrature over n G Gauss points (meaning approximating the integral by a sum over well-chosen and weighted points), we obtain:

Ωr,z G(r, z) drdz = nm n=1 n G l=1 G(l)J n (l).
(2.17)

Using the decomposition of G over the generic test functions then gives, with {n glob } n ⊂ [[1, n p ]] the global indices of the test functions at the nodes of the cell n:

Ωr,z G(r, z) drdz = nm n=1 n G l=1 J n (l) i∈{n glob }n G i ψi (l).
(2.18) For example, in order to evaluate A m,c ij , we use the sum:

A m,c ij = nm n=1 n G l=1 J n (l)r(l) (2.19) 3 2τ ψj (l) ψi (l) + K ∂ r ψj (l)∂ r ψi (l) + ∂ z ψj (l)∂ z ψi (l) + m 2 r(l) 2 ψj (l) ψi (l) ,
where the radius at the Gauss point l is r(l) = k∈{n W }n r k ψk (l) with r k the radius at the node k. The functions { ψi , ∂ r ψi , ∂ z ψi } are known as well as the Jacobian J n (l) when given the meshgrid.

. Derivatives

The last types of operation that need to be coded in the SFEMaNS code are derivatives. For efficiency and accuracy purposes, the derivatives are calculated on the Gauss points. In addition, the Fourier mode decomposition of the fields into sine and cosine parts must be taken into account. For example, one operation needed is to calculate j = ∇ × H, which decomposition in components (r, θ, z) gives:

j r = 1 r ∂ θ H z -∂ z H θ , (2.20
)

j θ = ∂ z H r -∂ r H z , (2.21
)

j z = 1 r ∂ r (rH θ ) - 1 r ∂ θ H r . (2.22)
Using the Fourier and finite element decompositions, this gives for cosine mode m at the Gauss point l:

j m,c r (l) = m r(l) i H m,s z,i ψ i (l) - i H m,c θ,i ∂ z ψ i (l), (2.23)
j m,c θ (l) = i H m,c r,i ∂ z ψ i (l) - i H m,c z,i ∂ r ψ i (l), (2.24)
j m,c z (l) = 1 r(l) i H m,c θ,i ψ i (l) + i H m,c θ,i ∂ r ψ i (l) - m r(l) i H m,s r,i ψ i (l), (2.25) j m,s r (l) = - m r(l) i H m,c z,i ψ i (l) - i H m,s θ,i ∂ z ψ i (l), (2.26)
j m,s θ (l) = i H m,s r,i ∂ z ψ i (l) - i H m,s z,i ∂ r ψ i (l), (2.27)
j m,s z (l) = 1 r(l) i H m,s θ,i ψ i (l) + i H m,s θ,i ∂ r ψ i (l) + m r(l) i H m,c r,i ψ i (l), (2.28)
where {ψ i , ∂ r ψ i , ∂ z ψ i } are already implemented using the previous Jacobian to go to the reference cell and using the reference test functions ψi .

That concludes a quick overview of the SFEMaNS framework and its various operations that need to be hard-coded when developing new terms in this code. Indeed, other simpler operations such as the scalar product and the cross product are already fully implemented.

. New mesh distribution on processors

This section is dedicated to improving the distribution of the meshgrid across processors. We need to do this in order to simulate fine meshes without being limited by the per-processor memory. This makes it possible to simulate higher Reynolds numbers. The problem was that each processor stored a copy of the entire grid in its memory. As each processor has only 4 to 8 GB of memory (depending on the computing cluster), it was impossible to obtain 2D meshgrids of more than 500,000 points, regardless of the number of processors used. This severely limited the highest Reynolds number achievable in 2D simulations and was very inefficient in terms of memory. It is now possible to simulate finer meshes by increasing the number of processors. This is very important for simulating magnetohydrodynamics on fine meshes, as several meshes are used (one for the velocity field, one for the pressure field, and one for the magnetic field).

To solve this problem, we have implemented a new way for each processor to save the meshgrid so that they only save in memory the cells that are allocated to them. This new implementation made it possible to have around 300,000 meshgrid points per processor (with 4Go of memory each) for 2D Navier-Stokes problems. Using multiple processors we were able to start 2D Navier-Stokes problems with hundreds of millions of points.

Two main difficulties were encountered. Firstly, the new implementation had to respect the already established convention and structure of the meshgrid. Secondly, due to parallelization, some processors needed information about bordering cells that are allocated to other processors.

As there are three levels of indexing (locally at the cell level, at the processor level, and at the global level), all the local numbering and indexation will refer to the ones at the processor level. I will also use directly the names of the Fortran tables and functions given in the code. Also, the notations used in this section are not correlated to the ones used before.

. Cell level convention

The convention that needs to be respected for the calculation to proceed accurately (or even at all) is schematized in Fig. 2.2 and 2.3. The vertex of the triangle are number for i ∈ [ [1,3]] will the other nodes are in [ [4,10]].

The first convention (Fig. 2.2) is that, for a cell, its vertices, edges and neighbours are cell-numbered such as, for a vertex with cell-numbering i ∈ [ [1,3]]:

• the edge at the opposite of i is cell-numbered i too,

• the cell neighbouring the opposite edge is also cell-numbered i.

The easiest way to memorize it is that all objects with the same cell number are 'aligned'. The second convention is that, for a cell, its finite element points are cellnumbered according to Fig. 2.3. This means that, for all P n finite elements, the points cell-numbered i ∈ [ [1,3]] correspond to the respective vertices and the edges are oriented from the lowest cell index of the vertex to the highest one.

The other points follow the convention in Fig. 2.3. As an example, in P 3 finite element points, the points 4 and 5 are on the cell-numbered edge 1 opposite to the cell-numbered vertex 1 such that 4 is closer to the cell-numbered vertex 2. These P 3 finite elements are the reason why we want the edges to be oriented according to a specific convention. For neighbouring cells to have the points on the edges in the same order, the orientation of the edges must be the same from one neighbor to the next. This means that, for an edge between two cells m 1 and m 2 with cell-numbering of the edge's vertices

(v 1,m 1 , v 2,m 1 ) in the cell m 1 and (v 1,m 2 , v 2,m 2 ) in m 2 : if v 1,m 1 < v 2,m 1 then v 1,m 2 < v 2,m 2 .

. Reading the global mesh

The first step is to read an already built P 1 mesh. This is done by using the function load_mesh_free_format_global_p1. For the mesh: These arrays are build so that, for i ∈ [ [1,3]], mesh_glob%jj(i,m) is the point opposite to the edge mesh_glob%jce(i,m) and the cell m has mesh_glob%neigh(i,m)

• mesh_glob%me = M • mesh_glob%np = mesh_glob%dom_np = N • mesh_glob%medge = M e For each cell m ∈ [[1, M ]] of
for neighbours along the edge mesh_glob%jce(i,m) (Fig. 2.2). These arrays are also built so that each edge is oriented in the same way. This is done so that the conventions discussed in section 2. It should be noted that there is no reason for the value of the indices v i and e i (for i ∈ [ [1,3]]) to be in any particular order. This property will be true for any other mesh that will be discussed in all the following sections.

For a point v ∈ [[1, N ]]:

• mesh_glob%rr(1,v) gives its r coordinate,

• mesh_glob%rr (2,v) gives its z coordinate.

In addition to this mesh structure, the external boundaries of the mesh must also be described in order to correctly take into account the different boundary conditions. They are thus separated in N interf ace interfaces making it possible to specify various boundary conditions (e.g. Dirichlet or Neumann or Robin conditions) on different borders of the mesh. The total number of edges at these interfaces is M es and is saved in mesh_glob%mes. For each edge

m s ∈ [[1, M es ]] : • mesh%sides(m s ) gives the number ∈ [[1, N interf ace ]] of the interface, • mesh%neighs(m s ) gives the associated cell ∈ [[1, M ]],
• mesh%jjs(m s ) gives the points ∈ [ [1, N ]] on the edge m s .

. Distribution of the P 1 global mesh on processors

Now mesh_glob will be distributed on N p processors, each processor having then its own mesh_p1_loc. This is done using the function reorder_mesh which does two things: first, it reorders the cells, edges, and vertices in a specific order; then, it distributes to each processor its part of the global mesh by calling the function create_local_mesh_with_extra_layer.

Attribution of cells and reordering

Using METIS_PartGraphRecursive, each cell m ∈ [[1, M ]] is attributed to a pro- cessor p c ∈ [[1, N p ]] (p c
is called the rank of the processor) (see example in Fig. 2.4). The points and edges are then allocated to the lowest-rank processor (p c ) whose cells contain this point or edge. This means that some processors can end up with no vertex (see processor 3 in Fig. 2.5, left panel). The same is done at the interfaces (i.e. on the external boundaries of the global mesh). The distributed mesh is called mesh_p1_loc. These numbers are saved by all processors in, for

i ∈ [[1, N p ]] : • mesh_p1_loc%domcell(i) = M (i), • mesh_p1_loc%domnp(i) = N (i), • mesh_p1_loc%domedge(i) = M e (i),
and the cumulative associated quantities are saved in:

• mesh_p1_loc%discell(i + 1) = 1 + i k=1 M (k), • mesh_p1_loc%disp(i + 1) = 1 + i k=1 N (k), • mesh_p1_loc%disedge(i + 1) = 1 + i k=1 M e (k).
By convention, mesh_p1_loc%discell(1) = mesh_p1_loc%disp(1) = mesh_p1_loc%disedge(1) = 1. These quantities are designed so that, for a point, cell, or edge, it is easy to know which processor it has been assigned to.

Each processor has a local indexation for its attributed cells, vertices, and edges, denoted by a subscript l. There is also a global indexation for all processors denoted by a subscript g. To facilitate the switch between global and local numbering of vertices, edges, and cells, they are reorganized as follows:

• for a cell with global number

m g ∈ [[1, M ]] attributed to processor p c ∈ [[1, N p ]] and whose processor number in p c is m l ∈ [[1, M (p c )]], the rela- tion reads: m g = m l + mesh_p1_loc%discell(p c ) -1,
• for a point with global number

v g ∈ [[1, N ]] attributed to processor p c ∈ [[1, N p ]] and whose processor number in p c is v l ∈ [[1, N (p c )]], the relation reads: n g = n l + mesh_p1_loc%disp(p c ) -1,
• for an edge with global number In other words, for cells, the M (p c = 1) first cells for the global index ∈ [ [1, M ]] are attributed to the first processor, the next M (p c = 2) cells in the global index are attributed to the second processor, and so on. The same goes for points and vertices. This makes it possible to define the function p(.), which uses the global index of an object to return the associated processor. For example, for a point of global index

e g ∈ [[1, M e ]] attributed to processor p c ∈ [[1, N p ]]
v g ∈ [[1, N ]], p(v g ) is such that p(vg)-1 k=1 N (k) < v g ≤ p(vg) k=1 N (k).
This reordering will change the values in mesh_glob%jj, mesh_glob%jce, mesh_glob%neigh, mesh_glob%rr, mesh_glob%jjs, mesh_glob%neighs and mesh_glob%sides as well as their ordering but it makes the actual distribution on processors much easier. As only the cell's global indices are changed and not their order at the cell numbering level, all the conventions reported in section 2.2.1 are still respected.

Distribution and dealing with external boundaries

Now that each cell, edge, and point is associated with a processor, the objects (%jj, %jce, %neigh, and so on) need to be created so that each processor knows the structure of its local mesh. In addition, other objects are created so that each p processor knows just enough about the cells that are not assigned to it, but which contain at least one point or edge assigned to it. In Fig. 2.4, for the processor 0 (dark blue), this would correspond to all the cyan, red, and yellow cells (associated with processors 1, 2, and 3) that contain one of the big blue dots. The aim is to add to each processor only the necessary information about its boundaries, in order to save memory.

Because of the borders between processors, for some processors:

• some of their cells have neighbours that are not attributed to those processors,

• some vertices of their cells are not attributed to those processors,

• some edges of their cells are not attributed to those processors.

However, to build the matrices needed to solve any problem (see section 2.1.2), each processor needs the following information: for each point v attributed to a processor, this processor needs to be able to construct the global index of all points having a cell in common with v. There are two scenarios in which this processor does not own all the points connected to its attributed point v:

• some points or edges of an attributed cell are not attributed to this processor,

• a cell of an attributed point or edge is not attributed to this processor.

The first scenario happens for a processor if one of its cells' neighbours is attributed to a lower-rank processor and the second scenario is just the other way around. As an example, in Fig. 2.4, the first scenario occurs for all the cyan cells (attributed to processor 1) that border the dark blue cells (attributed to processor 0): the processor 1 needs the global index of the big dark blue dots touching any cyan cell. In the same figure, the second scenario arises in the reverse situation, the processor 0 needs to know the information of all cells containing one of the big blue dots for a total of 12 cyan, 2 yellow and 12 red cells. Similarly, the processor 2 requires the information of only 10 red cells (the ones in contact with a big yellow dot).

The first scenario is the easiest to deal with: it is only needed to extend the number of vertices saved by each processor by N extra (p c ) vertices. Those new vertices are in an attributed cell of p c but not attributed directly to p c . The same procedure is used for the additional M e extra (p c ) edges that are in an attributed cell of p c but not attributed directly to p c .

In the second scenario, the processor must save the entire structure of the neighbouring cell. The structure of these M extra (p c ) cells will be saved in a new table. For example, in Fig. 2.5 (left), processor 1 requires 3 extra cells (2,3,4), processor 2 requires 2 extra cells (3,4), processor 3 requires 1 extra cell (4) and processor 0 none. For a processor p c ∈ [[1, N p ]], the local numbering goes:

• for a cell m l ∈ [[1, M (p c )]],
• for a point

v l ∈ [[1, N (p c )]
] if it is attributed to the processor, and

v l ∈ [[N (p c ) + 1, N (p c ) + N extra (p c )]] otherwise,
• for an edge e l ∈ [[1, M e (p c )]] if it is attributed to the processor,

• for an edge e l ∈ [[1, M e extra (p c )]] if it is not attributed to the processor.

This defines new useful constants:

• mesh_p1_loc%np = N (p c ) + N extra (p c ),
• mesh_p1_loc%mextra = M extra (p c ),

• mesh_p1_loc%medges = M e extra (p c ).

As an example, the values for Fig. 2.5 are, with p the processor number :

Unrefined (left panel) Refined (right panel)

p 1 2 3 4 1 2 3 4 M extra (p) 3 2 1 0 5 4 3 0 N extra (p) 0 1 3 2 0 1 5 3 M e extra (p) 0 0 2 1 0 0 4 2
Now that the local and global numbering are fixed, here are all the necessary structures: -mesh_p1_loc%jees(e l ) which is the edge's global index, -mesh_p1_loc%jecs(e l ) which is the local index of the cell containing that edge,

• for a cell at the borders with a higher rank processor and not attributed The tables %jj, %jce, %neigh, %extra_jj, and %extra_jce respect the cell conventions (see section 2.2.1) as they are extracted from the global mesh corresponding tables. Only the output indices are modified, not the actual structure of the cell. Now that the mesh is locally distributed, mesh_global is deallocated. The objective of ensuring that each processor retains only the necessary mesh information is achieved so that no new tables are created. All that remains is to refine the mesh and create the P 2 and P 3 vertices from the P 1 points.

to p c , cell m l ∈ [[1, M extra (p c )]] - mesh_p1_loc%extra_jj

. Refinement of the distributed P 1 mesh

Once the P 1 mesh has been distributed, it can be refined if necessary. The refinement process is carried out by subdividing any cell into four (see Fig. 2.6). The average cell size is thus halved. The refined mesh is called mesh_p1_loc_r. For the refined mesh, the numbers of cells, edges, and pointss attributed to a processor are given by:

v 1 v 3 v 2 m e 3 e 1 e 2 v 1 v 3 v 2 e 3 v n 3 e E 3 e 1 v n 1 e E 1
• mesh_p1_loc_r%dom_me = 4 * mesh_p1_loc%dom_me, • mesh_p1_loc_r%dom_np = mesh_p1_loc%dom_np + mesh_p1_loc%medge, • mesh_p1_loc_r%medge = 2 * mesh_p1_loc%medge + 3 * mesh_p1_loc%dom_me.
In other words, there are four times as many cells and one more point per edge. For the edges, every edge is divided into two new ones and there are three new ones to make the center cell (see Fig. 2

.6).

The updated cumulative quantities and per processor numbers of objects are saved in the same way as before in mesh_p1_loc_r%disxyz and mesh_p1_loc_r%domxyz (where 'xyz' should be replaced by cell, edge or np).

The notation for these numbers (M , M e , and so on) will include an exponent r if it is the corresponding number for the refined mesh. As an example, the formulas for the new numbers of cells, edges, and points attributed to a processor p c ∈ [[1, N p ]] will be written as:

• M r (p c ) = 4 * M (p c ), • M r e (p c ) = 2 * M e (p c ) + 3 * M (p c ) , • N r (p c ) = N (p c ) + 2 * M e (p c ).
Now the actual structure of the new cells needs to be created. Processor refinement and numbering are based on the following diagrams (see Fig. 2.6). An easy way to understand it is that, similarly to Fig. 2.2, every object with the same subscript number is aligned e.g. v 1 , e 1 , m 1 , e I 1 and e E 1 . For the cells, it means that the center cell m keeps the processor index of the unrefined one. For edges, the refined edge closest to the lowest cell index retains the processor index of the unrefined edge. And, for the points, the ones already existing keep their local index. For i ∈ {1, 2, 3} the cell index of in the unrefined cell, the local indexation for the others follows:

For the bulk

• for the cells:

-m r i l = M (p c ) + m l -1 + i
• for the edges (we call 'outside edge' the edges (e E i , e i ) and 'inside edge' the edge e I i as displayed in Fig. 2.6):

for the new outside edges, e r • for the points:

-v r n i l = N (p c ) + e i l
The cell-indices of the new cells are such that, for the center, the points are ordered as (v n 1 , v n 2 , v n 3 ) and, for the corner, the first point is an old one and the two others are in the order of the center cell. In Fig. 2.6, this can be seen by the orientation of the new edges. This is done in order to respect the convention displayed in Fig. 2.3 for P 1 finite elements.

The global indices are then constructed by adding

pc-1 k=1 M r (k) to the cells, pc-1 k=1 M r e (k)
to the edges and pc-1 k=1 N r (k) to the vertices. This means that, if the processor rank is greater than one, the ones whose local index did not change still have their global index shifted by the number of corresponding objects created in lower-rank processors. This is done in order to keep an easy conversion from local index to global index (see the end of section 2.2.3). This new indexation allows the construction of mesh_p1_loc_r%jj, mesh_p1_loc_r%jce, mesh_p1_loc_r%i_d, mesh_p1_loc_r%rr and mesh_p1_loc_r%loc_to_glob but only for cells, vertices, and edges attributed to p c . The position (r, z) of the new vertices v r i is simply the middle of the unrefined edge. For example, for v r mesh_p1_loc_r%rr(:,v r 1 l )= (mesh_p1_loc%rr(:,v 2 )+mesh_p1_loc%rr(:,v 3 ))/2. These structures now only need the indexation for cells and edges bordering a processor with a lower rank as they would not be attributed to p c (which will be presented in the following). Before that, one last thing can be built: the neighbours of the new cells. This is done by following again Fig. 2.6. The neighbours of the center cell m l are simply (m l 1 , m l 2 , m l 3 ) (this order does respect the cell level convention of neighbours, see Fig. 2.2). For a corner cell i ∈ {1, 2, 3}, the first neighbour is always m l and the others are determined by using mesh_p1_loc%neigh(:,m l ).

The two simplest cases are: if one edge of the refined cell belongs to an external boundary, the value is 0 and, if the neighbour of the unrefined cell m l is not attributed to the same processor, the value is -1. Lastly, if the cell is also attributed to the processor p c , it is necessary to reconstruct the cell index of this neighbour. With m n l the local index of the unrefined neighbours border- ing a new cell m i l and i n the cell index of the vertex shared by m i l and m n l , the new neighbour of m i l would have a local index of

M (p c ) + m n l -1 + i n .
It is the same formula as above but from the point of view of another cell. This is used to construct all of mesh_p1_loc_r%neigh (if the neighbours cell is not attributed to the processor, it is always assigned the value -1, so there is nothing special to do in these cases).

At the border with a lower-rank processor

The only indexation that is left is for the N extra (p c ) points and the M e extra edges that are in one of the cells attributed to p c and that are on the border with a lower-rank processor. Because of the refinement, new points and edges will be created on this border. The formulas for the new numbers of edges and points not attributed to the processor p c but belonging to one of the cells attributed to p c are:

• N r extra (p c ) = N extra + M e extra , • M r e extra = 2 * M e extra .
Each edge not attributed to p c will be divided in two and a new point will be created. They all will not be attributed to p c as the edge is already attributed to another processor.

The local index in p c of all the N extra unrefined points needs to be shifted by N r (p c ) -N (p c ) so that, for the refined mesh, we still get that, for v l ∈ [[1, N r (p c )]], it is attributed to the processor and, for

v l ∈ [[N r (p c ) + 1, N r (p c ) + N r extra (p c )]], it is not. In contrast, the M r
e extra extra edges do not need this shift.

The local index in p c of newly created extra objects simply follows the existing ones, in no particular order, as it has no impact.

However, it is very important that the global indices of the extra points and edges are well constructed. For an already existing point whose processor numbering in the unrefined mesh is

v l ∈ [[N (p c ) + 1, N (p c ) + N extra (p c )]],
it is possible to get its global index using v g =mesh_p1_loc%loc_to_glob(v l ). With that it is possible to find which processor that vertex is attributed to p(v g ) and, as for the bulk, the new global index only needs to be shifted by the number of vertices created in lower-rank processors: v r g = v g + p(vg)-1 k=1

N r (k). For an unrefined edge with local numbering e l (p c ) ∈ [[1, M e extra (p c )]] and global numbering e g =mesh_p1_loc%jees(e l ), two actions need to be done for the refinement: divide it by two and create a new vertex in the middle. Similarly to the vertex, it is possible to get the edge processor p(e g ) and thus get the edge's local numbering in this processor on the unrefined mesh using e l (p(e g )) = e g -p(eg)-1 k=1

M r e (k).

Hence, constructing the new point is easy. Using the same formula as for the bulk but from the point of view of another processor, the global index of this new vertex is v r newg = N (p(e g )) + e l (p(e g )) +

p(eg)-1 k=1

N r (k). The sum at the end is the shift on the global index due to the vertices on lower-rank processors.

The same is applicable for the refined edges. As for the bulk, the one closest to the lowest cell numbering vertex will have its global index shifted e r g = e l (p(e g )) +

p(eg)-1 k=1

M r e (k) while the other will use e r newg = M e (p(e g )) + e l (p(e g ))+ p(eg)-1 k=1

M r e (k). This makes it possible to finish building mesh_p1_loc_r%jees and mesh_p1_loc_r%jecs. All that remains is to refine the extra cells and the external boundaries.

At the border with an higher-rank processor

The refinement of the M extra (p c ) extra cells is done in the same way than for the bulk. However, there are two difficulties. Firstly, a recount of the extra cells is needed to keep it to a minimum to preserve memory usage. Secondly, as cells, points, and edges can belong to different processors, the local index in their attributed processor needs to be reconstructed for all objects similarly to what was done in the previous section 2.2.4. The only easy aspect is that the position of points does not need to be calculated.

The counting of the M r extra (p c ) new extra cells is similar to the global mesh distribution on the processors (see section 2.2.3). We count the number of refined cells not attributed to the processor p c but in contact with one of the refined vertices attributed to p c . For the refined mesh in Fig. For an extra unrefined cell of global index v g with edges global indices e ig and points global indices v ig , their local index in their attributed processor is given by:

• v l = v g - p(vg)-1 k=1 M r (k),
• for i ∈ {1, 2, 3}, e i l = e ig -p(e ig )-1 k=1

M r e (k),

• for i ∈ {1, 2, 3}, v i l = v ig - p(v ig )-1 k=1 N r (k).
This is done by simply using the conversion from global index to local index (see the end of section 2.2.3). Now the global index for the refined objects is built following the formulas for the bulk (see section 2.2.4) and adding the shift necessary for the global index. For i ∈ {1, 2, 3}, using the same notation as in Fig. 2.6: • for the refined extra cells:

-if it is a center cell, m r g = m l + p(mg)-1 k=1 M r (k), -if it is a corner cell, m r ig = M (p(m g ))+m l -1+i+ p(mg)-1 k=1 M r (k),
• for the edges of one refined cell (we call 'outside edge' the edges (e E i , e i ) and 'inside edge' the edge e I i as displayed in Fig. M r e (k),

• for the points:

-if it is an existing point, v r ig = v i l + p(v ig )-1 k=1 N r (k), -if it is a new point on the edge e ig , v r n ig = N (p(e ig ))+e i l + p(e ig )-1 k=1 N r (k).
For new objects, the processor p(.) used is in the above formulas is the one they will be attributed to. As an example, for the outside edges, the processor of the unrefined edges is used whereas, for inside edges, the processor of the unrefined cell is used. These new global indices are what is necessary and sufficient to build mesh_p1_loc_r%extra_jj, mesh_p1_loc_r%extra_jce, and mesh_p1_loc_r%extra_jcc.

External boundaries

The only structures left to implement are mesh_p1_loc_r%jjs(:,:), mesh_p1_loc_r%sides(:), and mesh_p1_loc_r%neighs(:). They correspond to the structure of the outside edge of the total mesh and are used to impose boundary conditions. As these edges are in contact with only one cell, their refinement is simple. Hence, the number of refined edges at the external boundaries is M r es (p c ) = M es (p c ).

However the edge's vertices (encoded by mesh_p1_loc%jjs(:,:)) can belong to another processor (for example, in Fig 2 .4, for the bottom and left most cyan left, one of its edges is part of the external boundary but have one point attributed to the blue processor).

For an unrefined edge

m s ∈ [[1, M es (p c )]], the interface number ∈ [[1, N interf ace ]
] of the refined edges will be the same as m s . However, in order to get the cells and points of the refined edges, it is necessary to know the cell number of the unrefined edge. Using that the cell of m s is given by m =mesh_p1_loc%neighs(m s ), the cell number i of m s is obtained by finding which vertex of m does not belong to points of m s (given by mesh_p1_loc%jjs(:,m s )).

We will now consider that i = 1 so that we can use the unrefined edge e 1 in Fig. 2.6 as an example (meaning that we consider e 1 as being part of the external boundary of the global mesh). Using Fig. 2.6, it is now simple to find the index of the points and cells bordering the refined edges by just copying the values of the already built cells m 2 and m 3 . The local indices of m 2 and m 3 are simple to get by using that they are neighbours of the cell m.

End result

Once all these steps are done, we end up with the refined mesh called mesh_p1_loc_r, which is fully structured and respects the convention discussed in section 2.2.1. The unrefined mesh_p1_loc mesh is now deallocated. The refined mesh is renamed mesh_p1_loc. This refinement process can be repeated as many times as necessary to create finer and finer meshes (see Fig. 2.1).

. Creation of the distributed P 2 /P 3 mesh

The mesh_p1_loc mesh previously created (and refined if necessary) is still only defined on the P 1 finite elements. It is fine as is for the pressure field if we use P 1 /P 2 Hood-Taylor elements for the Navier-Stokes equations. However, we still need at least P 2 elements for the velocity field. Creating P 2 and (optionally) P 3 elements only means adding new points to all of the cells according to Fig. 2.3. This is very simple, as all that is required is to extend the following structures:

• mesh_p1_loc%jj(:,:) by adding the new points to the cells,

• mesh_p1_loc%loc_to_glob(:) by creating the processor and global index of the new vertices,

• mesh_p1_loc%rr(:,:) by creating the (r, z) coordinates of the new points,

• mesh_p1_loc%jjs(:,:) by adding the new points to the edges at the external boundaries of the global mesh,

• mesh_p1_loc%extra_jj(:,:) by adding their global index to the extra cells.

The newly created mesh will be called mesh_p2_loc and mesh_p3_loc in the following sections. In the code, the created meshes are named pp_mesh (for the pressure field) and vv_mesh (for the velocity field). For the various quantities defined on the new mesh, a superscript 1, 2, or 3 will be used depending on the finite elements considered. Again, we will be working with a processor

p c ∈ [[1, N p ]].
The number of points will now be:

• for P 2 , N 2 = N 1 + M 1 e ,
• for P 3 ,

N 3 = N 1 + 2 * M 1 e + M 1 .

For the bulk

For the rest of this section, we will be working with a cell with local index The local indices of the vertices do not change. However, we still need to create the other nodes and their indices. The cell numbering of the new nodes will follow the diagram shown in Fig. 2.3. For an edge e i l with cell number i ∈ {1, 2, 3}, the local index of the new nodes will be:

m l ∈ [[1, M (p c )]],
• for P 2 , v 2 n i l = e i l + N 1 (p c ),
• for P 3 , -for the one closest to the lowest cell number vertex (4,6,8 in Fig. 2.3),

v 3 n 1i l = e i l + N 1 (p c ),
for the one closest to the highest cell number vertex (5,7,9 in Fig. 2.3),

v 3 n 2i l = e i l + M 1 e (p c ) + N 1 (p c ),
for the point in the middle of the cell,

v 3 c l = m l +2 * M 1 e (p c )+N 1 (p c ).
and their position is created as follows, using x 1(i) (resp. x 2(i) ) the position of the vertex of lowest (resp. highest) cell number of e i l and x i the position of the vertex opposite of the edge e i l :

• for P 2 , mesh_p2_loc%rr(v 2

n i l )= (x 1(i) + x 2(i) )/2,
• for P 3 , -for the one closest to the lowest cell number vertex, mesh_p3_loc%rr(v 3

n 1i l )= (2 * x 1(i) + x 2(i) )/3,
for the one closest to the highest cell number vertex, mesh_p3_loc%rr(v 3

n 2i l )= (x 1(i) + 2 * x 2(i) )/3,
for the point in the middle of the cell, mesh_p3_loc%rr(v 3

n i l )= (x 1(i) + x 2(i) + x i )/3.
As with refinement, the global index of all points (old and new) is then created (or updated) by adding pc-1 k=1 N n (k) to their local index. This makes it possible to almost fully create mesh_p(2/3)_loc%jj, mesh_p(2/3)_loc%loc_to_glob and mesh_p(2/3)_loc%rr. Only the points at the border with a lower-rank processor now need to be created or updated.

At the border with a lower-rank processor

Creating the P 2 and P 3 meshes generates new nodes on edges that are not attributed to the processor p c but that are part of one of p c 's cells. Hence, the number of extra nodes becomes:

• for P 2 , N 2 extra (p c ) = N 1 extra + M 1 e extra ,
• for

P 3 , N 3 extra (p c ) = N 1 extra + 2 * M 1 e extra .
Firstly, the local index

v l ∈ [[N 1 (p c ) + 1, N 1 (p c ) + N 1 extra (p c )]
] of the node whose cell number is i ∈ {1, 2, 3} needs to be shifted by N As with refinement, the edge global numbering is e g =mesh_p1_loc%jees(e l ) and its attributed processor will be p(e g ). Hence the edge's local numbering in that processor p(e g ) is e l (p(e g )) = e g -p(eg)-1 k=1

M r e (k). The global indices of the new points are now easily created:

• for P 2 , v 2 ng = e l (p(e g )) + N 1 (p(e g )) + p(eg)-1 k=1 N 2 (k),
• for P 3 , -for the one closest to the lowest cell number vertex, v 3 n 1g = e l (p(e g ))+ N 1 (p(e g )) + p(eg)-1 k=1

N 3 (k),
for the one closest to the highest cell number vertex, v 3 n 2g = e l (p(e g ))+ M 1 e (p(e g )) + N 1 (p(e g )) +

p(eg)-1 k=1

N 3 (k).
The positions of these points are created in the same way as for the bulk by using the positions of the end points of this edge. This finishes the creation of mesh_p(2/3)_loc%jj, mesh_p(2/3)_loc%loc_to_glob and mesh_p(2/3)_loc%rr.

At the border with an higher-rank processor

The extra cells are also needed to use P 2 and P 3 finite elements. The same process as for the bulk is used but it is not necessary to calculate the position of the points. As with refinement, only the global index needs to be calculated.

We will now work with an extra cell with global index m g .

Similarly to the refinement process, every vertex already created (points (1,2,3) in Fig. 2.3) has its global index shifted by p(vg-1 k=1

N 2/3 (k)- p(vg)-1 k=1 N 1 (k)
where p(v g ) is its attributed processor.

For the points created on the edges, it is the exact same process as explained in the previous section. For an edge of global index e ig with i ∈ {1, 2, 3}, attributed processor p(e g ) and processor index in p(e g ) being e l (p(e g )), the global indices of the points created on this edge are:

• for P 2 , v 2 ng = e l (p(e g )) + N 1 (p(e g )) + p(eg)-1 k=1 N 2 (k),
• for P 3 , -for the one closest to the lowest cell number vertex, v 3 n 1l = e l (p(e g ))+ N 1 (p(e g )) + p(eg)-1 k=1

N 3 (k),
for the one closest to the highest cell number vertex, v 3 n 2l = e l (p(e g ))+ M 1 e (p(e g )) + N 1 (p(e g )) +

p(eg)-1 k=1

N 3 (k).

For the P 3 mesh, the global index of node 10 is given using the processor index of the cell: m l = m g -p(mg)-1 k=1 M r (k). The global index of this node will be :

• v 3 ng = m l + 2 * M 1 e (p(m g )) + N 1 (p(m g )) + p(mg)-1 k=1 N 3 (k).
It is thus possible to fully create mesh_p2_loc%extra_jj and mesh_p3_loc%extra_jj.

External boundaries

All that remains is to add the newly created nodes to the external boundaries of the total mesh and update the global index of the vertices. For every edge

m s ∈ [[1, M es (p c
)]] at the border, using that m l =mesh_p1_loc%neighs(m s ) is its cell local index, we just need to copy the global index of the relevant nodes given by mesh_p(2/3)_loc%jj(:,m l ) into mesh_p(2/3)_loc%jjs(:,m s ). The relevant nodes are found using the cell index i of the edge. As an example, for i = 2, using Fig. 2.3, the relevant nodes are (1,6,7,3) in that order.

. Benchmark on fine meshes for solving the Schur complement for the discrete Stokes problem

This section takes the form of an article currently submitted to Computer Methods in Applied Mechanics and Engineering (October 2023).

The aim is to implement and test several preconditioning methods for solving the Stokes problem. The end goal is to see if these methods can compete with the prediction-correction method currently used in the SFE-MaNS code in order to accelerate the simulation of the magnetohydrodynamics equations. This could potentially save a lot of computational time making it possible to see long-term evolution in the dynamo mechanisms. The conclusion is that preconditioning methods are still not competitive with predictioncorrection methods.

The meshes used for solving the Stokes problem in the following article are successive refinements of Fig. 2.1.

1. Introduction. The most popular time-stepping techniques for solving the time-dependent Navier-Stokes equations are the so-called projection methods initially proposed by Chorin Chorin [9] and Temam Temam [27]. Although significant progresses have been made since the original work of Chorin and Temam (see e.g., [19] for a brief overview of the literature on the topic), all these methods have however limited accuracy in time. Moreover, when combined with space discretization, they do not allow for a good control of the weak divergence of the velocity field. A better control on the accuracy and on the weak divergence of the velocity can be obtained by using more traditional time-stepping techniques combining implicit and explicit Runge-Kutta techniques. The downside of these techniques though is that one has to solve at each time step a generalized Stokes problem with a saddle point structure. Many iterative techniques to solve this problem have been proposed in the literature (see e.g., Elman et al. [11,Chap. 8] for a review of the literature). They are all more or less based on Uzawa iterations on the pressure Schur complement. The key problem with this type of methods is that the condition number of the Schur complement grows unboundedly as the time step and the viscosity go to zero. This in turn implies that performance of the method strongly depends on how it is preconditioned.

The objective of this paper is to revisit some preconditioning techniques for solving the generalized Stokes problem in the context of time stepping techniques. We revisit in particular the use of the augmented Lagrangian formulation introduced by Fortin and Glowinski [15] and we show that when combined with a preconditioning technique proposed by Cahouet and Chabard [8], one obtains an iterative method that has converging properties that are uniform with respect to the time step, the mesh size and the viscosity. The analysis is done with the assumption that the time step is proportional to the mesh size to reflect that in practical situations the nonlinearities are treated explicitly and therefore the standard CFL restriction on the time step has to be enforced to achieve stability. Thorough numerical tests are done using PETSC's parallel framework and mixed Hood-Taylor finite elements P 2 /P 1 and P 3 /P 2 . The tests are also done with the generic Newtonian form of the viscous tensor inducing nontrivial couplings between the Cartesian components of the velocity. This allow us to test the robustness of the method and its capability to handle nonlinear expressions for the stress tensor.

The solution of the generalized Stokes problem has a long history that started with the landmark papers of Fortin and Glowinski [15] and Cahouet and Chabard [8]. In [15] the authors introduce the augmented Lagrangian method and show that the conditioning of the pressure Schur complement improves as the scalar multiplier of the augmented Lagrangian term increases (cf. e.g., Golub and Greif [17, Prop. 2.1] Benzi and Olshanskii [3,Lem. 4.1] and Proposition 2.2 below). In [8], the authors propose a preconditioner of the Schur complement that behaves well in the two limits µτ h -2 → 0 and µτ h -2 → ∞ where µ is the viscosity, τ is the time step, and h is the meshsize. The above two approaches are combined in Benzi and Olshanskii [3], Benzi et al. [5], and Benzi and Wang [4]. In these papers the authors show that combining the Cahouet&Chabard preconditioner and the augmented Lagrangian technique gives a good preconditioning method of the Generalized Stokes problem. This method has been numerically evaluated in many papers like Farrell et al. [13], Moulin et al. [23], and Shih et al. [26], where it is indeed shown that the number of MINRES or GMRES iterations used to solve the pressure Schur complement problem behaves well with respect to the mesh size and the viscosity.

The key conclusion of the paper is that although the number of outer iterations to solve the Schur complement is uniform with respect to all the parameters (meshsize, time step, and viscosity), the time efficiency of the method is nowhere near that of projection methods. The main reason is that the method takes too many inner iterations to solve the velocity problem associated with the augmented Lagrangian formulation. The loss of efficiency cannot be reasonably offset by arguing that larger time steps can be used since the nonlinear advection can be made implicit in this context. Our conclusion is that, at the time of this writing, methods based on the augmented Lagrangian and grad-div stabilization methods cannot reasonably compete with projection methods or variations thereof on very large problems requiring iterative solutions of the velocity problem. Of course, some of these problems can be offset on small size problems by using direct solution methods.

The paper is organized as follows. First we formulate the problem in §2. We describe the discrete setting and propose a preconditionner combining ideas from Cahouet and Chabard [8] and Fortin and Glowinski [15]. We discuss in §3 key practical difficulties that arise when implementing this preconditioner. We particularly focus on the parallel aspects of the question. The full preconditioner and some of its variants are tested in §5. Its robustness with respect to the viscosity, the time step and the mesh size are illustrated.

2. Formulation of the problem. We formulate the problem and introduce the notation in this section.

2.1. Semi-discrete problem. We consider the generalized Stokes problem arising from the approximation in time of the time-dependent Navier-Stokes equations in a domain D ⊂ R d . Here D is an open bounded Lipschitz polyhedron. We denote by D the diameter of D. The exact form of the time discretization that is adopted to approximate the problem in time is irrelevant for our purpose. Let V be a closed subspace of H 1 (D) and Q be a closed subspace of L 2 (D). We assume that we are given at each time step a linear form f ∈ L(V ; R), and the problem is reduced to finding a pair (u, p) ∈ V ×Q so that

(2.1) a(u, p) -b(v, p) + b(u, q) = f (v), ∀(v, q) ∈ V ×Q,
where the bilinear form a and b are defined by

(2.2) a(v, w) := D 1 τ v•w + 2µe(v):e(w) dx, b(v, q) := D q∇•v dx.
Here e(v) := 1 2 (∇v + (∇v) T ) is the strain rate tensor, µ is the shear viscosity, and τ is the time step. The boundary conditions of the problem are encoded in the definition of the space V and the linear form f . This problem has a unique solution (see e.g., Girault and Raviart [16], Boffi et al. [6, § 4.2]). Our objective is to construct an iterative method for the solution of a discrete version of (2.1) that performs well with respect to the product τ µ and the meshsize.

Discrete problem.

In the rest of this paper we solely focus on a discrete version of (2.1). We consider two sequences of finite-dimensional vector spaces {V h } h∈H and {Q h } h∈H where the index set H is countable and has 0 as unique accumulation point. The velocity is approximated in V h and the pressure is approximated in Q h . To simplify some arguments we assume that the velocity approximation is conforming, i.e., V h ⊂ V . We also assume that the pair (V h , Q h ) is inf-sup stable, i.e., there is β h so that (2.3)

β h := inf q h ∈Q h \{0} sup v h ∈V h \{0} b(v h , q h ) v h H 1 (D) q h L 2 (D) > 0.
The discrete version of (2.1) consists of seeking u h ∈ V h and p h ∈ Q h so that the following holds true for all v h ∈ V h and all q h ∈ Q h :

(2.4) a(u h , v h ) -b(v h , p h ) = f (v h ), b(u h , q h ) = 0.
Let {φ i } i∈V and {ψ k } k∈Q be bases of V h and Q h (say, finite-element-based global shape functions). We define the velocity mass matrix M V , the velocity stiffness matrix E V , and the pressure mass matrix M Q by setting

M V,ij := D φ i •φ j dx, ∀i, j ∈ V (2.5a) M Q,ij := D ψ i •ψ j dx, ∀i, j ∈ V (2.5b) E V,kl := D 2e(φ k ):e(φ l ) dx ∀k, l ∈ Q. (2.5c)
To formulate the algebraic version of (2.4), we introduce the matrices associated with the bilinear forms a and b:

A ij := 1 τ M V,ij + µE V,ij = a(φ j , φ i ), ∀i, j ∈ V, (2.6) B kj := b(φ j , ψ k ), ∀k ∈ Q, ∀j ∈ V. (2.

7)

Denoting u h := i∈V U i φ i and p h := k∈Q P k ψ k the approximations of the velocity and the pressure, the discrete problem (2.4) is then equivalent to solving the following linear system:

(2.8)

AU -B T P = F, BU = 0.

The Schur complement.

As the coercivity of a implies that the matrix A is invertible, a standard way of solving (2.8) consists of constructing the Schur complement with respect to the pressure. Letting S := BA -1 B T , the problem (2.8) is equivalent to finding P and U so that (2.9)

SP = -BA -1 F, AU = F + B T P.
As this problem is in general solved iteratively, it is important to have some information on the 2 -condition number of S. The following result can be found in Elman et al. [11,Prop. 5.24] (see also [12,Prop. 50.14] for other informations on the topic).

Proposition 2.1 (Spectrum of S). Let σ(S) be the spectrum of S. Let µ min , µ max be the smallest and largest eigenvalues of the pressure mass matrix M Q . Let a and b be the norms of the bilinear forms a and b. Let α h be the coercivity constant of a. Then σ(S) ⊂ [µ min

β 2 h a , µ max b 2 α h ]. As a ∼ µ 2 D + 1 τ and α h ∼ µ 2 D
, the 2 -condition number of S scales like (1+

2 D τ µ ) b β h . This number is larger than (1 + 2 D
τ µ ) and grows unboundedly as

2 D
τ µ → ∞. Hence, if not properly preconditioned, iterative methods for solving SP = -BA -1 F have convergence rates that degrade as the viscosity decreases and the time step goes to zero.

The loss of convergence rate mentioned above can be mitigated by using the augmented Lagrangian technique introduced by Fortin and Glowinski [15]. This method replaces the original system (2.8) by the following equivalent system:

(2.10) (A + λµB T M -1 Q B)U -B T P = F, BU = 0,
where λ is a non-dimensional positive parameter. Here again, one can solve the problem by constructing the Schur complement of the pressure

S λµ := B(A + λµB T M -1 Q B) -1 B T , (2.11a) S λµ P = -B(A + λµB T M -1 Q B) -1 F, AU = F + B T P. (2.11b)
The main proporties of this method are summarized in the following standard result (cf. e.g., Golub and Greif [17, Prop. 2.1] Benzi and Olshanskii [3, Lem. 4.1]) Proposition 2.2 (Augmented Lagrangian). Let s and s be the smallest and largest eigenvalues of S. Then the following holds true:

S -1 ρ = ρM -1 Q + S -1 , (2.12a) σ(M -1 Q S ρ ) ⊂ [(ρ + s -1 ) -1 , (ρ + s -1 ) -1 ] (2.12b)
The above estimate of the spectrum of S λµ shows that the 2 -condition number of S λµ converges to unity as λµτ / 2 D → ∞. Of course the difficulty is to find a balance, because as the condition number of S converges to unity when λµτ / 2 D → ∞, the condition number of (A + λµB T M -1 Q B) becomes unrealistically large. In the rest of the paper we discuss preconditioning techniques for S λµ .

Preconditioning.

Ideas to precondition S λµ can be found by investigating how S λµ degenerates in the limit µτ /h 2 → ∞ and in the limit µτ /h 2 → 0.

In the limit µτ /h 2 → 0, we have

S λµ → τ B(M V ) -1 B T . The inverse of S λµ in this limit is then τ -1 (B(M V ) -1 B T ) -1 .
In the limit µτ /h 2 → ∞, we have

S λµ → µ -1 B(E V + λB T M -1 Q B) -1 B T . As E V and B T M -1 Q B have similar spectra, it seems natural to approximate S λµ in the limit µτ /h 2 → ∞ by µ -1 (1 + λ) -1 B(B T M -1 Q B) -1 B T . Hence we posit that the inverse of S λµ approximately behave like µ(1 + λ)M -1 Q when µτ /h 2 → ∞.
In conclusion, reasoning as in Cahouet and Chabard [8] (see also Bramble and Pasciak [7]), the preconditioner we consider for S λµ is

(2.13) S -1 λµ ∼ C λµ := τ -1 (B(M V ) -1 B T ) -1 + µ(1 + λ)M -1 Q .
In light of the identity (2.12a), the above expression is consistent with the approximation of S -1

0 by τ -1 (B(M V ) -1 B T ) -1 + µM -1
Q , which is the preconditionner of the Schur complement S 0 considered in Cahouet and Chabard [8, Eq. ( 44)]. In this paper, the authors also propose another preconditioner using the discrete pressure Laplacian instead of B(M V ) -1 B T , but we are going to show that it is important to use B(M V ) -1 B T to be uniform with respect of the three parameters τ , µ and h (see preconditioner CaCh 2 defined in (5.2b) and left panel in Figure 2).

3. Solving the velocity problem. We discuss in this section implementation and performance issues associated with the solution to the velocity problem which consists of solving (A + λµB T M -1 Q B)X = F. We use Dirichlet boundary conditions over the entire boundary of the computational domain; that is, V := H 1 0 (D). 3.1. Numerical details. We focus on approximation settings with a continuous representation of the pressure. All the test reported in this section and in §4, §5 are done with continuous Lagrange finite elements in two space dimensions on nonuniform triangular meshes. The parallel implementation is done with PETSC [2]. All the errors that are reported are relative.

Depending on the situation we use either the direct linear solver MUMPS (see Amestoy et al. [1]) or the iterative algebraic multigrid solver BoomerAMG (see Henson and Yang [21]). For completeness and reproducibility we give in Table 1 the options of BoomerAMG that we use.

-pc type hypre -pc hypre type boomeramg -pc hypre boomeramg strong threshold 0.7 or 0.1 -pc hypre boomeramg coarsen type Falgout -pc hypre boomeramg relax type all Chebyshev Table 1: BoomerAMG options used in the paper.

In all the tests reported in the paper, the velocity and pressure mass problems involving the matrices M V and M Q are systematically solved by using BoomerAMG with the strong threshold set to 0.1.

Preconditioners for the velocity problem. It is very difficult to assemble the matrix coefficients of A+λµB

T M -1
Q B when the approximation of the pressure is based on continuous finite elements (like for the P k+1 /P k Hood Taylor elements, k ≥ 1). The first reason is that the inverse of the pressure matrix M -1

Q cannot be computed. It can be approximated by using some lumping technique though. The second reason is that the sparsity pattern of the matrix is doubled. This drastically increases the memory footprint as well as the communications needed for parallel computing. Therefore, it is significantly easier to solve iteratively this problem by invoking only matrix vector multiplications not requiring any assembling. We now test this idea.

As we use Dirichlet boundary conditions, the bilinear form D e(u):e(v) dx is also

equal to D ∇u:∇v dx + D ∇•u∇•v dx because two integrations by parts show that D (∇u) T :∇(v) dx = D ∇•u∇•v dx for all u, v ∈ H 1 0 (D).
Hence defining L V to be the stiffness matrix associated with the bilinear form D ∇v h :∇w h dx and D to be the stiffness matrix associated with the grad-div bilinear form

D ∇•v h ∇•w h dx; i.e., L V ,ij := D ∇φ i :∇φ j dx, ∀i, j ∈ V, (3.1)
D ij := D ∇•φ i ∇•φ j dx, ∀i, j ∈ V, (3.2)
we obtain the following equivalent representation of the matrix A:

(3.3) A = τ -1 M V + µ(L V + D),
where M V is the velocity mass matrix. This leads use to consider the following three candidates to precondition the matrix

(A + λµB T M -1 Q B): τ -1 M V + µL V + µ(1 + λ)D = A + λµD, (3.4) τ -1 M V + µL V + µD = A, (3.5) τ -1 M V + µL V . (3.6)
In the first preconditioner, we replace B T M -1 Q B by D. In the second preconditioner, we remove

B T M -1 Q B.
In the third preconditioner, we remove both D and B T M -1 Q B. The stencil for these three preconditioners is standard when using Lagrange finite elements. The matrix coefficients can be easily assembled by looping over the mesh cells. These coefficients can in turn be used to construct incomplete LU factorizations or any other approximations.

Choice of velocity preconditioner.

We now test the three preconditioners (3.4)-(3.6) with continuous Lagrange finite elements in two space dimensions. Setting k := 4π, we consider the divergence-free vector field u and source term f

(3.7) u(x) = sin(kx) sin(ky) cos(kx) cos(ky) , f (x) = u τ -2µe(u).
Recalling that the velocity shape functions are denoted {φ i } i∈V , we construct the vector F with entries

F i = D f (x, t)•φ i (x)
dx and solve the linear system

(3.8) (A + λµB T M -1 Q B)U = F.
Setting u h := i∈V U i φ i , the relative errors in the L 2 -norm reported below are defined to be equal to u -

u h L 2 (D) / u L 2 (D) .
We perform tests with P 2 Lagrange elements. The tests are done on five nonuniform meshes composed of 118837, 474065, 1893697, 7569665, 30268417 velocity grid points, respectively. In order to test the weak scalability of the preconditioning, the computations are done on those fives meshes using 2, 8, 32, 128, 512 processors, respectively. We use the parameters µ ∈ {1, 10 -2 , 10 -4 }, with λ = 1. We also use τ = 1 as our goal is just to test the preconditioners defined above.

The problem (3.8) is solved by using the GMRES solver in PETSC with threshold 10 -10 . The linear systems induced by the action of the preconditioners are solved by using the iterative algebraic multigrid solver BoomerAMG. The options of Boomer-AMG are given in Table 1. We take the strong threshold in BoomerAMG to be equal to 0.7 for the first and second preconditioners (see (3.4)-(3.5)) and equal to 0.1 for the last one (see (

of A + λµB T M -1 Q B by τ -1 M V + µL V , (with τ = 1, λ = 1
, and P 2 Lagrange elements for the velocity), using BoomerAMG, strong threshold 0.1.

The results obtained with the preconditioner (3.4) are shown in Table 2. The results with the preconditioner (3.5) are shown in Table 3. The results with the preconditioner (3.6) are shown in Table 4. We show in these tables the number of GMRES iterations that are necessary to reach the assigned (relative) threshold, which we recall is 10 -10 . We also show the L 2 -norm of the relative error on u and the wall-clock time in seconds to perform each computation.

In all the cases, we observe that the number of GMRES iterations are small. This number clearly converges to 1 for the third preconditioner as the number of grid point goes to infinity. This is because the tests are done with a manufactured velocity field that is divergence-free. We also observe that the L 2 -norm of the relative error scales like O(h 3 ) which is the expected rate of convergence for P 2 approximation. Notice though that for the firt preconditioner, the error has a strange behavior on the last two meshes. We have verified that to get the proper behavior one has to lower the GMRES threshold to 10 -14 . This behavior is not observed for the other two preconditioners.

Notice also that, although the number of GMRES iterations for all the preconditioners are roughly the same and very small, the wall-clock time for the first and second preconditioners do not scale well. This is due to the difficulty to coarse-grain the matrix D in the BoomerAMG algorithm. Entirely getting rid of the contribution of the grad-div operator as done in the third preconditioner gives a methods that scales perfectly well (see Table 4). Note that the accuracy and wall-clock time for the third preconditioner are robust with respect to µ, which is not the case of the other two. This series of tests shows in passing that just looking at the number of GMRES iterations, as often done in the literature, is not fully informative. The wall-clock time is also an important factor to differentiate preconditioners.

As the third method works optimally, we did an additional test on a mesh composed of approximately 120 millions grid points using 700 processors. If we had 2048 processor available, the wall-clock times would have been 700

4×512 × 12.4s = 4.24s for µ = 1, 700 4×512 × 11.3s = 3.86s for µ = 10 -2 , and 700 4×512 × 7.2s = 2.46s for µ = 10 -4 . We therefore conclude that the scaling is almost perfect. This is also well illustrated in Figure 1 (see the horizontal green curves). Note that for the last two meshes, the L 2 -norm of the relative error on u is of the same order as the GMRES threshold. To get back the proper O(h 3 )-scaling of the error one would have to decrease the GMRES threshold to 10 -14 .

The main conclusion of the series of tests done in this section is that although the grad-div operator has very interesting properties emphasized in the literature (Benzi and Olshanskii [3], Olshanskii et al. [24], Guermond and Minev [18], de Frutos et al. [10]), solving problems requiring the solution of linear systems involving the matrix D or variations thereof is difficult (see e.g., Heister and Rapin [20], Jenkins et al. [22], Fiordilino et al. [14]). This conclusion aligns with the recommendations made in Benzi et al. [5] that dropping the coupling blocks in D greatly improves the efficiency of the velocity preconditioner.

BoomerAMG vs. MUMPS.

In this section we compare the use of iterative and direct solvers for the solution of the linear systems induced by the velocity preconditioners. We only test the third preconditioner, τ -1 M V + µL V , as we have shown in the previous section that it is the most efficient one. We compare the direct solution method MUMPS with BoomerAMG. We use the same parameters as in §3.3; that is, λ = 1 and τ = 1.

The results for MUMPS are shown in Table 5. We recall that those for Boomer-AMG are in Table 4. We observe that the number of GMRES iterations and the L 2 -norm of the relative error are the same as when using BoomerAMG. A key difference with BoomerAMG is that MUMPS performs first the LU factorization and then does the backward substitution. None of these two operations scale very well. We

of A + λµB T M -1 Q B by τ -1 M V + µL V , (with τ = 1, λ = 1
, and P 2 Lagrange elements for the velocity), using MUMPS.

only report here the backward substitution time. We observe that MUMPS' backward substitution time is faster than BoomerAMG for meshes composed of less than 2 million points. But, the memory footprint of MUMPS grows very fast with the size of the problem, making it impossible for us to solve the problem on the mesh composed of only 30 million velocity grid points. A+λµD, BoomerAMG A, BoomerAMG ); µ = 10 -2 ( ), µ = 10 -4 ( ).

τ -1 M V +µL V , BoomerAMG τ -1 M V +µL V , MUMPS
In Figure 1, we compare the efficiency of MUMPS using the third preconditioner with the efficiency of BoomerAMG using the three preconditioners discussed in §3.3. The X-axis shows the number of velocity grid points. The Y -axis of the figure shows the elapsed time multiplied by the number of processors and divided by the total number of degrees of freedom (two times the number of velocity grid points plus the number of pressure grid points). The result is shown in millisecond. In the ideal case, this number should be independent of the number of processor and the number of degrees of freedom. We observe that the only method for which this is indeed the case is the third preconditioner solved with BoomerAMG (see the horizontal green curves).

Finally, as already mentioned in §3.2, Figure 1 shows that the wall-clock time per degree of freedom and time step of the other two preconditioners does not behave optimally, whereas it does for the third preconditioner.

3.5. Tests with non divergence-free solution. We now confirm the conjecture made in §3.3 that the number of GMRES iterations to solve the problem (3.8) converges to 1 as the number of grid point goes to infinity because the solution defined in (3.7) is divergence free.

In this section we test the third preconditioner with a right-hand side F constructed with a velocity field that is not divergence-free. More precisely, we take

(3.9) u(x) = sin(2kx) sin(ky) cos(kx) cos(ky) , f (x) = u τ -2µe(u) -λµ∇∇•u,
with k = 4π, λ = 1, and τ = 1. Here again we solve (A + λµB T M -1 Q B)U = F and the entries of the vector F are given by

F i = D f (x, t)•φ i (x) dx.
The tests are done with P 2 and P 3 finite elements for the velocity. For the P 3 elements the tests are done with meshes composed of 266902, 1065685, 4258897, 17027905, and 68096257 velocity grid points, respectively. Table 6: Non divergence-free solution. Preconditioning of A + λµB T M -1 Q B by τ -1 M V + µL V , (with τ = 1, λ = 1) using BoomerAMG, strong threshold 0.1. Top: P 2 Lagrange elements for the velocity. Bottom: P 3 Lagrange elements for the velocity The results are shown in Table 6. At variance with what is observed in Tables 234, the number of GMRES iterations is bounded from below away from 1. Depending on the value of µ, the number of GMRES iterations needed to reach the 10 -10 threshold varies between 12 to 16. This number is remarkably constant when µ is fixed. The L 2 -norm of the relative error properly scales like O(h 3 ) for the P 2 approximation and like O(h 4 ) for the P 3 approximation until the 10 -10 threshold is reached. Finally, we observe that the weak scaling is satisfactory for both the P 2 and P 3 approximations.

That the velocity preconditioner takes more time to solve a problem where the solution has a nonzero divergence will impact the performance of the preconditioner of the pressure Schur complement.

3.6. Solving the velocity problem with discontinuous pressure. In the case of an approximation of the pressure based on discontinuous finite elements (like the Scott-Vogelius elements [25] and their generalizations) the coefficient of the matrix B T M -1 Q B can be computed. This is possible because the support of the discontinuous finite elements for the pressure are localized on each mesh cell; as a result, the pressure matrix M Q can be locally inverted on each cell. Hence, it is possible to construct the coefficients of the matrix associated with the velocity problem. Note in passing that, contrary to B T M -1 Q B, the stencil of the pressure matrix BM -1 V B T is not simple, which means that computing the entries of this matrix in parallel is non trivial (see discussion in §4.3).

4. Implementation of the preconditioner C λµ . We discuss in this section the implementation of the preconditionner (2.13). We mainly focus our attention on the matrix BM -1 V B T as the linear algebra associated with the pressure mass matrix M Q is simple.

4.1. The matrix BM -1 V B T .
Here again, the matrix coefficients of BM -1 V B T cannot be easily computed when using continuous elements for the pressure. Even if the velocity mass matrix is lumped, the stencil of the resulting simplified matrix is double. Hence, the only practical way to solve the linear system BM -1 V B T X = Y is to use a preconditioned iterative method. A natural preconditioner for this matrix is the pressure Laplacian L Q whose entries are given by (4.1)

L Q,kl := D ∇ψ k •∇ψ l dx ∀k, l, ∈ Q.
However, this matrix is singular when Dirichlet boundary conditions are enforced on the velocity over the entire boundary of the domain D. As we only want to construct a preconditioner, one can consider instead the matrix with entries Table 7: Tests for the preconditioning of BM -1 V B T by M Q + L Q for P 1 and P 2 Lagrange elements, using BoomerAMG, strong threshold 0.1.

M Q + L Q with > 0.
As in §3.2, we illustrate this idea on continuous Lagrange finite elements in two space dimensions. We consider P 2 /P 1 and P 3 /P 2 elements. Recall that these finite elements are legitimate for the Navier-Stokes equations. For each mesh we build the right-hand Y := BM -1 V B T X by applying the matrix BM -1 V B T to some given pressure field. Here X is the vector whose entries are the value of the following pressure field at the pressure grid points: (4.2) p(x, y) := sin(4π(xy)).

We use = 1 in the preconditioner M Q + L Q . We use BoomerAMG with strong threshold equal to 0.1 to solve the linear system induced by the preconditioner. We show in Table 7 the L 1 -norm of the relative error when the GMRES threshold is reached. As the GMRES threshold is 10 -10 , we observe that the errors are of order 10 -10 , (these are not approximation errors since there is noting to approximate). The tests are done on five different meshes with mixed Lagrange elements. This test shows that M Q + L Q is an excellent preconditioner of BM -1 V B T as the number of GMRES iterations is small and decreases as the mesh is refined.

Using a lumped mass matrix instead M -1

V . Even though Table 7 shows that the number of GMRES iterations for solving BM -1 V B T P = F is small, the time for solving this problem may seem a bit too high for such a simple problem. The analysis of the problem reveals that most of the computation time is spent inverting the velocity mass matrix M V . As the matrix BM -1 V B T is only used as part of the preconditioner for S λµ , replacing the mass matrix M V by a lumped one could save some computation time. To explore this idea, we then propose to replace M V by a diagonal matrix Λ V with entries defined by

(4.3) Λ V ,ij := δ ij D |ϕ i (x)| dx, ∀i, j ∈ V.
Taking the absolute value of the shape function is important here as P 2 and P 3 velocity shape functions are not uniformly positive. Then instead of solving BM -1 V B T P = F, we propose to solve BΛ -1 V B T P = F while again preconditioning the system with Table 8: Tests for the preconditioning of BΛ -1 V B T by M Q + L Q for P 1 and P 2 Lagrange elements, using BoomerAMG, strong threshold 0.1.

M Q + L Q . P1 Nb.
We show the results of this test in Table 8. We observe that even though using Λ V instead of M V takes more GMRES iterations, the wall-clock time is much smaller and the method scales very well. Once again the tests are purely algebraic; as a result, the errors are proportional to the GMRES threshold.

Discontinuous finite elements for pressure. The matrix coefficients of

BM -1 V B T are given by i,j∈V D ψ k ∇•φ i dx(M -1 V ) ij D ψ l ∇•φ j dx, for all k, l ∈ Q.
If the velocity mass matrix is lumped or replaced by a diagonal matrix with the right scaling, the computation of the matrix's coefficients are possible. However, we observe that the velocity degrees of freedom supported at the mesh interfaces interact with all the pressure degrees of freedom on the cell touching this interface. As a result, for each cell K, the stencil of the pressure degrees of freedom located on K are composed of all the pressure degrees of freedom located on all the cells touching K (either by a face, an edge or a vertex). Hence the pressure stencil is significantly larger than that of the velocity. This type of stencil is very inconvenient for parallel processing. In conclusion, it is again more practical to use an iterative method preconditioned by a discontinuous Galerkin approximation of the Laplacian.

5. Numerical illustration of the performance of C λµ . We now illustrate the performance of the method (2.11b) using the preconditioner C λµ defined in (2.13). We compare various realizations of the method. We demonstrate numerically that (2.13) is indeed a preconditioner that is robust with respect to µ, τ and the meshsize. To be representative of actual situations where the time step deceases like the mesh size due to nonlinearities being made explicit in time, we set

(5.1) τ = N -1 2 ,
where N is the total number of grid points for the velocity. Here N -1 2 scales like the mesh size because the tests are done in dimension two.

5.1. Tests without augmented Lagrangian. We start by illustrating the method proposed in Cahouet and Chabard [8, Eq. ( 44)]. This method does not use the augmented Lagrangian, which is equivalent to take λ = 0 in (2.11b) and (2.13). The authors consider two preconditioners for S λµ which we call CaCh 1 and CaCh 2 :

CaCh 1 := τ -1 (B M -1 V B T ) -1 + µM -1 Q , (5.2a) CaCh 2 := τ -1 (L Q ) -1 + µM -1 Q . (5.2b)
It is observed in [8] that CaCh 1 outperforms CaCh 2 when µ is large. The test reported therein also suggest that these two preconditioners are somewhat equivalent when the viscosity µ is small, (see Figure 6 in [8]). We show in this section that this is not the case on fine meshes. We report in Figure 2 the GMRES residual (reported by PETSC) as a function of the iteration count for various small values of the viscosity (µ ∈ {1, 10 -2 , 10 -4 }) and various mesh refinements. The simulations are done with mixed P 2 /P 1 continuous finite elements. We test four meshes with the following velocity/pressure grid point counts: Mesh 1 (118837/29870); Mesh 2 (474065/118837); Mesh 3 (1893697/474065); Mesh 4 (7569665/1893697). As the velocity field is two dimensional, the total number of degrees of freedom is for each mesh: Mesh 1 (267544); Mesh 2 (1066967); Mesh 3 (4261459); Mesh 4 (17033027). We show in the left panel of Figure 2 the performance of the method with the preconditioner τ

-1 (L Q ) -1 + µM -1
Q . The performance of the method with the preconditioner CaCh 1 is shown in the right panel of Figure 2. We observe that the preconditioner CaCh 2 does not behave properly. The convergence rate strongly depends on the viscosity and the mesh size. The method seems to converge well on coarse meshes, but its performance becomes erratic when the meshes are refined. The method behaves well over the ten to twenty first iterations and then the convergence rate flattens. On the other hand, the method based on the preconditionner CaCh 1 works reasonably well. Notice though that the performance is not uniform with respect to µ and there seems to be a slight dependency on the mesh size.

We report in Table 9 the performance of the method using the preconditioner CaCh 1 with the stopping threshold equal to 10 -10 . We show in this table the number of GMRES iterations that are required to reach the assigned threshold. We also show the relative L 1 -norm of the error on the pressure, the relative L 2 -norm of the error on the velocity, and the L 1 -norm of the weak divergence. We observe the expected second-order convergence rate on the pressure and third-order convergence rate on the velocity. We also observe that the weak divergence is well below the assigned threshold. Table 9: Preconditioner CaCh 1 with λ = 0.

Tests with augmented

Lagrangian. We now test the method with the augmented Lagrangian term, i.e., we solve (2.11b) with the preconditionner (2.13) where M V is replaced by Λ V , i.e., (5.3)

C λµ := τ -1 (BΛ -1 V B T ) -1 + µ(1 + λ)M -1 Q .
We test the method with P 2 /P 1 and P 3 /P 2 elements with λ ∈ {1, 10} and compare the results with those obtained with the preconditioner CaCh 1 (which we recall corresponds to setting λ = 0). For the P 2 /P 1 approximation we use the meshes already considered above with the following velocity/pressure grid point counts: Mesh 1 (118837/29870); Mesh 2 (474065/118837); Mesh 3 (1893697/474065); Mesh 4 (7569665/1893697). For the P 3 /P 2 approximation we use meshes with the following velocity/pressure grid point counts: Mesh 1 (66937/29857); Mesh 2 (266785/118785); Mesh 3 (1065217/473857); Mesh 4 (4257025/1892865).

We report the results in Figure 3. Here again we show GMRES residuals versus the GMRES iteration count. We observe now that the convergence rate improves as λ increases and is almost uniform with respect to the mesh size and the viscosity for λ = 10. It takes between 7 to 15 iteration to reduce the residual by the factor 10 -10 when λ = 1, whereas it only takes between 7 to 9 iteration for λ = 10. Hence, the method seems to be close to be optimal when one only look at the GMRES iteration count.

λ = 0 λ = 1 λ = 10
Fig. 3: GMRES residual vs. iteration count for the preconditioner

τ -1 (B M -1 V B T ) -1 + µ(1 + λ)M -1
Q with µ = 1 (blue), µ = 10 -2 (red), and µ = 10 -4 (green). Mesh 1 ( ); Mesh 2 ( ), Mesh 3 ( ); Mesh 4 ( ). Top: Hood-Taylor P 2 /P 1 finite elements. Bottom: Hood-Taylor P 3 /P 2 finite elements. Leftmost column: Cahouet&Chabard preconditioner (λ = 0). Center column: Augmented Lagrangian with λ = 1. Rightmost column: Augmented Lagrangian with λ = 10.

Finally, we compare the actual efficiency of the Cahouet&Chabard method CaCh 1 (blue lines) with that of the preconditioned augmented Lagrangian methods with λ = 1 (red lines) and λ = 10 (green lines) . For these three methods, we show in Figure 4 the wall-clock time multiplied by the number of processors divided by the total number of degrees of freedom in millisecond: (5.4) wall-clock time×nb. processor nb. vel.+press. degrees of freedom .

We note also that all the methods scale well. But, quite surprisingly, we observe that all the methods have roughly the same efficiency; that is, between 0.4 to 2 millisecond is spent per degree of freedom on average to solve the problem, irrespective of the size of the problem and the number of processors. Even though, the number of GMRES iteration count is larger for the Cahouet&Chabard method, it takes much more time to solve the velocity problem in the augmented Lagrangian methods as discussed in section 3. Hence, the Cahouet&Chabard method seems to be the most efficient on average. The efficiency being between 0.4 to 2 millisecond can be compared to the implementations of the preconditioning method proposed by Benzi et al. [5]. In the first implementation reported in Benzi and Wang [4], tests done on a two-dimensional uniform Cartesian mesh (256×256 grid with Q 1 /Q 1 finite elements) and reported in Table 3.12 therein, give an efficiency ranging between 0.2ms and 0.8ms per degree of freedom for the stopping threshold 10 -6 (which would give 0.33 to 1.33 millisecond with the stopping threshold 10 -10 ). Similarly using the data from Table 5.5 in Farrell et al. [13], we obtain 0.2ms milliseconds for a two-dimensional mesh composed of 1.6×10 5 degrees of freedom with the stopping threshold of 10 -6 (which would give M. Creff, J.-L. Guermond 0.33ms with the stopping threshold 10 -10 ). For a two-dimensional test reported Figure 7 in Moulin et al. [23] using 14828 triangles with P 2 /P 1 elements (roughly 37000 degrees of freedom) and threshold 10 -3 , we get around 4ms per degree of freedoms (which would give 13ms with the threshold 10 -10 ). Using the three-dimensional test reported in Figure 9 therein with 75×10 6 degrees of freedom and threshold 10 -6 , over 2048 processors, we obtain 20ms. Fig. 4: . Y-axis:

wall-clock time×nb. processor nb. vel.+press. degrees of freedom in millisecond; X-axis: number of velocity grid points. µ = 1 ( ); µ = 10 -2 ( ), µ = 10 -4 ( ).

Conclusion.

We confirmed in the paper that the preconditioned augmented Lagrangian method requires less GMRES iterations than the traditional Cahouet&Chabard method. However, too much time is spent solving the velocity problem in the augmented Lagrangian method. The main difficulty consists of inverting the discrete version of the grad-div operator. Overall, the wall-clock efficiency of the augmented Lagrangian method is of the same order as that of the Cahouet&Chabard method.

For the CahouetCahouet&Chabard preconditioning, we showed that, using the discrete pressure Laplacian instead of the full matrix B M -1 V B T gives a methods that converges slowly on fine meshes and when the thresholds is stringent. Hence, we recommend using B M -1 V B T in the preconditioner even if this means solving another linear problem using another iterative method.

From a numerical perspective, properly using an algebraic multigrid solver like BoomerAMG yields very good parallel scaling and makes it possible to solve the problem on very fine meshes with no memory problem (for instance, 120 millions velocity grip points with 700 processors). In contrast, direct linear solvers are more performant on coarser meshes (under 2 millions velocity grid points).

Finally, the results about the time efficiency reported in this paper show that solving the time-dependent Stokes and Navier-Stokes equations by means of a pressure Schur complement technique does not yet realistically compete with pressureand velocity-correction methods, as the time efficiency for projection methods is usually of the order of 10 -2 ms per degree of freedom per time step with the current hardware technology. It might be possible to fine tune the preconditioned augmented Lagrangian method to reach 0.1ms per degree of freedom per time step, but gaining two orders of magnitude seems unfortunately problematic at the moment without a genuine algorithmic breakthrough. Of course one can argue that Schur complement techniques are more accurate in time than projection methods as they do not induce any time splitting error. Moreover, they can be used with larger time steps as they can be implemented with unconditionally stable time stepping technique.

-Post-processing tools

The SFEManS code by itself only gives access to the velocity field u, the pressure field p, and the magnetic field H. Hence, I had to implement new tools to study the energy transfers between the velocity and magnetic fields. First, the various quantities we are interested in (see 4) need to be calculated in SFEMaNS. Afterward, they have to be transformed into a usable file format for either visualisation or statistical study. Finally, as we are dealing with a great amount of data, all visualisation and analysis have to be automatized. This automation was done using Python and was coupled with Paraview for visualisation.

In order to easily transfer and use those tools, they are available for our research group in our private GitLab as different projects. A README.md file is included in each of them with a simple tutorial. The goal of this chapter is not to make a complete documentation of these tools but to present key points of each step of the workflow.

. Post-processing calculation in SFEMaNS

While the time simulation in SFEMANS is done, only the magnetic, velocity and pressure fields are stored. However, further computation is needed to get the different energy transfers happening in our system. In order to keep a high accuracy and efficiency for the transfer terms, computations of these new terms are still done using the SFEMaNS framework. The corresponding project is called 'DR_MHD' in GitLab. For ease of use, all new options needed are written in a 'data_example' file that just needs to be copy-pasted into the 'data' file used for the temporal simulation. There is no need to modify the various options previously used for the temporal simulation.

. Filtering process

In order to analyze what happens at different scales, we need to differentiate between large and small structures. This is important in order to find rare events in the flow that could indicate possible singularities. It also, make it possible to study how the energy transfers from one structure to another and hence find where and when important events for dynamo and turbulence happens. To do this separation, we use a filter. It is the most important part of the computation and the most costly in resources in terms of memory and computation cost.

The SFEMaNS grid is unstructured and is build to solve linear differential equations easily. Hence we define the filtered field g = u ℓ as the solution of:

g -ℓ 2 ∆g = u, (3.1)
with u a vector or scalar field and ℓ ∈ R + .

The filtering process is then transformed into a linear problem using the same method as for the evolution of the temperature (see section 2.1.2). Equation (3.1) becomes for the Fourier mode m of a cosine component of a scalar or vector field:

A m,c ij g m,c j = u m,c i , (3.2)
with:

A m,c ij = r,z r(ψ j ψ i + ℓ 2 (∂ r ψ j ∂ r ψ i + ∂ z ψ j ∂ z ψ i + m 2 r 2 ψ j ψ i )). (3.3)
As it takes a lot of time to build this matrix to solve the system, one matrix is saved in memory for each scale ℓ. However, it is a simple way to construct a filter that can be computed on unstructured meshes efficiently by using the powerful tools of linear algebra. This filter was already programmed by Hugues Faller [START_REF] Guermond | SFEMaNS documentation[END_REF] but a few modifications have since been done. Firstly, the same filter has been duplicated for the magnetic mesh as it only existed for the velocity mesh. This at least doubles the usage of memory for saving the matrices. Secondly, the boundary conditions have been changed from Dirichlet's to Neumann's in order to correct a problem that had occurred for ℓ = 0. Thirdly, the filter has been encapsulated in its own function name 'filter' for ease of use.

Currently, in order to choose at which scale to perform the calculation, a list of floats is asked from the user. Those floats are then multiplied by h max , the size of the biggest cell in the mesh. This makes it difficult to specify filters at characteristic lengths of the flow (Kolmogorov's scale for example). Modifications are underway to make it more practical.

. Computation of the various energy transfers

The biggest part of the implementation was the computation of each individual term of energy transfer. As it is only a succession of scalar and vector products and computations of a few derivatives, most of the framework was already in place. One of the difficulties was that some fields are defined on the magnetic mesh while others are on the velocity mesh. Hence, a few functions needed to be revamped in order to work with both meshes. As calculating all the energy transfers is costly, particular attention has been put toward memory usage in order to save computation time and memory space.

In order to minimize the cost of operations (filtering, deriving fields, and so on), simple optimization strategies are used, other than parallelization. When calculating different transfers, the quantities that are independent of the scale ℓ are kept in memory. This is important when calculating the energy transfers for different scales as it reduces the number of operations done. As an example: ∇ × 1 σ , ∇ × b, j × b, and others are stored when doing the calculation for the first scale ℓ so that they can be used for all other scales. However, storing quantities puts a strain on computer memory and can reduce performance. As such, when possible, generic tables for vectors and scalars are used. The aim is to minimize the number of these storages by cleverly designing the order of all operations. To make it practical to follow, a file has been created (named 'memory_usage.ods') that tracks which table contains which quantity at each operation step. Figure 3.1 shows the way this usage was planned for the scalar tables on the velocity mesh. To illustrate the practical implementation of the calculation, here is the order of operations for the calculation of two main terms:

T u→b = 1 4 (u ℓ • (b × j) + u • (b × j ℓ ) + (u × b ℓ ) • j + (u × b) • j ℓ ), (3.4)
D b = 1 4 (u ℓ • (b × j) + u • (b × j ℓ ) -(u × b ℓ ) • j -(u × b) • j ℓ ). (3.5)
The fields recorded from a previous simulation are u the velocity field, b the magnetic flux (or induction), and H the magnetic field on the nodes.

• Operations only done once:

-Memory allocation for all used structures.

-Create ghost vectors and scalar structures for PETSC parallelization.

-Initialize MPI parallelization for the Fourier modes.

-Assemble matrices for transfer between nodes and Gauss points on both the magnetic meshgrid and the velocity one.

-Calculation of the radius of the Gauss points (r(l)) on both meshgrids.

• For the first scale ℓ:

-Send b from the magnetic meshgrid to the velocity meshgrid and save it.

-Send H on the Gauss points and calculate j = ∇ × H on these points.

-Transfer j from the Gauss points to the nodes and save it.

-Calculate j × b and u × b using fast Fourier transform and save it.

-Send j × b from the magnetic meshgrid to the velocity meshgrid and save it.

• For all scales ℓ:

-Prepare the matrix to solve the filtering problem for the scale ℓ and save it.

-Calculate u ℓ in a generic vector n°1 (it will be used to compute other quantities before being used in T u→b ).

-Calculate j ℓ on the magnetic meshgrid and send it on the velocity meshgrid in a generic vector n°2.

-Calculate u ℓ • (b × j) and (u × b) • j ℓ and store them in generic scalars n°1 and n°2.

-Calculate b × j ℓ and store it in the generic vector n°1 as u ℓ will not be used anymore.

-Calculate u • (b × j ℓ ) and store them a generic scalar n°3.

-Sum generic scalars n°1 and n°3 into n°1 to get

u ℓ • (b × j) + u • (b × j ℓ ).
-Calculate u × b ℓ and store it in the generic vector n°1.

-Calculate (u × b ℓ ) • j and store them a generic scalar n°3.

-Sum generic scalars n°2 and n°3 into n°2 to get

(u × b ℓ ) • j + (u × b) • j ℓ .
-Sum generic scalars n°1 and n°2 and multiply by 1 4 into generic scalar n°3 to get T u→b .

-Subtract generic scalars n°1 and n°2 and multiplying by 1 4 into generic scalar n°3 to get D b .

For simplicity's sake, it might not match the numbers in 'memory_usage.ods' (Fig. 3.1) as the savings on the magnetic mesh are not discussed here.

. Output files

The computed fields are saved in the SFEMaNS format which makes it easily readable by SFEMaNS. The naming convention is as follows:

suite_${field_name}_nDR${filter_number}_S${domain_number}_I${iteration_number}.$${mesh_name}.FEM
where 'field_name' is the name of the field (for example 'Ec' for the kinetic energy), 'filter_number' is associated to the scale's number for the filter, 'do-main_number' is the number associated to the processor in the meridian plane, 'iteration_number' is the temporal iteration and 'mesh_name' is the name of the SFEMaNS mesh used for the temporal simulation. Fields that do not have filtering (as for example the helicity H) have the filtering number set to 1.

The last output file is named 'fort.75' and contains all the spatial averages of all fields for each iteration number. It also has additional information such as the simulation time at each iteration number and the scale of filtering. A description of 'fort.75' is written in 'means info' (Fig. 3.2) which tells which column corresponds to which quantity.

. Problem with chaining products

Simple tests have been put in place in order to verify that all fields have been correctly computed. For example, for ℓ = 0, a few terms are expected to be equal to zero. This works well for all except for D b (Fig. 3.3). The cause is the product of fields. Indeed, as the fields are decomposed in Fourier space in the azimuthal direction, it is necessary to project back into the physical space for computing products. However, if not enough Fourier modes are used when projecting, aliasing happens. Indeed, the product of two vector having N ∈ N + Fourier modes would contain Fourier modes up to 2N (cos(N θ) * cos(N θ) = 1 2 (cos(2N θ)+1)). Hence, if not enough Fourier modes are used in the resulting vector, the modes for n > N would pollute the ones for n < N . This is already taken into account in SFEMaNS: anti-aliasing is done by using N + N/2 modes when projecting. Still, this is not enough for D b as it contains two products. As such, it creates Fourier modes up to 3N and would require anti-aliasing to 2N when projecting. As simple as it sounds, it is not currently possible to do so because it requires the three fields used in the products to be projected at the same time in the physical space. Hence, there is currently an error of a few percent on the terms requiring double products of fields (Fig. 3.3).

. Transfer SFEMaNS fields to usable file formats

The files used by SFEMaNS are of a specific format that contains additional information in order to restart the temporal simulation. Moreover, all Fourier modes are saved in a complex manner depending on the mesh structure and its repartition on the processors. To make those files usable for statistical analysis and visualisation, a simple code has been developed to change the formatting of those files. Again those computations are done inside SFEMaNS and the project name is suites_to_bins in GitLab.

This project takes all the temporal iterations previously computed and outputs either binary files for statistical analysis using Python or pvd/vtu files for visualisation with Paraview. Once again, it is only needed to add the new options from 'data_example' into the data file used for the simulation and select the wanted options.

. Output file formats and name convention

The pvd/vtu files created for visualisation are named as follows:P2/P3 finite elements of

${field_name}_${iteration_number}_${filtering_number}$

The name for the binary files needs more specifications. Indeed, these files can be the representation of a field either in the Fourier or in the physical space. They are saved in a new directory called 'binaries_out'. For each field and each filter, a subdirectory is created with the name:

D${filter_number}_${field_name}

Depending if they are written in the physical or the Fourier space, the files will then have different names: phys_${field_name}_S${domain_number}_I${iteration_number} fourier_${field_name}_S${domain_number}_I${iteration_number} For vector fields, a number is added at the end of the field's name according to which component is written: (1, 2, 3) for (r. θ, z). For Fourier modes, a c or s is specified if the modes in the files are the cosine or sine ones. The domain's number is kept for simplicity's sake: it would take a lot more time to write the files on the total domain since all processors of each domain would wait for each other to write.

The files are structured such that the mesh points are in the same order in all of them. Inside the files, the points are written such as the 2D domain associated to the processor is repeated for each angle (or Fourier mode) and ordered in increasing order. This makes it possible to save the (r, z) positions of the points only once, the angles being deduced by knowing the number of Fourier modes used. The (r, z) positions as well as the statistical weight of each point are saved in files named:

${mesh}rr_S${domain_number} ${mesh}zz_S${domain_number} ${mesh}weight_S${domain_number}
where 'mesh' is the type of the mesh. As this code is currently only implemented for magnetohydrodynamics and we do not study the pressure field, it is either 'vv' for the velocity mesh or 'H' for the magnetic one.

Again, fields that do not have filtering have the filtering number set to 1.

. Temporal means and renormalization

Two other interesting features present in this code are the possibility to create a temporal average of a field or to renormalize a field in time. The renormalization is done by dividing a field by the average spatial value of the magnetic energy. It is very useful in order to study the mode responsible for the exponential growth of the magnetic field (see Chapter 4 section 5.1 of the article).

Making a temporal average will create a new field with an added 'mean' in front of its name and give it the iteration number '0'. The renormalization will just create a new temporal field with 'renorm' before the field's name. By combining both, it is possible to get the time average of a renormalized field whose new name will start with 'meanrenorm'.

. Selecting the fields to be processed

SFEMaNS post-processing outputs a lot of different fields and sometimes not all need to be visualised. As such, a dedicated file needs to be filled in to select the fields to study. This file is called 'data_for_suites' and, on each line, it contains:

${field_name} ${bool_1} ${bool_2} ${bool_3}

See an example in Fig. 3.4. Each boolean is either '.t.' or '.f.', the first is true if the field is on the velocity mesh, the second if the field is a vector, and the third if the field is filter dependent. The code will then only process the fields present in that file. Note that, if some renormalization or temporal average is done, 'field_name' needs to have 'renorm' or 'mean' in front.

. Visualisation

Visualisation of data is a key point for the analysis of physical phenomena by highlighting where and when they take place. All that process of transforming the data into usable pictures has been automatized because of the large quantity of fields and iterations to process. For that, we take advantage of Python's ability to run within Paraview, allowing it to operate without its graphical user interface (GUI). The fields need to be written in the pvd/vtu format in order to be read by Paraview. The project name in GitLab is par-aview_python. The Paraview version used is 5.8.0. There can be problems with some other versions as Paraview does not guarantee that the names of their Python functions are kept the same.

. Types of visualisation

As SFEMaNS is in cylindrical coordinates, three types of visualisation have been implemented (Fig. 3.5):

• 'animation_translation': cross sections at varying height z,

• 'animation_rotation': cross sections at varying angle θ,

• 'animation_DDD': isosurfaces in 3D.

As these animations only work for scalar fields, a desirable component needs to be specified for vector fields in 'component' (see the function call below). Also, it takes a Paraview source as an entry to avoid loading in memory multiple times the same field for different animations. As an example, here is one of the function's definition: def animation_rotation ( source , field , p , screen_path , screen_name , nb_frames , component = ' ' , id_colormap =0 , resolution = None , opt_field_tex = ' ' , opt_contour = ' ' , opt_view_norm = 1 . , opt_blade_source = None )

The necessary inputs are the Paraview source 'source', the name of the corresponding field 'field', 'p' which sets the maxima and minima of the color scale, the directory path for the pictures 'screen_path', the name for the output pictures 'screen_name' and the number of angles (resp. heights) designated by 'nb_frames' at which visualisations are performed for rotation/3D (resp. translation). There are a lot of other options. The main ones will be explained in the following as others are here for higher-level automation. The camera distance to the field is automatically rescaled in order to fit all of it into the frame by checking its maximal and minimal coordinates. However, it can still be scaled using 'opt_view_norm'. It is very useful when fields for the same simulation are on meshes of different dimensions. For example, in our case, the velocity mesh has a radial extension r ∈ [0, 1] whereas the magnetic mesh has r ∈ [0, 1.6]. By setting 'opt_view_norm' to 1.6 for fields that are on the magnetic mesh we can easily compare all quantities (see for example Fig. 3.6). ones but with the color scale being used for the perpendicular coordinate to the plane and arrows for the in-plane components. However, they did not get much use during the thesis as all energy transfers are scalar fields.

. Automation of visualisations

We have defined the visualisation functions. We now need to prepare a Paraview source. This is done to avoid loading the field multiple times and repeating computations on that field. For example, animations for vector fields need the computation of perpendicular components (r, θ). This is done in the 'make_animations' function in 'make_animations.py' which definition is: For basic usage, this function only needs the file path of the field to visualize (file name in 'file_names' and directory path in 'directory_path'). The field's characteristics need to be added beforehand to 'fields_dictionnaries.py' which requires its name in Tex and booleans to specify if it is vector, if it is on the magnetic mesh, and if it depends on the filtering. A plethora of parameters makes it possible to select different visualisations, the number of frames for each, the colormap, the number of isosurfaces, and so on.

However, further computation and combination of fields can also be done by giving multiple fields in 'file_names' and by indicating the computation to be done via 'computation_name'. Calculators are defined in 'calculators_dictionnaries.py' and require a lot of different pieces of information: their input and output type (vector (2), scalar (1) or same as input (0) ), how to construct the output field's name in both Paraview and Tex as well as how the calculator is written in Paraview. Here is an example of the pieces of information for calculating the z component of a field:

calculators_intype = { 'X ': [ 2 ] } calculators_outtype = { 'X ': 1 } calculators_name_constructor = { 'X ': [ ' ' , 0 , 'X '] } calculators_calc_constructor = { 'X ': [ ' ' , 0 , ' _X '] } calculators_tex_constructor = { 'X ': [ r '$ ' , 0 , r ' _x$ '] }
For the ones ending with 'constructor', the number 0 will be filled by the name of the field automatically at computation time. Sometimes, those calculators require fields from different meshes and thus a projection from the magnetic mesh to the velocity mesh has also been automatized in those cases.

Still, this function only does one field at a time. To visualise all fields in a directory, the function 'make_all_pvd_fields' is used. It is very similar to 'make_animations' but it only needs the directory's path containing all the fields as an input. Once again, various options make it possible to select the visualisations and their parameters. It is also possible to set the extrema of the colormap for each field via 'opt_ps' which is very useful when we do not want jumps in colorscale through multiple temporal iterations. An example of its usage is given in the project 'main.py' which is simplified below: from meta_func import * # directory where to put the . pvd and . vtu files data_directory = ' data_para ' # directory where to save the pictures screen_directory = ' screen_para ' # number of frames wanted for each animation nb_frames = 100 make_all_pvd_fields ( data_directory , screen_directory , nb_frames , tr = False , rot = False , DDD = True )

This code will output the 3D animations (for 100 angles of view) for all fields present in 'data_para' and write the pictures in the directory 'screen_para'.

As there can be a lot of fields to process, the computation can time out on supercomputer nodes by exceeding the time limit. In that case, an integer is written in the program's output file for each field done. The paths of those files need to be given as an argument when running this Python code to avoid doing the same field twice (see 'list' in the next section).

. Parallelization

The most efficient way to run this code on multiple processors is to give one field to visualize to each processor. This is done by specifying the rank and total number of processors in the first and second arguments of the Python call. An example of a complete command to run the code on multiple processors is given in 'job_para.slurm' in the GitLab project as well as simplified below: nb_tasks = 20 # number of paraview to run in parallel nb_proc_per_task = 2 # number of processeur given to each paraview for (( i = 0 ; i < nb_tasks ; i + + ) ) ; do srun -n1 --exclusive -c $nb_proc_per_task --unbuffered pvbatch --mpi --force -offscreen -rendering " $PARAVIEW_POST_PRO / main . py " " $i " $nb_tasks " $list " & done$ where 'nb_tasks' is the number of visualisation run in parallel and 'list' is the list of the previous output files. However, when using that kind of parallelization, one processor might not have enough RAM memory for a field, hence multiple processors need to be allocated to one Python call. This is done by using 'nb_proc_per_task'.

. Accessing the pictures

The quantity of pictures produced using the post-processing code can be very important. In our case, a run with around 20 different scalar fields on 120 iteration numbers with 3 different filters will output tens of thousands of frames only for 3D isosurfaces. In order to avoid transferring that huge amount of data (around 500G), a simple GUI has been implemented in Python. By just launching the code in the project 'image_viewer_sfemans' and giving it the path to the repository containing the pictures, it is possible to glimpse through all fields, iteration numbers, and filters (see Fig.

. Outlook

It was essential to develop these different programs so as to have an ironed-out flowchart from time simulation to final visualisation. They were also important to quickly test and view any modification implemented to the filter and to the computation of our various fields. It also made it very easy and fast to produce data for various internships.

Statistical analysis tools in Python have also been implemented early in the thesis but have since been completely revamped by Nikita Allaglo (L3 intern, spring 2023) and Rémi Bousquet (M2 intern, spring 2023). They have not since been used to do analysis for the subject of this thesis.

There is still room for a lot of improvements to this workflow. The main one is being able to run an instance of Paraview on multiple processors at once. This would help when multiple processors are allocated to one field to avoid memory overflow issues.

-Energy transfers in the von Kármán Sodium experiment

This chapter is dedicated to the analysis of the energy transfers in SFE-MaNS simulations of the von Kárman Sodium experiment. It presents itself in the form of an article currently submitted to Physics of Plasmas (August 2023). It contains the analytical demonstrations and equations for the energy transfers between fields and between scales for the magnetohydrodynamics equations for varying magnetic permeability and electrical conductivity. These energy transfers are then analyzed using two numerical simulations which generate a magnetic field but with different materials for the impellers (steel and soft iron). The conclusion is that, depending on the material of the impellers, the dynamo mechanisms have very different origins as well as spatial locations. Abstract Motivated by the observation that dynamo is a conversion mechanism between magnetic and kinetic energy, we develop a new approach to unravel dynamo mechanism based on local (in space, scale, and time) energy budget describing dissipation and scale-by-scale energy transfers. Our approach is based upon a new filtering approach that can be used effectively for any type of meshes, including unstructured ones. The corresponding formalism is very general and applies to any geometry or boundary conditions. We further discuss the interpretation of these energy transfers in the context of fast dynamo and anomalous dissipation. We apply it to results from direct numerical simulations of the Von Kármán Sodium setup using a finite element code, showing dynamo action for two types of impellers (steel or soft iron) in the magnetic field growth and saturation phases. Although the two types of dynamo hardly differ from the mean-field theory point of view (the velocity fields are the same in both cases), the locality of our formalism allows us to trace the origin of the differences between these two types of dynamo: for steel impellers, the dynamo is due to the transfer of velocity energy both in the bulk and in the vicinity of the impellers, whereas for soft iron impellers, the dynamo effect mainly comes from the rotation of the blades. We finally discuss possible signatures of precursors to anomalous dissipation and fast dynamo, that could become relevant in the inviscid limit. In many natural systems, such as the Earth or the Sun, a large-scale magnetic field can be generated by the so-called dynamo process. Such a process is described via the magnetohydrodynamics (MHD) equations coupling the equations modeling incompressible viscous fluids (Navier-Stokes equations) with the equations modeling electromagnetism (Maxwell equations in the quasi-static approximation).

Given the non-linearities of the equations and the enormous range of scales involved in the process, it is very difficult to identify the main mechanisms allowing the generation of a magnetic field by the motion of an electrically conducting fluid. From a theoretical point of view, exact analytical results are generally only obtained in very simple geometries or use symmetries to produce antidynamo conditions. Simplified models of dynamo generation can be built using mean field or stochastic approaches [1,2]. These approaches highlight the importance of helicity (via the so-called α-effect) and differential rotation (via the so-called Ω-effect) in the dynamo process. Such mean field approach provided several possible dynamo scenarios, depending on which combination is chosen [2]. However, confronting such simple scenarios with numerical or experimental observation turns out to be very difficult. From a numerical point of view, simulations of MHD equations are only possible for a limited range of parameters, excluding the very turbulent situations observed in the liquid core of the Earth or in the plasma making up stars like the Sun. Moreover, the large amount of data and information makes it difficult to identify the relevant mechanisms at work in the dynamo action. From an experimental point of view, the situation is the opposite: using liquid metal experiments, it is possible to reach very turbulent regimes. However, owing to the difficulty of measurements, it is only possible to produce sparse data that are not detailed enough to enable the identification of dynamo mechanisms, except in some special situations where the geometry is very constrained and the flow of liquid sodium reproduces analytical examples [3,4].

The von Kármán flow experiment using liquid sodium (hereafter VKS) is a good example of such difficulties. The facility uses 150 liters of liquid sodium that are set in turbulent motion by two counter-rotating impellers in a cylindrical container [5]. During its 16 years of operations, the VKS facility produced dynamo onset only with (i) impellers made of soft iron [5], and not with cooper or steel, (ii) for fluid motions with limited fluctuations [6]. This last observation rules out the possibility that the VKS dynamo is a mere rotating magnet. It also suggests that velocity fluctuations may impede the dynamo process, as already proved by dedicated numerical simulations in a simple geometry [7] or asymptotics [6,8]. Various attempts have been made to explain such results within the mean-field framework. For example, Fauve and Petrelis suggested a possible α -Ω mechanism for the dynamo, based on the combination of eddies generated near the impellers, and differential rotation [9,10], a scenario confirmed by numerical simulations of a simplified model of the velocity field [11] or RANS simulations of the von Kármán flow [12]. This model however does not explain the important influence of the iron impellers on the dynamo onset. Additional physical processes were then invoked such as the interplay between the Ω effect and a new conversion effect linked with spatially varying electrical conductivity [13,14,15] or magnetic permeability [16,17,18,19]. Finally, enhancement of the α effect via eddy collimation by the impeller magnetic field was discussed in [20,21].

Detailed verification of these hypotheses is difficult in the context of the mean-field model since it is a large-scale model where all the small-scale action is parameterized in the coefficients α and Ω. Here, we propose a new approach to unravel these mechanisms, motivated by the observation that dynamo is a conversion mechanism between magnetic and kinetic energy. The core of our approach is based on local (in space, scale and time) energy budget describing dissipation and scale-by-scale energy transfers.

Identifying the locations (in scale, space, and time) where kinetic and magnetic energies are exchanged supplies information about the location of the dynamo effect. Comparison with local magnetic or energy dissipation then tells us whether the conversion is sufficient to overcome damping, and thus, sustain dynamo action. The formalism is very general and applies to any geometry or boundary conditions. Then, by applying it to the case of the VKS dynamo, we show how it can be used to detect where and when the energy conversion between kinetic and magnetic energy takes place, thereby tracing the processes at work in the dynamo. Specifically, we use the massively parallel multiphysics code called SFEMaNS (for Spectral/Finite Element code for Maxwell and Navier-Stokes equations). This code solves the magnetohydrodynamics equations (MHD) describing the time evolution of the velocity and magnetic fields in the von Kármán sodium experiment. With this code two types of dynamo action have been explored: one using soft iron impellers and the other using steel impellers [22]. Although the two types of dynamo hardly differ from the mean-field theory point of view (the velocity fields are the same in both cases), the locality of our formalism allows us to trace the origin of the differences between these two types of dynamo.

The paper is organized as follows. The coupled equations describing the motion of a conducting fluid are presented in section §2. We first recall elements of basic mean-field theory, focusing on typical α and Ω dynamo mechanisms. In section §3, we present our new formalism, based on the derivation of local energy budgets from the MHD equations. We first build coarse-grained energy budgets including inertial dissipation, anomalous dissipation and transfer terms between velocity and magnetic fields.

Taking the limit of infinite resolution, we obtain exact energy budgets with no filter that are discussed. The generalization of the famous Kolmogorov 4/3rd law is presented in this context. This formalism is applied to the numerical simulations of the dynamo in the von Kármán geometry that are presented in sections §4.1-4.2-4.3. A crucial ingredient for application of the formalism is the notion of scale filtering. Traditional definitions of scale filtering, through convolution with a filter, are not suitable for unstructured grids. We, therefore, define a suitable filtering process that applies to any type of grid and boundary conditions, including the unstructured grid used in our numerical simulations in section §4.4. We discuss the standard dynamo mechanisms based on the mean-field theory in section §4.5. In section §5, we use the energy budget to unravel the dynamo mechanisms at work in the VKS experiment for impellers with different magnetic permeabilities. The same tools are used to study the saturation of the dynamo. This is done in section §6. Using these results, we discuss in section §7 the mechanisms for dynamo action in the VKS experiment and compare them with the literature. Our conclusion follows in section §8.

2 Background material and general formalism

Magnetohydrodynamics equations

The magnetohydrodynamics equations (MHD) for a neutral conductive incompressible Newtonian fluid in the quasi-static approximation are:

∂ t u + ∂ i (u i u) -ν∆u + ∇p ρ = - b × j ρ + f, (1)
∂ t b = ∇ × (u × b) -∇ × j σ , (2)
∇ • u = 0, (3)
∇ • b = 0, (4) j = ∇ × H, (5)
b = µH, ( 6
)
where H is the magnetic field, b the magnetic flux density (i.e. magnetic induction1 ), j the electric current density, u the fluid velocity (with u i one of its components with i ∈ {1, 2, 3}), p the fluid pressure, f the force per unit mass, µ the magnetic permeability, σ the electric conductivity, ν the kinematic viscosity and ρ the fluid density. The space and time dependence of all fields has been removed in the notations for convenience except when necessary.

Mean field theory: α and Ω effects

We briefly summarize classical mean-field concepts that are traditionally used to interpret numerical results and to understand the basic dynamo mechanisms. For the sake of simplicity, it is assumed here that σ and µ do not vary in space nor in time. Therefore equation ( 2) reads:

∂ t b = ∇ × (u × b) + 1 µσ ∇ 2 b. (7)
We consider b and u as superimpositions of mean and fluctuating parts such that b = b + b and u = u + u . Injecting these decompositions into equation ( 7) and taking its average, we get:

∂ t b = ∇ × ( u × b + emf) + 1 µσ ∇ 2 b with ∇ • b = 0 (8)
where emf = u × b is the mean electromotive force due to the fluctuations of motion and magnetic field. These equations together with proper initial and boundary conditions determine b if u and emf are given. The determination of the mean electromotive force as a function of u , b , u is the crucial step in the development of mean-field theory. For the sake of simplicity, we focus only on a low magnetic Prandtl situation νµσ 1, relevant to a fluid such as liquid metal. In such a case, we can consider that the dynamics of the magnetic fluctuations b is dominated by the interaction of the velocity fluctuations u with the mean magnetic field b only. This means that:

∂ t b ∇ × (u × b ) = ( b • ∇)u -(u • ∇) b . (9)
This fluctuating magnetic field is used to compute the electromotive force as detailed in [2]:

emf i = α ij b j + β ijk ∂ b j /∂x k ( 10
)
where the tensors α ij and β ijk are mean quantities determined by u . The α-effect is the occurrence of a mean electromotive force with a component pointing in the direction of the mean magnetic field b . Based on results from [12], we will focus on this effect and discard the β tensor. Approximating the fluctuations by a short in-time field (supposing a short turbulence correlation time), one can link the α tensor to the helicity tensor as

α ij = τ h ij = τ ε ikn u k (∇u) nj ( 11
)
where τ is the correlation time of the velocity perturbations, ε ikn is the Levi-Civita tensor and ∇u is the velocity gradient expressed in cylindrical coordinates (see appendix F).

Another conversion mechanism is due to differential rotation. Indeed, using the cylindrical coordinate system, the induction equation ( 7) can be written as:

∂ ∂t + u r ∂ ∂r + u θ r ∂ ∂θ + u z ∂ ∂z b r = b r ∂ ∂r + b θ r ∂ ∂θ +b z ∂ ∂z u r + 1 µσ ∆ * b r - 2 r 2 ∂b θ ∂θ (12)
∂ ∂t + u r ∂ ∂r - 1 r + u θ r ∂ ∂θ + u z ∂ ∂z b θ = b r r ∂ ∂r + b θ ∂ ∂θ + b z r ∂ ∂z u θ r + 1 µσ ∆ * b θ + 2 r 2 ∂b r ∂θ (13)
∂ ∂t + u r ∂ ∂r + u θ r ∂ ∂θ + u z ∂ ∂z b z = b r ∂ ∂r + b θ r ∂ ∂θ +b z ∂ ∂z u z + 1 µσ ∇ 2 b z (14)
with

∆ * = ∇ 2 -1 r 2 = ∂ 2 ∂r 2 + 1 r ∂ ∂r + 1 r 2 ∂ 2 ∂θ 2 + ∂ 2 ∂z 2 -1 r 2 .
This form of the induction equation clearly shows that, if the poloidal magnetic field, i.e. {b r , b z }, is non-zero, its shearing by the differential rotation (u θ /r) will always generate a toroidal magnetic field, i.e. b θ . This mechanism is called the Ω-effect. Denoting the two main terms in the rhs of (13):

Ω r = ∂ ∂r u θ r , Ω z = ∂ ∂z u θ r ,
one can see that the source terms for the toroidal field are b r rΩ r and b z rΩ z .

The dynamo effect is due to a combination of poloidal and toroidal magnetic field generation mechanisms, in the form of closing loops. Therefore different dynamo mechanisms can be conceived: there is a possibility for the poloidal magnetic field {b r , b z } to be generated from the toroidal field b θ through the α-effect. To close the loop and generate further toroidal field from such poloidal field, there are two possibilities (possibly concomitant): either through the α-effect or (and) the mean differential rotation rΩ z . The corresponding dynamo loops are then respectively called α 2 dynamo or α -Ω dynamo.

Local Energy budget

In order to evaluate the energy budgets describing relevant transfers of energy between the magnetic and velocity fields, we use a filtering approach and consider products of resolved and non-resolved quantities. This strategy is inspired by the mathematical derivation of anomalous dissipation due to the irregularities of (i) velocity field in the context of the Navier-Stokes equations [23] and (ii) velocity and magnetic fields in the context of the magnetohydrodynamics equations [24,25]. Our work extends these studies by introducing possible variations of the material properties (µ, σ). This is important in order to deal with the physical quantities resulting from the simulations presented in sections § 4.1-4.2-4.3 and to evaluate the impact of jumps in electrical conductivity σ and magnetic permeability µ.

Filtered equations

To characterise the energy transfers between scales of the same vector field, we use a filtering process in order to distinguish between the energy held by the different scales of the field. Given a vector field u and ∈ R + , we denote its filtered field u . At the present time, we do not specify the shape of the filter that is just assumed to obey the following properties: linearity, and such that u →0 --→ u for a smooth field u.

Applying the filter to the MHD equations ( 1) and ( 2), we get:

∂ t u + ∂ i u i u -ν∆u + ∇p ρ = - b × j ρ + f , (15)
∂ t b = ∇ × (u × b ) -∇ × ( j σ ) , (16)
∇ • u = 0, (17)
∇ • b = 0. ( 18
)
Because the density and the viscosity of the fluid are constant in our simulations, they do not impact the filtering and can be taken out of the convolution.

Local MHD coarse-grained energy budget at finite scale

Summing the scalar product of ρu/2 with equation ( 15) and the scalar product of ρu /2 with equation (1) and using analogous products for b, the equations for the filtered kinetic energy E c = ρuu 2 and filtered magnetic energy E m = bb 2µ can be respectively written as:

∂ t E c + ∇ • J N S + D ν + D u + D b = -T u→b + P, (19)
∂ t E m + ∇ • J M + D σ + D σ + D b + D µσ + D µT = T u→b + M µ , (20)
with:

J N S = ρ 2 ((uu )u -ν∇(uu )) + 1 2 (u p + up ), (21)
J M = 1 2µ { j σ × b + ( j σ ) × b + b × (u × b ) + b × (u × b)}, (22)
P = ρ 2 (u • f + u • f ), (23)
D ν = ρν∂ i u j ∂ i u j , (24)
D σ = j • j σ , D σ = ( j σ ) - (j ) σ • j 2 , (25)
D u = ρ 2 (uu i ∂ i u -uu i ∂ i u ), (26)
D b = 1 4 (u • (b × j) + u • (b × j ) -(u × b ) • j -(u × b) • j ), (27)
T u→b = 1 4 (u • (b × j) + u • (b × j ) + (u × b ) • j + (u × b) • j ), (28)
M µ = - b • b 2µ 2 ∂ t µ, (29)
D µσ = j 2σ • ∇ × (b ) µ -H , D µT = - (u × b) 2 • ∇ × (b ) µ -H . (30)
All terms involved in equations ( 19)-( 20) measure the transferred energy from all scales larger than to the scale . By convention, we do not use a subscript for all terms in contrast with [24,26,27].

The demonstration for (19) (respectively (20)) is in appendix A (respectively in appendix B).

We first discuss the different terms appearing in equation (19). On the left-hand side (lhs), J N S represents an energy current responsible for the spatial movement of energy and thus does not participate in any kind of dissipation or transfer of energy between the different scales. Integrated throughout space with the proper boundary conditions, it should disappear. D ν is the dissipation due to the kinematic viscosity: it corresponds to a loss of energy of the velocity field. D u is the anomalous/Duchon-Robert dissipation. This term is the local energy transfer at scale and thus characterises where the irregularities, and possible singularities, appear in the fields. It is called anomalous dissipation because, in the limit → 0, D u might not vanish [23]. These two different energy transfers (D ν and D u ) have already been derived and implemented in the SFEMaNS code (see section §4.2) in order to study their effects in a pure hydrodynamic problem modeled by the Navier-Stokes equations [28]. D b is a transfer term between different scales due to the interplay between velocity u and magnetic induction b fields. Like in the kinetic case, such a term could also produce an anomalous contribution to the magnetic energy budget in case the velocity field or the magnetic induction field are sufficiently irregular. If the contribution is positive, the net effect would then be an anomalous dissipation, like in the kinetic case. However, there is also the possibility that the contribution is negative, resulting in an increase of the magnetic energy due to irregularity. This would correspond to a scenario by which a dynamo is produced via a sufficiently rough velocity field [29], a mechanism known as fast dynamo mechanism [30].

The first right-hand side (rhs) term (-T u→b ) corresponds to the power of the Lorentz force and reflects an energy transformation between the fields u and b. It naturally appears with the opposite sign in the rhs of the magnetic energy budget (20) and is therefore responsible for the emergence and maintenance of the magnetic field. In other words, T u→b is a source term for the dynamo effect. The second term, P, displays the mechanical power injected by the forcing.

Aside the T u→b term, new interesting energy transfer terms appear in (20). In the lhs part, as before, J M is a purely local term that describes how magnetic energy is transported across the flow, and it should vanish after integration over space with the appropriate boundary conditions. D σ is the Joule dissipation, hence a loss of energy of the magnetic field, and can be used to trace, locally and across scales, the rate of resistive energy dissipation. D σ has contribution only where the electrical conductivity undergoes jumps, i.e. at {r = 1.4, 0 ≤ θ < 2π, -1 ≤ z ≤ 1} between the liquid sodium steady layer and the copper wall of the container (see fig. 4.1a). D µσ and D µT are two terms associated with the irregularity of the magnetic permeability µ which occurs only at the surface of the impellers. They are respectively denoted Joule induced in-between scales magnetic transfer and Lorentz induced in-between scales magnetic transfer. Lastly, in the rhs part, another source term for the dynamo effect appears, namely M µ , which we will call the magnet effect. It originates from the variation in time of the magnetic permeability and is a direct energy transfer from the impellers to the magnetic field.

Exact local MHD energy budget

The exact local energy budget can be obtained by taking the limit → 0. A technical difficulty may arise if the field is sufficiently irregular. In that case, terms corresponding to energy transfer may not need to vanish, producing anomalous dissipation of kinetic energy [23] or magnetic energy [31,32], or even dynamo action [29].

In the case where the field is sufficiently regular, so that these terms vanish, the time evolution for the kinetic and magnetic energies (per unit volume) simplifies to:

∂ t ρu 2 2 + ∇ • J N S 0 + D ν 0 = -T u→b 0 + P 0 , (31)
∂ t b 2 2µ + ∇ • J M 0 + D σ 0 = T u→b 0 + M µ 0 , (32)
with:

J N S 0 = ρ 2 (u 2 u -ν∇u 2 ) + up, (33)
J M 0 = 1 µ { j σ × b + b × (u × b)}, (34)
D ν 0 = ρν∂ i u j ∂ i u j , D σ 0 = j 2 σ , (35)
T u→b 0 = u • (b × j), M µ 0 = - b 2 2µ 2 ∂ t µ, P 0 = ρu • f. ( 36
)
As before, the J terms represent energy current responsible for the spatial displacement of energy and do not participate in transfers. The terms D ν 0 and D σ 0 are the standard viscous and resistive dissipation terms respectively. The right-hand side of (32) shows the two dynamo source terms, namely (minus) the Lorentz force power and the magnet effect and the rhs of (31) shows the Lorentz force power and the mechanical power due to the forcing.

Note that in the case where the velocity and/or magnetic induction fields are irregular, other terms should be considered in the energy budget, due to the non-vanishing of terms like D u , D b , D µσ and D µT in the limit → 0 in equations ( 19) and (20).

Definition of a local dynamo growth rate

We can use the energy budgets to define a global dynamo growth rate at several scales. Indeed the rhs of equation ( 20) evidences two possible source terms for the increase of the magnetic energy: the source term already discussed T u→b and the M µ term that exists only in case where permeability varies in time. 2 This last term is called in the sequel the "rotating magnet effect". Physically, the sum of these two terms must exceed the sum of all dissipative terms in the lhs of equation (20) for the dynamo to operate.

This can be formalized with the local growth rate λ m ( x, t) of the magnetic field at scale given by:

λ m = T u→b + M µ -D σ -D σ -D b -D µσ -D µT E m (37)
where . denotes the spatial average. Therefore the evolution of the filtered magnetic energy (20) can be recast into:

∂ t E m + ∇ • J M = λ m E m . ( 38
)
For → 0, the growth rate can have two types of behavior. For the first one, the field is sufficiently regular, in which case all transfer terms between scales vanish, and the local growth rate is simply:

λ m 0 = T u→b 0 + M µ 0 -D σ 0 E m . ( 39
)
This case will correspond to numerical simulations as no singularity can develop there because of the finite resolution. For the second type of behavior, in natural or experimental cases, there may be a contribution from the terms like D σ , D b , etc, resulting in a possibility of dynamo generation by singularity if these terms are negative. This is the case discussed in [29].

Connection with previous works

Our formalism can be connected to previous works by introducing δ u, the longitudinal velocity increment at scale , defined as

δ u = ˆ • (u(x + , t) -u(x, t)) (40)
where ˆ is the increment unit vector.

In [24], Galtier derives the local energy conservation in incompressible Hall MHD equations when the velocity and magnetic fields are possibly non-regular. He focuses on the time evolution of the total energy E c + E m for a plasma with constant electrical conductivity and magnetic permeability, corresponding to the sum of equations ( 19) and (20). As a result, the source term T u→b disappears together with all terms depending on the variation of σ and µ (namely M µ , D σ , D µσ , D µT ). The energy dissipation by viscous and resistive effects is then reduced to one term D ν + D σ . In the same way, the energy transfer term becomes lumped into mixed terms, involving both u and b, which does not allow separating the various source terms. The final energy budget can then be written using the notations of [24] (see equation (22) therein with the ion skin depth d I = 0) as:

D ε = 1 4 (∇ϕ ) • (δuδu 2 + δuδb 2 -2(δu • δb)δb)d r = S N S + SMag = 1 ρ D u + 2µD b + 1 2 ∇ • (J S M ag + J S N S ) (41)
where S N S (resp. SMag ) and J S N S (resp. J S M ag ) are defined in appendix C (resp. appendix D).

Therefore our analysis has a direct link with previous works [24,32], with the advantage that tracking E c and E m separately allows to get insights into the dynamo effect, as we now show.

4 Application to the VKS dynamo

Numerical setup

A standard setup for investigating turbulence in constrained geometry is the von Kármán swirling flow. The symmetries of this setup are presented in appendix F. It has been extensively studied using various fluids such as glycerol, water, air, helium gas and superfluid helium [26]. It was also used to study the dynamo action using liquid sodium in the so-called von Kármán Sodium experiment.

This experiment was successful in 2006 when a dominantly dipolar axisymmetric magnetic field with a toroidal magnetic field concentrated in the impellers was observed [5] for a magnetic Reynolds number larger than R c m = 34. We have performed MHD simulations of this setup on various supercomputers (up to about 2000 cores on the IDRIS supercomputer), allowing us to reach kinetic Reynolds number up to 10 5 and magnetic Reynolds number of the order of a few hundred, resulting in a magnetic Prandtl number around 10 -3 [22]. This is larger than the experimental value, which is of the order 5 10 -6 . Despite this, several key features of the magnetic field dynamo observed in the VKS experiment were recovered. The MHD simulations show the key role played by the ferromagnetic material (with high magnetic permeability) constituting the impellers. When soft iron impellers are modelled, the magnetic field is dominated by an axisymmetric axial dipole and azimuthal components located near the impellers in agreement with the experimental data (see figure 4.1b adapted from [22]).

SFEMaNS code

The SFEMaNS code integrates the system of coupled equations ( 1)-( 4). It uses a hybrid spatial discretization combining finite elements and spectral decomposition. For the hydrodynamic part, the approximation in space is done by using a Fourier decomposition in the azimuthal direction and the continuous Hood-Taylor Lagrange element P 1 -P 2 in the meridian section (quadratic approximation for the velocity field and linear approximation for pressure). All the discrete scalar functions f are written in the generic form:

f (r, θ, z, t) = f c,0 h (r, z, t) + M m F =1 f c,m F h (r, z, t) cos(m F θ) + M m F =1 f s,m F h (r, z, t) sin(m F θ), (42)
with (r, θ, z) the cylindrical coordinates, t the time and M the number of Fourier azimuthal modes considered. The functions f c,m F h and f s,m F h belong to a finite element space with h the typical mesh-size. The approximation in time is done by using a pressure-correction method [33]. The moving counterrotating impellers are accounted for by using a pseudo-penalty technique described in reference [34]. Modulo the computations of nonlinear terms with the fast Fourier transform, the linear problems at each time step for each Fourier mode in the meridian section are uncoupled and are thereby parallelized by using the message passing interface. The solution of each linear problem in the meridian section is further parallelized by using graph partitioning techniques from the METIS library [35] and subroutines from the portable extensible toolkit for scientific computation library (PETSc) [36], for the linear algebra. For the magnetic part, the algorithm solves the problem using the magnetic induction b in the conducting region (after standard elimination of the electric field) and the scalar magnetic potential in the insulating exterior. The fields in each region are approximated by using H 1 -conforming Lagrange elements with a penalty technique to control the divergence of b in a negative Sobolev norm that guarantees convergence under minimal regularity (see details in [37, §3.2], [38]). The coupling between conducting and insulating media is done by using an interior penalty method. SFEMaNS has been thoroughly validated on numerous manufactured solutions and against other MHD codes (see e.g. [39,40]).

Description of the numerical runs

In this paper, we focus on two series of runs described in [22] (Table 3). We use non-dimensional units such that the reference length L ref is set to the inner cylinder radius R cyl . The computational domain for the hydrodynamic study is

Ω = {(r, θ, z) ∈ [0, 1] × [0, 2π) × [-1, 1]}. The computational domain for the MHD study is the larger cylinder Ω ∪ Ω out with Ω out = {(r, θ, z) ∈ [1, 1.6] × [0, 2π) × [-1, 1]} (see schematic on figure 4.1a).
Denoting by σ 0 the electrical conductivity of the liquid sodium, ρ its density, and µ 0 the magnetic permeability of vacuum, the magnetic induction is made non-dimensional by using the Alfvén scaling B = U √ ρµ 0 , with U = ωR cyl where ω is the angular velocity of the impellers. The two governing parameters are R m = µ 0 σ 0 R 2 cyl ω, the magnetic Reynolds number, and R e = R 2 cyl ω/ν, the kinetic Reynolds number, with ν the kinematic viscosity of the fluid. We define a relative electrical conductivity σ r and a relative magnetic permeability µ r . These quantities are not constant since the walls and the impellers are made of different materials like copper, steel and soft iron. Specifically, we We use two series of runs described in detail in Table 3 of [22] at a moderate kinetic Reynolds number R e = 1500 but with a different relative magnetic permeability for the impellers: µ r = 1 for steel impellers, µ r = 50 for soft iron impellers. As indicated, for pseudo-vacuum boundary conditions (the tangential components of H are zero on the outer cylinder), the critical magnetic Reynolds numbers for dynamo action are R c m (µ r = 1) = 190 ± 10 and R c m (µ r = 50) = 90 ± 5. Therefore we will use R m values above the respective threshold, respectively R m (µ r = 1) = 300 and R m (µ r = 50) = 150 and in two regimes of the dynamo, i.e., in the growth and saturated phases. Each regime is sampled with 128 snapshots (every one-eighth of a turn to stroboscope the impellers' motion). We denote the average on all the snapshots as {.}.

The computations are performed with 128 Fourier modes and a non-uniform meridian mesh, with minimal spacing h min = 2.5 × 10 -3 and maximal spacing h max = 10 -2 leading to an averaged meridian mesh h mean = 7.92 × 10 -3 . We can estimate the Kolmogorov scale η = (ν 3 / ) 1/4 using the kinematic viscosity ν = 1/R e and the energy dissipation rate = 2R cyl K p /πH (with K p the time-averaged torque applied by the impellers to the fluid and H = 2 the total height of the container). Table 2 of [22] provides K p = 5.08 × 10 -2 at R e = 1500 resulting in η = 1.16 × 10 -2 .

Scale filtering

As the construction of local energy budgets relies heavily on the notion of filtering, we first present a general definition that is appropriate for any type of grid or boundary conditions. Classical definitions of scale filtering involve the convolution of a quantity with a filter of appropriate width (generally a Gaussian field). This definition is not well adapted to data on unstructured grids such as those used in SFEMaNS. We therefore consider another construction, that can be adapted to any type of grid, or boundary conditions. Given a vector field u and ∈ R + , we define the filtered field g = u as the solution of:

g -2 ∆g = u. ( 43
)
This definition ensures that it can be constructed for any type of grid, and is easy to implement inside a numerical program. Moreover, it is endowed with all properties of a classical filter: it is linear and, for a continuous field u, we have the limit u →0 --→ u. In fact, this filter corresponds to a convolution with a special filter. Indeed, in Fourier space, equation ( 43) reads:

g = ũ 1 + 2 k 2 , ( 44
)
where g is the Fourier transform of g. Thus, u is the convolution of u with the inverse Fourier transform of the Lorentzian function 1 1+ 2 k 2 .

Analysis using the mean-field framework

For later discussions, we first analyze the two different dynamos within the classical mean-field framework in the saturated regime of the dynamo. For this, we compute in each case the local helicity tensor using the velocity field averaged on the 128 snapshots. They are shown on figures 4.2-4.3. Their spatial distribution is complex and exhibits fine scales. Note that the dominant patterns of each local helicity component follow the symmetry rules detailed in appendix F (see equation ( 84)). This means that the large structures of the velocity field are mainly R π symmetric at R e = 1500 for both µ r values.

We can see that the maxima in absolute value are always localized near the impellers like in the RANS simulation of [12] whereas the minima are dispersed over the whole fluid domain. The maxima are associated with the swirling vortices attached to each blade and occupying part of the space between the blades. These vortices are thought to be a key ingredient of the dynamo mechanism [9,41,42].

The maximum absolute non-dimensional value of around 10 for both magnetic permeability values corresponds to the {h rr }, {h zz } components. These terms may be associated with an α-effect that would transform a poloidal field (b r , b z ) into a toroidal field b θ (first phase of the dynamo loop). To close the dynamo loop, we need the {h θθ } term that would transform the toroidal field into the poloidal field (second phase). This would lead to an α 2 dynamo mechanism, e.g.

(b r , b z ) {hrr},{hzz} -------→ b θ {h θθ } ---→ (b r , b z ).
Another possibility to transform a poloidal field into a toroidal field (first phase of the dynamo loop) is to resort to Ω terms as seen in equation ( 13). This would lead to an Ωα dynamo mechanism, These two dynamo mechanisms may be present in the VKS setup as it is underlined in [12,20,21]. Our results are not in contradiction with these explanations. However, as we now show, the local energy budget procedure allows for a much finer analysis of the dynamo mechanisms.

Dynamo growth mechanisms from local energy budgets

We now turn to the analysis of the local energy budgets and show how they can be applied to identify important dynamo mechanisms that are absent in the mean-field formulation adapted to statistically stationary problems. The filtering is computed at three different scales η = (0, 1, 16). The case η = 0 is used to ascertain the validity of the calculations, i.e. it is used to verify that D σ = D u = D b = D µσ = D µT = 0 in this case, where the fields are regular by construction. The other two filtering scales are chosen so as to identify the the mechanisms at small and large scale. Similarly, M µ = D µσ = D µT = 0 is expected for µ r = 1 (impellers with vacuum magnetic permeability such as the steel ones).

To proceed, we first analyze the spatial average of each quantity of the energy budgets (19), (20) as a function of time in the growing phase. We next describe the spatial repartition of each term of the energy budget averaged on a time range belonging to the exponential growth regime for the two types of dynamo action.

Time evolution of energy terms in the dynamo growth phase

Figure 5.1 shows that, in the growth regime of the dynamo, all transfer terms in the kinetic energy budget (19) are very similar for µ r = 50 and µ r = 1 and for each filtering scale. The reason is that the magnetic field is too weak to affect the velocity field. The only relevant terms are the viscous dissipation D ν and the anomalous dissipation D u (for > 0). D ν decreases with as expected since viscosity acts on small scales. In contrast, D u increases with because the larger as a result of energy cascade from large to small scale, see [26].

Figure 5.2 shows the exponential growth in time (in absolute value) of all terms involved in the magnetic energy budget (20) for the two magnetic permeabilities. In order to better discriminate the respective role of these terms, we renormalize them on figure 5.3 by the average magnetic energy E m (t).

We can note that the two dynamo regimes are different. For µ r = 1 and = 0, the growth rate is given by:

λ m 0 = T u→b 0 -D σ 0 E m . ( 45
)
Therefore T u→b 0 must exceed D σ 0 for the magnetic field to grow. The situation is similar for = η where the terms associated with irregularity such as D b and D σ are quasi-zero due to the fact that the fields are somehow regular at this scale. Increasing enhances the irregular transfer term D b (which is positive) and diminishes T u→b and D σ while D σ stays negligible. This may be an indication that in the inviscid limit, the dynamo could be sustained by irregularities of the velocity field, so that a fast dynamo could be achieved. However, in the present moderate Reynolds number regime, the dominant term is T u→b which compensates for all dissipative terms.

For µ r = 50, the first instants (t < 3) show the same behavior as that for µ r = 1. This means that T u→b is crucial in "launching" the dynamo process. After the launch is achieved, however, the dynamo trajectory becomes very different, as the rotating magnet effect M µ takes over and maintains the dynamo action, while T u→b is negative acting therefore like an anti-dynamo Lorentz power. Increasing enhances the modulus of the irregular terms: D µT relates the jump in µ and the electromotive field (u × b) and is positive; D b follows the variations of T u→b and is negative, i.e. it is an inverse cascade for the energy transfer between scales for both the velocity and magnetic fields. The terms D µσ (correlation between the jump in µ and the current) and D σ (dissipation due to the jump in σ) are negligible and can be ignored.

In summary, the term T u→b is essential to launch the dynamo in both cases, highlighting the crucial role played by the velocity field structure in the lightning of the dynamo. However, the further development of the dynamo is very different in the two cases, as M µ becomes the driving source term for the dynamo when using soft iron impellers, while T u→b becomes anti-dynamo. To further characterize these different mechanisms, we now consider the spatial localization of the local growth rate of the magnetic field.

T ub /<E m > D σ /<E m > M µ /<E m > D b /<E m > D σ _ /<E m > D µT /<E m > D µσ /<E m > 0

Spatial distribution of linear growth rate

We use equations ( 37) and ( 39) to determine the local dynamo growth rate. The spatial distributions are different depending on the value of the magnetic permeability: for µ r = 50, the growth rate is localized near the impellers for each while, for µ r = 1, the main localizations are near the impellers and the shear-layer. In this region, the growth rate is positive inside the bulk and negative near the lateral wall at r = 1 for each . Therefore the source of the growing magnetic field is very different depending on the impellers material.

Kinetic and magnetic energy distributions in the growth regime

The spatial distribution of the kinetic energy (fig. 5.5) is the same for both permeabilities for each scale because the magnetic field is too weak to impact the velocity field in the growth phase.

The linear mode of the magnetic energy exponential growth (fig. 5.6) is mostly axisymmetric and localized in the impellers' region for µ r = 50 while it is distributed in the bulk for µ r = 1. The structures are smoother with the increase of .

Energy dissipation in the growth regime

The viscous dissipation (fig. 5.7) is the same for both permeabilities and for each scale because the magnetic field is too weak to impact the velocity field in the growth phase. A small amount of dissipation occurs in the bulk at small scales but it is mainly concentrated near the impellers for all filtering scales.

The localization of the Joule dissipation (fig. 5.8) is different: for µ r = 50, {D σ / E m } is dominant in the impellers for all filtering scales . For = 16η, some negative regions appear near the lateral wall but the spatial average stays positive as indicated by figure 5.3; for µ r = 1, {D σ / E m } shows complex structures for all . The Joule dissipation in the bulk is mainly positive for small while it becomes negative for large . However, an interesting feature appears in the shear layer and stays positive for all . As shown in figure 5.9, the dissipation due to the irregularity of the current density {D σ / E m } is localized on the lateral wall where the conductivity is discontinuous, i.e., r = 1.4, 0 ≤ θ < 2π, -1 ≤ z ≤ 1. It is similar for both permeabilities for each scale and negligible.

Local energy transfers between scales in the growth regime

The anomalous dissipation (fig. 5.10) {D u / E m } is similar for both permeabilities at each scale in the growth stage. At = 0, it reaches the zero machine as expected. For = η and = 16η, it is mainly positive in the bulk where large scales provide energy to small scales and it is negative near the impellers where small scales feed large scales through the motion of the impellers blades. This behavior highlights the energy cascade at work in the von Kármán flow, where energy is injected at the blades, and transferred towards smaller scales in the bulk of the flow.

In contrast, the anomalous magnetic dissipation {D b / E m } (fig. 5.11) is very different for both per-meabilities. For µ r = 50, for all , it is localized near the impellers and it becomes more negative when increases. Its spatial average is negative in this case as shown in figure 5.3 (left) which makes it an inverse energy cascade.

For µ r = 1, blobs of alternate signs are localized near the impellers while larger areas of positive {D b / E m } are localized around the shear layer. The fact that {D b / E m } is negative (pro-dynamo sign) near the impellers and not in the shear-layer is in agreement with the mean-field arguments based on helicity [9].

Conversion of energy in the growth regime

We now examine the two source terms already identified. At the early stage t < 3, the (minus) power of the Lorentz force (fig. 5.12) {T u→b / E m } is positive for both permeabilities, but beyond, its sign and spatial repartition are very different depending on the value of µ r .

For µ r = 1, the shapes of {T u→b / E m } are similar for all : alternate signed zones near the impellers and positive areas in the bulk are visible. The latter contributes to a positive average which is compatible with the standard dynamo mechanisms linked with the presence of differential rotation and vortices in the shear layer and within the disks.

For µ r = 50 and all , we can observe two behaviors: at sufficiently small scales, it is positive within the shear layer, showing that the kinetic dynamo mechanism is still present. However, this mechanism is countered by an anti-dynamo mechanism occurring near the impellers, where the term alternates sign such that its average is negative. Another term is required for the dynamo to be sustained: the rotating magnet effect {M µ / E m }. As shown in figure 5.13, this magnet term is also only localized in the blades and disks with alternate signs, its averaged contribution being positive.

Dissipation and energy transfer through variation of magnetic permeability

Two dissipative terms are associated with the irregularity of the magnetic permeability µ r and exist only for µ r > 1: {D µσ / E m } (fig. 5.14) and {D µT / E m } (fig. 5.15). These two terms are therefore localized on the blades and disks on the impellers. They reach the zero machine as expected at = 0. For > 0, their shapes are similar for all . However, their magnitudes are different with {D µT / E m } being one order of magnitude greater than {D µσ / E m }. m

Saturation of the dynamo from local energy budgets

To determine the main mechanisms for dynamo saturation (when ∂ t E c = 0 and ∂ t E m = 0 on average), we now consider the local energy budget given by ( 19)-( 20) and evaluate the spatial average of all its components as a function of time. We want to determine what is the dominant balance. We next describe the spatial distribution of each term of the energy budget averaged on a time range belonging to the saturated regime for the two types of dynamo action.

Time evolution of energy terms in the dynamo saturated phase

Concerning the kinetic energy budget displayed in figure 6.1, for µ r = 1, D ν and D u follow the same behavior as in the growth regime displayed in figure 5.1: D ν is positive and decreases with , D u is positive and increases with . The transfer from the velocity field to the magnetic field is strong enough to be shown by the term -T u→b which is negative for all filtering scales. The irregular term D b is positive and acts as a direct energy cascade. Consequently, the velocity field maintains the magnetic field via the -T u→b term and against all dissipative terms.

In contrast, for µ r = 50 in the saturated regime, the term -T u→b is positive and D b is negative. Hence, -T u→b is a source term for the velocity field which extracts energy from the magnetic field and D b is an inverse cascade. The other terms, D ν , D u , are positive and behave like in the case µ r = 1. For µ r = 50, the time evolution of all terms of equation ( 20) follows that of the renormalized terms

M µ / E m , T u→b / E m , D µσ / E m , D b / E m , D µT / E m of

Spatial distribution of the local growth rate in the saturated regime

Note that, in the saturated regime, the space-average of the local growth rate in the saturated regime {λ m } is zero. However, we can look at its spatial distribution. For µ r = 50, the local growth rate is dominant in the neighborhood of the impellers. For µ r = 1, the localization is in three regions: the two impellers and the shear layer. Note that there is an asymmetry between the top and the bottom halves due to the event in the time range [970 : 990] when the magnetic energy decreases. This event induces a bias in the time average of {λ m }.

Kinetic and magnetic energy repartition

We first examine the kinetic energy (figure 6.4). One can note that, for = 0 and = η, the kinetic energy is distributed almost identically in space. That means that the kinetic energy has typical scales larger than the Kolmogorov scale. However, it presents more fine structures in the bulk for the µ r = 1 case than for the µ r = 50 case. The magnetic energy (figure 6.5) shows very different shapes depending on the relative magnetic permeability. For high µ r , it is concentrated near the impellers for all filtering scales and is mainly axisymmetric. In contrast, at µ r = 1, the magnetic energy spreads across the entire volume and presents non axisymmetric features. These results are coherent with the dominant m F = 0 magnetic mode for µ r = 50 and the mixing of m F = 0, 1 magnetic modes for µ r = 1 as indicated in Table 3 of reference [22].

Energy dissipations

These terms refer to the dissipation of energy from usual phenomena such as viscous dissipation and Joule heating.

Most of the viscous dissipation happens near the blades and the rim of the disks which correspond to the areas with the greatest gradient of velocity (see figure 6.6). This dissipation is strictly positive everywhere for all , it does not really change with µ r since the Reynolds number is identical.

In contrast, Joule dissipation depends strongly on µ r : when µ r = 50, for all , it occurs mostly in the impellers and along the cylinder axis, indicating that electric current is concentrated there. This corresponds to the spatial distribution of the current lines shown in figure 22 of [22]. For = 16η, a few negative spots appear between the blades.

When µ r = 1, Joule dissipation is distributed in three zones, the two impellers and the equatorial shear layer. For large scales ( η), the filter reveals regions of negative values corresponding to areas where the resistivity smoothes the field: it transfers energy from small scales to larger ones.

= 1.4, 0 ≤ θ < 2π, -1 ≤ z ≤ 1}
where the conductivity has a jump from the stagnant liquid sodium layer to the copper vessel. For all filter scales, it remains low (maximum absolute values around 10 -4 ) and is therefore not really relevant for the dynamo action.

Local energy transfers between scales

The most interesting feature of the filtering process is to highlight the energy transfer that happens in-between different scales via the D u and D b terms. These terms vanish in the limit → 0 in our simulation, but they may provide finite contributions for sufficiently rough velocity fields, resulting in anomalous dissipation and fast dynamo mechanism [23,24]. Figure 6.9 displays {D u }. This term has only a kinetic origin linked to the Navier-Stokes equations [28]. In the inviscid limit, it may contribute to the non-vanishing dissipation observed in the von Kármán flow [43,44]. In our simulation, this term is zero at = 0, as expected. For = η, it is mainly localized near the impellers. For = 16η, more negative regions appear. They correspond to areas where the scales smaller than 16η inject energy into the system. They are found near the impellers but also near the lateral wall, highligthing the non-homogeneity of the cascade and the importance of the boundaries in the transfer mechanism. Due to the presence of the magnetic induction, a new transfer channel is expected allowing transfers from kinetic to magnetic energy through the term {D b }, that can provide a fast dynamo mechanism in the inviscid limit [30], for rough enough velocity fields. As seen in figure 6.10, this term is dominantly localized near the impellers for µ r = 50. Patterns with alternate signs are stuck to the blades but their spatial average is negative as shown in figure 6.2 (bottom left column panel). That means that an inverse cascade involving the velocity and magnetic fields exists from small scales to large scales near the impellers.

Similarly, for µ r = 1, spots of alternate signs are close to the impellers, but complex positive structures also form close to the shear layer as increases. These positive regions are dominant and contribute to an overall positive spatial average, as shown in figure 6.2 (bottom right column panel). In this case, the cascade involving the velocity and magnetic fields is direct from large scales to small scales.

Conversion of energy

The conversion of kinetic energy to magnetic energy relies on the T u→b and M µ terms (as seen in equation ( 20)). The former is the standard source term: when it is on average positive, the velocity field provides energy to the magnetic field; the opposite occurs when it is negative. Figure 6.11 shows that the spatial variation of {T u→b } depends on the value of µ r . For µ r = 50, it is localized near the impellers with alternate signs. Its spatial average is negative in the saturated regime as shown in figure 6.2 (left column). It corresponds to a loss of magnetic energy. For µ r = 1, there are more transfers in the bulk and smaller values in the impellers. Note that it is the only source term for the maintenance of a magnetic field when µ r = 1. Therefore it has to be positive on average for a magnetic field to saturate as it is seen in figure 6.2 (right column).

In contrast, for µ r = 50, {M µ } is another source of magnetic field increase. This term only exists where the blades rotate as expected, and it clearly follows the azimuthal variation of the blades with m F = 8. It is on average positive and greater than the negative average of {T u→b } (see figure 6.2 left column). This allows qualifying this dynamo as a magnet effect in the saturated regime. However, at the beginning of the linear regime, it is dominated by the {T u→b } term as shown in figure 5.3 for t < 3.

Dissipation and energy transfer through variation of magnetic permeability

The variation of µ r due to the impellers is responsible for the terms {D µσ } and {D µT }. These are localized at the surface of the blades and disks. They exist only for > 0 and show alternate signs. They are specific energy transfers between scales linked to jumps in magnetic permeability. The first one is induced by the Joule dissipation while the second one is induced by the Lorentz force power.

Mechanisms for dynamo action

Inspection of the linear regime has proved that the standard mechanism for the growth of the dynamo acts for both values of the impellers permeability at early instants: as soon as the threshold in magnetic Reynolds number is exceeded, the (minus) Lorentz force power overrides the dissipative terms. For steel impellers (µ r = 1), at larger times, a balance between T u→b (which is in space-average positive) and the dissipative terms (whose space-average values are positive) is reached and allows a saturated regime to be achieved. In contrast, for soft iron impellers (µ r = 50), the rotating magnet effect due to the motion of µ r and the anomalous magnetic dissipation (which is negative on space-average) take over, and become the key players for the dynamo mechanism in the saturation phase. In light of our findings, we can now discuss the different scenarios proposed in the literature to explain the VKS dynamo.

The mechanism based upon spatial variations in electrical conductivity has first been proposed by [13]. The authors have considered the magnetic field generation by an electrically conducting fluid flowing parallel to a rigid plate presenting spatial variations in the streamwise direction. It was shown that the magnetic Reynolds number for dynamo process decreases with increasing conductivity of the plate and with increasing amplitude of its modulation. The growth mechanism is the coupling between poloidal and toroidal components of the magnetic field through the high conductivity zones. This mechanism was also recovered in [14,15]. Spatial variations in magnetic permeability have been already put forward as a source of dynamo action [16,17,18,19]. In [16], by combining the effects of the soft-iron impellers and a small uniform homogeneous α θθ -effect modeling the induction effects of unresolved small-scale flow fluctuations, an axisymmetric magnetic field mode was obtained in a kinematic dynamo setup. However, a uniform homogeneous α-effect is not realistic as shown in figure 4.2. In [17], the growth of the magnetic field is due to the shear of the velocity localized on the frontier between fluid and the top of the blades. The jump conditions of section 2.2 (iii) and (iv) in [17] can be interpreted as a M µ effect localized only on the blades' top. In our case, the helical flow in between the blades distorts the magnetic field and bends it along the surface of the blades and disks making the M µ term large there (as shown in figure 6.12). Therefore it is not the dynamo mechanism described in [17] as it would occur only on the top of the blades. In [18], the solid body rotation of an anisotropic conductive or an anisotropic permeable cylindrical rotor (the anisotropy of material properties follows a logarithmic spiral shape) immersed in a static cylinder can generate an axisymmetric magnetic field. The kinematic dynamo mechanism there is due to differential rotation coupled with anisotropic diffusion. However, the generated magnetic field has a horizontal structure following spirals and a Fourier vertical structure as exp(ikz) with k the vertical wavenumber of the corresponding eigenmode, far from the VKS dynamo mode shown in figure 4.1 (right). In [19], the effect of the rotating disk is studied: by adding opposite α θθ and Ω effect above and under the disk with an higher magnetic permeability than the liquid, a decrease of the magnetic threshold is observed. However, an α θθ -effect with opposite sign above and under the disk is not what is observed in figure 4.2. Moreover, any α-or-Ω-driven dynamo mechanism would only participate through T u→b which is globally negative in the current case µ r = 50.

Here, the main dynamo mechanism at large times for µ r = 50 is simply a self-amplifying phenomenon due to the variation of the magnetic permeability in time. Re-writing (20) by keeping only M µ gives:

∂ t E m ∼ E m ∂ t 1 µ . (46)
This implies that the magnetic field grows (resp. decreases) behind (resp. in front of) the blades as can be seen in figure 6.12, which amplifies the phenomenon: stronger M µ happens where the magnetic field is amplified. Hence, for soft-iron impellers (µ r = 50), the dynamo effect is due to rotating magnets.

In the saturated regime, this is counterbalanced by both the Joule dissipation, D σ , and the energy transfer toward the velocity field, -T u→b , (see figure 6.2).

The dynamo mechanism for steel impellers (µ r = 1) is very different as it is only due to the energy transfer from the velocity field. At the same time, it couples dynamo processes in the shear layer and the in-between blades zones (see figure 5.12).

Conclusion

The main result of this article is the expression on the different local transfer terms appearing in equations ( 19)- (20). They express all possible exchanges between velocity and magnetic fields. Obtaining these terms was made possible by the development of a new filtering approach that can be used effectively for unstructured grids such as the finite element mesh used in SFEMaNS. We believe that this filtering method offers a general framework for the analysis of data obtained on unstructured meshes.

The viscous D ν and resistive D σ dissipations correspond to the classical terms. The lack of smoothness of the solutions leads to the appearance of D b that generalizes the anomalous dissipation D u [23] (which quantifies the kinetic energy dissipated by non-viscous means) to the anomalous magnetic energy dissipation. It is rather remarkable that the combination of these anomalous dissipative terms can be written in the form of structure functions already derived in the literature [24,26] (as shown in appendices C-D).

The novelty is that the separation of the kinetic and magnetic energy equations allows highlighting the main terms responsible for the dynamo effect, namely T u→b and M µ . Using two illustrative numerical dynamo simulations with different impellers materials, we were able to show that the dynamo mechanisms are totally different when using steel or soft iron impellers. In the case of steel impellers, the only source term is the (minus) Lorentz power T u→b that needs to compensate for all dissipative terms. In the case of soft iron impellers, the ferromagnetic nature of the impellers plays a crucial role in maintaining the magnetic field through M µ (and D b for a filtering scale > 0). This second mechanism cannot be described by a mean-field approach, which does not include variations in magnetic permeability nor electrical conductivity.

MHD and Navier-Stokes equations share a lot in common. However, from a mathematical perspective, it is unknown whether the viscous and resistive mechanisms that tend to smooth out potential irregularities in the velocity and magnetic field are effective enough to constrain the fields to stay smooth at all times. The analytical tools we have developed will allow us to dynamically study the rare events associated with extreme values of the anomalous dissipative terms in future work.

Acknowledgement

The HPC resources for SFEMaNS were provided by GENCI-IDRIS (grant 2022-0254) in France.

space with the proper boundary conditions, the latter should disappear but it can be responsible for localized large values on the boundaries. If we are interested in the energy injected in the small scales, the important term is D u .

On the other hand:

4D b = -u • (j × b -j × b) -j • (u × b -u × b) = u • (b × j ) + j • b × u + j • b × u + u • (b × j) (71)
Hence, by regrouping all terms with a divergence in a current, and using ∇ × b = µj and ∇ × b = µj (because µ = const.):

2µD b = SMag - 1 2 ∇ • J S M ag (72)
with:

J S M ag = b ×(b×u)+b×(b × u )-(u • b)b -b(u • b )-b(b•u )-b(u•b )+ 1 2 ub 2 + 1 2 u b 2 + 1 2 u(b •b) (73)
And so, in this special case :

S M ag = SMag - 1 2 ∇ • J S M ag (74)

E Filtering issues

We examine the issues related to the filtering process concerning the boundaries and its non-commutation with derivatives. We start with equation (??) in a domain D ⊂ R 3 with boundary ∂D = ∂D:

(1 -2 ∆)f i = f i (75)
The weak formulation is obtained using test functions φ i :

D [(1 -2 ∆)f i ]φ i = D f i φ i ⇔ D f i φ i + D 2 (∇ k f i )(∇ k φ i ) - ∂D 2 n k (∇ k f i )φ i = D f i φ i
This is the formulation used in SFEMANS to solve these equations.

E.1 Discontinuous Galerkin method

This section explains how to solve the filtering equation at the interface between the velocity mesh and the magnetic field mesh (that includes the velocity mesh). This is needed in order for both meshes to communicate, otherwise the filter would create a discontinuity in the filtered field at the interface between the meshes.

Let us consider a domain D ⊂ R 3 with boundary ∂D. The domain is assumed to be partitioned in two subdomains

D 1 and D 2 such that D = D 1 ∪ D 2 , D 1 ∩ D 2 = ∅. The interface between both domains is denoted by Σ = D 1 ∩ D 2 . D 1 ∪D 2 [(1 -2 ∆)f i ]φ i = D 1 ∪D 2 f i φ i ⇔ D 1 ∪D 2 f i φ i + D 1 ∪D 2 2 (∇ k f i )(∇ k φ i ) - ∂D 2 n k (∇ k f i )φ i - Σ 2 n 1k (∇ k f 1i )φ 1i + 2 n 2k (∇ k f 2i )φ 2i = D 1 ∪D 2 f i φ i
The vector n i is the outward normal vector to the interface Σ with respect to D i .

Using

(a 1 b 1 + a 2 b 2 ) = 1 2 [(a 1 + a 2 )(b 1 + b 2 ) + (a 1 -a 2 )(b 1 -b 2 )] : Σ n 1k (∇ k f 1i )φ 1i + n 2k (∇ k f 2i )φ 2i = Σ 1 2 (n 1k ∇ k f 1i + n 2k ∇ k f 2i )(φ 1i + φ 2i ) + (n 1k ∇ k f 1i -n 2k ∇ k f 2i )(φ 1i -φ 2i ) = Σ 1 2 n 1k ∇ k f i {φ i } + n 1k {∇ k f i } φ i using n 2k = -n 1k and defining f i = (f 1i -f 2i ) and {f i } = (f 1i + f 2i ).
In cylindrical coordinates, we have:

n k (∇ k f i )φ i = (n r ∂ r + n z ∂ z )f r (n r ∂ r + n z ∂ z )f θ (n r ∂ r + n z ∂ z )f z • φ r φ θ φ z
For regularization purposes, the following term is added:

Σ β h (f 1i -f 2i )(φ 1i -φ 2i ) = Σ β h f i φ i
where h represents the mesh size and β a tunable penalty parameter that is always set to one in our numerical investigations.

The final implemented equation is:

D 1 ∪D 2 f i φ i + D 1 ∪D 2 2 (∇ k f i )(∇ k φ i ) - ∂D 2 n k (∇ k f i )φ i - Σ 2 2 n 1k ∇ k f i {φ i } + n 1k {∇ k f i } φ i + Σ β h f i φ i = D 1 ∪D 2 f i φ i (76)
Note that, when = 0, we get f i 0 = f i as needed. When = 0, the impact of the boundaries (∂D, Σ)

scales like 2 .

E.2 Non-commutation of filter and derivatives

In this section we study commutation between filtering and derivatives. As an example, we study the case where we derive according to one direction X. The equation to filter such a derivative would be:

(1 -2 ∆)∂ X f i = ∂ X f i , (77)
with its weak formulation being:

D [(1 -2 ∆)∂ X f i ]φ i = D (∂ X f i )φ i . (78)
Working on the rhs gives:

D (∂ X f i )φ i = - D f i (∂ X φ i ) + ∂D f i φ i = - D [(1 -2 ∆)f i ]∂ X φ i + ∂D f i φ i = D [∂ X (1 -2 ∆)f i ]φ i - ∂D [(1 -2 ∆)f i ]φ i + ∂D f i φ i = D [(1 -2 ∆)∂ X f i ]φ i + ∂D 2 (∆f i )φ i + ∂D (f i -f i )φ i Hence: D [(1 -2 ∆)(∂ X f i -∂ X f i )]φ i = ∂D (f i -f i )φ i + ∂D 2 (∆f i )φ i . (79)
When = 0, the filter and derivatives commute if Dirichlet boundary conditions are used (vanishing of the first rhs term). In contrast, when = 0, there are errors which scale like 2 on the boundaries ∂D as indicated by the second rhs term in 2 . This means that, as long as we are interested in derivatives away from the edges ≥ 2 , the error introduced by commutating derivatives and filtering is negligible.

F Symmetries of the von Kármán flow and implications for the Ω and α terms

The von Kármán flow with exactly counter-rotating impellers has one particular symmetry that can provide information about the spatial distribution of instantaneous and averaged quantities based on the velocity field. It is called the R π symmetry [45] and corresponds to a symmetry of rotation of π about any horizontal axis:

R π   u r u θ u z   (r, θ, z) =   u r -u θ -u z   (r, -θ, -z) (80)
where u = (u r , u θ , u z ) is the velocity vector in cylindrical coordinates with r the radial distance, θ the azimuthal angle and z the height. Another symmetry, apart from the impellers, is axisymmetry about the z-axis defined as:

S θ 0   u r u θ u z   (r, θ, z) =   u r u θ u z   (r, θ + θ 0 , z). (81)
Based on [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], one can deduce that, for a velocity field symmetric under the R π symmetry, the Ω-terms satisfy: Ω r (r, θ, z) = -Ω r (r, -θ, -z) and Ω z (r, θ, z) = +Ω z (r, -θ, -z)

For the α-tensor (11), one needs the expression for the velocity gradient tensor in cylindrical coordinates:

Using (83) in the helicity tensor, one can infer that, for a velocity field symmetric under the R π symmetry, the helicity tensor satisfies:

  h rr h rθ h rz h θr h θθ h θz h zr h zθ h zz   (r, θ, z) =   +h rr -h rθ -h rz -h θr +h θθ +h θz -h zr +h zθ +h zz   (r, -θ, -z). (84)
The α-tensor will follow the same rules.

-Conclusion

In this thesis, the goal was to better understand the dynamo effect in turbulent conducting flows. For this purpose, we simulated the von Kármán Sodium configuration using the SFEMaNS code. The choice of this configuration was motivated by the possibility of comparing the numerical results with the ones of a real experiment carried out at CEA Cadarache. To shed more light on the dynamo mechanism, we decided to study energy transfers between magnetic and velocity fields, as well as between scales. This is an extension to the magnetohydrodynamics equations of the approach used in a previous thesis by H.Faller [START_REF] Guermond | SFEMaNS documentation[END_REF]. These choices led to the need to implement various processing tools as well as updates to the SFEMaNS code.

A major improvement we implemented to the SFEMaNS code was the distribution of the simulation meshgrid across the processors. This has saved an enormous amount of memory, reaching around 300,000 grid points per processor in Navier-Stokes simulations which makes it possible to simulate higher kinetic Reynolds numbers. Thanks to this, we plan in the near future to compare DNS and LES computations performed at high Reynolds numbers to see the impact of the turbulence model.

In parallel, we developed a new method based on preconditioning the Stokes problem in SFEMaNS for potential future use in a time-stepping technique and we compared it against previously published preconditioning methods. The various preconditioning techniques have achieved similar performance and good scaling for parallel computing. However, they are not powerful enough to compete with the prediction-correction method currently in use in the SFEMaNS code. The limiting factor in these preconditioning methods is the time needed to invert the velocity problem. Hence, we see no reason to use these methods until further theoretical work has been carried out to make them more efficient.

In order to efficiently calculate and process the quantities of interest, we implemented numerous post-processing tools as there was no automation for the analysis and visualisation of SFEMaNS' outputs. Particular attention has been paid to saving memory space and processing time when calculating these quantities. A general tool to transfer any output of SFEMaNS into simple file formats for statistical treatment or visualization has been created. An automatized Paraview workflow with multiple functionalities has been implemented to easily visualize structures in 3D or slices. All these tools were necessary to implement, given the quantity of fields to be analyzed. They are also independent of each other and can therefore be used for any future simulations using SFEMaNS. The only tools that remain to be developed are those of statistical analysis, such as high-order moments and structure functions.

The derivation of energy balances in magnetohydrodynamics has led to the discovery of a variety of energy transfers. Above all, the influence of magnetic permeability variations on energy transfer at different scales is now well described. In our two case studies, the observed dynamo mechanism is highly dependent on the magnetic permeability of the turbines. In the case of steel impellers, the Lorentz force power of the fluid manages to compensate for all dissipation terms, whereas in the case of soft-iron turbines, it becomes anti-dynamo. For soft-iron impellers, in the saturated regime, rotation of the high-permeability impellers becomes the sole source of dynamo. However, in both cases, the ignition of the dynamo is purely due to the nature of the flow as the Lorentz force power is the only contributor.

A more complete analysis of the various energy transfers can be pursued. An analysis at a higher kinetic Reynolds number of the energy transfers responsible for anomalous dissipations may lead to traces of quasi-singularities or fast dynamo mechanisms. Simulations for other magnetic permeabilities of the impellers should provide a better understanding of the balance between the Lorentz force power and the rotating magnet effect. For all these simulations, it should be possible to track rare events and interesting coherent structures thanks to the implementation of all visualization tools. Finally, a statistical treatment of these quantities should be carried out to highlight correlations enabling a better understanding of their interactions. L'effet dynamo est un processus permettant de générer un champ magnétique grâce au mouvement de matière conductrice de l'électricité. Pour les planètes, il s'agit d'un liquide (du fer dans le coeur liquide de la Terre) mis en mouvement par de nombreux mécanismes (effet Coriolis, convection thermique et solutale, etc) et, pour les étoiles, de gaz ionisé. Dans tous les cas, l'effet dynamo requiert la présence de matière conductrice et une source d'énergie la maintenant en mouvement. La génération d'un champ magnétique a, par ailleurs, été reproduite en laboratoire en utilisant du sodium liquide mis en mouvement grâce à différents dispositifs [1,2,3]. Cependant la seule expérience ayant reproduit une telle génération sans contraintes fortes sur le mouvement du fluide est celle appelée von Kármán Sodium (VKS) [4], les expériences précédentes ayant utilisé un champ de vitesse mathématiquement prédéterminé pour générer de la dynamo.

-Résumé en français

L'expérience de von Kármán Sodium est une installation impressionnante mise en place grâce à la collaboration de l'ENS Paris, l'ENS Lyon, du CEA Saclay et du CEA Cadarache. L'objectif de cette expérience est de créer et maintenir un champ magnétique grâce au transfert d'énergie provenant du fluide conducteur, et ainsi reproduire l'effet dynamo présent à l'intérieur des planètes et des étoiles. Le fluide utilisé est du sodium, maintenu à 120°C, et mis en mouvement turbulent grâce à deux turbines en contra-rotation.

Le maintien d'un champ magnétique a été mesuré pour la première fois en 2006. Le champ a été généré seulement lors de l'utilisation de turbines en fer doux [4] (et non en acier ou en cuivre) et avec un écoulement contenant peu de fluctuations [3]. Ce dernier point prouve alors que le champ magnétique généré ne provient pas uniquement de la rotation d'un aimant dans un milieu conducteur. Cependant, cela semble indiquer que les fluctuations nuisent aux mécanismes dynamos comme déjà remarqué dans des simulations numériques pour des géométries plus simples [5] ou en régime asymptotique [6]. Le champ obtenu est alors en moyenne dipolaire, aligné avec l'axe de symétrie, avec une composante azimutale concentrée au niveau des turbines. Expérimentalement il n'est pas possible la vitesse du fluide ou les variations non axi-symétriques du champ magnétique. Mais, pour certains paramètres, l'inversion du champ magnétique (Nord-Sud vers Sud-Nord et inversement) a été observée comme c'est le cas pour le champ magnétique terrestre.

La principale difficulté pour trouver les conditions permettant d'obtenir l'effet dynamo provient de la nature turbulente des écoulements des fluides dans les étoiles et les planètes. Un écoulement turbulent correspond à un champ de vitesse désordonné variant rapidement et contenant de nombreux vortex. Ce type d'écoulement apparaît lorsque les forces visqueuses (correspondant à la capacité pour le fluide de dissiper son mouvement) sont faibles, en comparaison à la force inertielle du système. Cela génère alors un champ de vitesse caractérisé par de nombreuses échelles spatiales : du mouvement global généré par les forces externes aux plus petits tourbillons dissipant l'énergie cinétique. L'effet dynamo provenant de cet écoulement, il devient alors difficile de trouver quels évènements provoquent la génération d'un champ magnétique. En particulier, cela peut entraîner la génération d'un champ magnétique par un champ de vitesse très irrégulier [7] via des mécanismes de dynamos rapides [8,9]. Cela arrive quand le champ magnétique croit exponentiellement en un temps caractéristique de l'ordre de grandeur du temps advectif de l'écoulement. Ce phénomène est provoqué par des écoulements très turbulents pour lesquels les lignes de champ magnétique sont étirées de manière exponentielle par les vortex.

Pour pallier ces difficultés, une approche possible est de construire des modèles de dynamo simples en utilisant les théories de champs moyens et la théorie stochastique [10,11]. Ces approches révèlent l'importance de deux éléments de l'écoulement du fluide : son hélicité (à travers l'effet α) et sa rotation différentielle (à travers l'effet Ω). Pour l'expérience VKS, de nombreuses combinaisons de ces effets ont été considérées pour expliquer l'émergence du champ magnétique [12,13]. Cependant, ces modèles ne peuvent qu'analyser l'impact du champ de vitesse sur le champ magnétique. Il leur est donc impossible d'expliquer les différences de comportement observées dans l'expérience pour les différents matériaux utilisés pour les turbines. Cela requiert alors l'introduction de modélisations supplémentaires pour les interactions entre l'effet Ω et les variations spatiales ou temporelles de la conductivité électrique [14,15,16] ou de la perméabilité magnétique [17,18,19,20]. D'autres effets d'amplification de l'effet α par collimation proche des turbines ont aussi été analysés [21,22].

De plus, la théorie des champs moyens rend difficile la compréhension de la phase initiale de croissance du champ magnétique car elle nécessite la distinction entre un champ moyen et un champ rapidement variable pour la vitesse et le champ magnétique. Finalement, comme tous les évènements de petites échelles se retrouvent paramétrés dans les coefficients des effets α et Ω, il n'est pas possible d'observer les structures potentiellement responsables pour le mécanisme de dynamo rapide.

Dans cette thèse, nous proposons une nouvelle méthode pour comprendre les mécanismes dynamo permettant de générer un champ magnétique dans les fluides conducteurs et turbulents. L'idée principale est d'étudier les bilans d'énergie locaux (temporellement et spatialement) en décrivant toutes les dissipations et transferts d'énergie entre champs et échelles. Cette approche a été proposée par Duchon et Robert [23] pour l'étude des écoulements turbulents décrits par les équations de Navier-Stokes. Cette méthode se base sur un principe de filtrage des champs permettant de séparer les différentes échelles de l'écoulement. Comparée à d'autres techniques de filtrage [24], celle-ci permet de détecter de potentielles singularités grâce à l'introduction de dissipations anormales correspondant à des transferts d'énergie vers des échelles de plus en plus petites. Cette technique a ensuite été étendue à la magnétohydrodynamique [25,26]. Cependant, dans ces références, les variations de perméabilité magnétique et de conductivité électrique ne sont pas prises en compte alors qu'elles sont cruciales pour notre étude de l'expérience VKS et dans nos simulations. De plus, comme ces articles sont spécialisés dans l'étude de phénomènes astrophysiques (en particulier les vents solaires), les énergies magnétique et cinétique sont combinées. Comme nous souhaitons étudier les transferts d'énergie entre le champ de vitesse et le champ magnétique, nous devons revisiter cette approche pour séparer les deux. Ainsi, dans cette thèse, nous étendons cette méthode pour la magnétohydrodynamique dans le cas de conductivité électrique et de perméabilité magnétique variables (spatialement et temporellement) en décrivant tous les transferts d'énergie (d'un champ à l'autre et d'une échelle à l'autre). En combinant des simulations numériques déjà produites pour VKS avec cette approche de transferts locaux en énergie, il devient possible de déterminer où et quand les effets dynamos se produisent.

. Magnétohydrodynamique et effet dynamo

Pour commencer, les équations de la magnétohydrodynamique (MHD) sont utilisées pour modéliser les interactions entre le fluide conducteur et le champ magnétique. Ces équations proviennent de la combinaison des équations de Maxwell prises dans la limite non-relativiste, couplées aux équations de Navier-Stokes à travers la force de Lorentz et la loi d'Ohm. Dans cette section, nous proposons une courte description de ces équations ainsi que des quantités importantes pour l'étude de la turbulence et de l'effet dynamo. En particulier, nous introduisons les nombres de Reynolds et quelques théorèmes anti-dynamo importants pour comprendre les choix et limites des configurations expérimentales et numériques. Pour des raisons de simplicité, nous utilisons les conventions d'Einstein pour les sommes, c'est-à-dire pour des vecteurs f, g comprenant trois composantes 3 i=1 f i g i sera écrit f i g i .

Les équations de la MHD sont les suivantes pour un fluide incompressible, conducteur et électriquement neutre dans une approximation quasi-statique 131 (donc non-relativiste) : une dynamo provenant d'un champ de vitesse déterminé (appelé dynamo cinématique), plusieurs champs ont pu être analytiquement dérivés (dynamo de Ponomarenko [27], dynamo de Roberts [28]) et reproduits expérimentalement [1,2]. Cependant, les effets dynamo les plus intéressants proviennent du couplage non-linéaire entre l'équation du mouvement et celle d'induction. Pour ce système couplé, il n'existe pas de théorème donnant les circonstances permettant d'obtenir un champ magnétique. Pire encore, il existe de nombreux théorèmes anti-dynamo (des conditions pour lesquelles une dynamo ne peut pas apparaître peu importe la valeur de R m ) [10]. Par exemple, le théorème de Cowling [29,30] stipule qu'un champ magnétique purement axi-symétrique ne peut jamais être produit quel que soit le champ de vitesse.

∂ t u + ∂ i (u i u) -ν∆u + ∇p ρ = j × b ρ , (6.1)
∂ t b = ∇ × (u × b) -∇ × j σ , (6.2)
∇ • u = 0, (6.3
∇ × (u × µH) ∇ × ((∇ × H)/σ) ∼ σµvL = R m > 1.
Le second paramètre de contrôle intéressant pour l'expérience VKS est le nombre de Reynolds (cinétique) pour un fluide turbulent. Il caractérise le degré de turbulence de l'écoulement en comparant le terme inertiel (∂ i (u i u)) avec celui de dissipation visqueuse (ν∆u) dans l'équation (6.1):

∂ i (u i u) ν∆u ∼ vL ν = R e . (6.8)
Plus le nombre de Reynolds est grand, plus il est possible d'observer des structures petites en taille comparées à la taille du système. Ces deux nombres de Reynolds sont liés par le nombre de Prandtl magnétique qui ne dépend alors que des propriétés physiques du fluide :

P m = R m R e = µσν.
(6.9)

Comme le fluide dans l'expérience VKS est fixé (sodium liquide), la seule façon de maximiser les chances d'obtenir un champ magnétique est d'augmenter la vitesse du fluide. Pour le sodium liquide à 120°C : ν = 7 10 -7 m 2 s -1 , σ = 9.6 10 6 S m -1 et µ = µ 0 = 4π 10 -7 H m -1 . Cela donne P m = 0.8 10 -5 et impose donc R e ≳ 10 5 pour avoir une chance d'observer un champ magnétique. Cependant, dans l'expérience, il a été nécessaire de monter à R e ≳ 6 10 6 pour générer un champ magnétique. Ce grand nombre de Reynolds implique alors que quelque chose se passe au niveau des petites échelles pour l'obtention d'un champ magnétique.

. Simulations numériques

Pour étudier les équations de la magnétohydrodynamique et les mécanismes dynamo, de nombreux codes numériques existent car chaque configuration nécessite son propre cadre numérique. Ainsi, en fonction de l'objet d'étude, un système de dynamo cinématique ou les équations complètes de la MHD sont calculés. Quand l'objectif est d'étudier le seuil de l'instabilité, le champ de vitesse peut être donné (issu de résultats expérimentaux, de modèles analytiques ou de solveurs de Navier-Stokes) et seule l'équation (6.1) est résolue en faisant varier le nombre de Reynolds magnétique. Le nombre de Reynolds critique R c m est alors obtenu quand la solution b = 0 devient instable. Quand l'objectif est d'étudier les phases de croissance et de saturation du champ magnétique généré, il est nécessaire de simuler les équations complètes de la MHD (1.15-1.16).

Pour l'étude de la dynamo de Ponomarenko, un code spectral 1D pour la dynamo cinématique a été développé [31]. Dans [32], un code pseudospectral parallélisé pour les équations de la MHD dans une boite périodique est utilisé pour mieux comprendre la dynamo de Roberts.

Pour des écoulements plus complexes, des codes possédant au moins deux directions périodiques ont été développés : dans [33] pour l'étude de la dynamo provoquée par un écoulement cisaillé au sein des disques d'accrétion des galaxies et dans [34] pour la compréhension de l'action dynamo lors de la rotation de deux couches convectives. Pour la simulation de cylindres infinis (conditions périodiques en θ et z), différents codes ont été développés [35,36,37] dont le dernier a permis l'étude de la dynamo de Taylor-Couette. Il existe aussi de nombreux codes dédiés à la géométrie sphérique : codes pseudospectraux utilisant en général une décomposition poloidale-toroidale de u et b [38,39,40,41], ou des algorithmes utilisant des volumes finis [42], ou des éléments finis [43] et même des codes commerciaux comme COMSOL utilisés dans [44]. Cependant, pour simuler des domaines finis, seulement quelquesuns existent : pour des sphéroïdes [45], pour un tore [START_REF] Wu | On a dynamo driven by topographic precession[END_REF], et pour des cylindres [START_REF] Morales | Magnetohydrodynamically generated velocities in confined plasma[END_REF][START_REF] Giesecke | Generation of axisymmetric modes in cylindrical kinematic mean-field dynamos of vks type[END_REF].

Pour mieux comprendre l'expérience VKS et les mécanismes dynamo y prenant place, de nombreuses simulations numériques ont été faites. Des codes périodiques simples de MHD complets dans des géométries Cartésiennes ont permis d'estimer le nombre de Reynolds magnétique critique et d'amener à l'observation de cycles d'hystérésis liés à la force de Lorentz ainsi que les conditions pour de la dynamo non-linéaire [START_REF] Iskakov | An integro-differential formulation for magnetic induction in bounded domains: boundary element-finite volume method[END_REF][START_REF] Ponty | Dynamo action at low magnetic prandtl numbers: mean flow versus fully turbulent motions[END_REF][START_REF] Dubrulle | Bifurcations and dynamo action in a taylor-green flow[END_REF][START_REF] Ponty | Subcritical dynamo bifurcation in the taylor-green flow[END_REF]. Un code cinématique pour des cylindres infinis a permis d'optimiser le choix des turbines pour l'expérience VKS et d'obtenir les premiers résultats sur la structure du champ magnétique généré [START_REF] Krstulovic | Axial dipolar dynamo action in the taylor-green vortex[END_REF]. Des codes similaires mais pour des cylindres finis ont alors montré l'impact négatif sur la dynamo des couches de fluide derrière les turbines [START_REF] Ravelet | Bifurcations globales hydrodynamiques et magnetohydrodynamiques dans un ecoulement de von karman turbulent[END_REF][START_REF] Stefani | Ambivalent effects of added layers on steady kinematic dynamos in cylindrical geometry: application to the vks experiment[END_REF]. Cependant, la structure du champ magnétique obtenue était complètement différente de celle qui a été observée dans l'expérience un peu plus tard : le champ magnétique numérique était dominé par une composante perpendiculaire à l'axe du cylindre et une autre parallèle à cet axe (soit par un mode azimutal m = 1 compatible avec le théorème antidynamo de Cowling).

Ensuite, en introduisant les effets α et Ω dans des simulations de dy-namo cinématique, la structure moyenne du champ magnétique observée dans l'expérience a été reproduite [START_REF] Laguerre | Effects of conductivity jumps in the envelope of a kinematic dynamo flow[END_REF][START_REF] Laguerre | Impact of impellers on the axisymmetric magnetic mode in the VKS2 dynamo experiment[END_REF]17]. Cependant, elle nécessite soit des valeurs très grandes de α ou un effet α homogène qui représentent des hypothèses irréalistes. Finalement, des simulations numériques directes (DNS) et avec modèle de turbulence (LES) intégrant l'ensemble des équations de la MHD et imposant des turbines en contra-rotation ont été développées [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF].

Dans cette thèse nous effectuons de nouvelles simulations avec des paramètres déjà utilisés dans [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF]. Elles sont réalisées avec le code SFEMaNS [START_REF] Nore | Numerical simulation of the von Kármán sodium dynamo experiment[END_REF] pour les raisons suivantes : il utilise les coordonnés cylindriques permettant de représenter simplement le cylindre, les turbines sont correctement modélisées grâce à une méthode de pénalisation (voir [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] p.5), des données de simulations précédentes [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] sont accessibles et certains outils pour l'étude de la turbulence y sont déjà implémentés [START_REF] Guermond | SFEMaNS documentation[END_REF]. La simulation d'un tour complet des turbines nécessite en moyenne 10 heures de simulation sur 2040 processeurs du supercalculateur Jean-Zay fournis par l'IDRIS. La géométrie de la simulation (fig. Dans ce cadre de simulation il est alors possible d'observer la génération d'un champ magnétique très proche des résultats expérimentaux [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] en utilisant des DNS pour les faibles nombres de Reynolds ou des LES pour les grands nombres de Reynolds [START_REF] Faller | Dissipation, cascades et singularités en turbulence[END_REF] (le plus grand R e simulé étant de 10 5 et correspondant alors à un nombre de Prandtl magnétique de P c m ≃ 7 10 -4 au seuil dynamo). Dans cette thèse, nous utilisons les résultats obtenus pour R e = 1500 dans le cadre d'une DNS pour pouvoir étudier les transferts d'énergie jusqu'aux plus petites échelles sans modélisation. Ce nombre de Reynolds est assez faible comparé au cas expérimental mais donne tout de même une idée concrète des ingrédients nécessaires à l'apparition de l'effet dynamo. En particulier, il a été observé que l'utilisation de turbines avec une haute perméabilité magnétique diminue le nombre de Reynolds magnétique critique pour lequel un champ magnétique est produit. Ainsi, comme observé dans l'expérience, le seuil pour lequel un champ magnétique est généré dépend du matériau utilisé pour les turbines. Cependant, l'origine des mécanismes dynamo reste un mystère.

. Plan de la thèse

Pour mieux comprendre ces résultats numériques et donc d'où provient l'effet dynamo dans l'expérience VKS, nous implémentons les bilans d'énergie locaux dérivés des équations de la MHD dans le code SFEMaNS.

Avant de proposer un résumé de chacun des chapitres de cette thèse, nous décrivons sa structure.

Dans le premier chapitre-2 résumé dans la section 6.2, nous décrivons le programme utilisé ainsi que son contexte numérique. Nous présentons aussi quelques améliorations faites pour la répartition de la grille de simulations sur les processeurs ainsi qu'une étude sur une méthode potentiellement plus rapide pour simuler les équations.

Dans le second chapitre-3 résumé dans la section 6.3, nous décrivons les différents outils de post-traitement permettant de calculer, analyser et visualiser les transferts d'énergie locaux. Ces outils sont essentiels pour que l'étude des ces termes soit automatisée et efficace.

Finalement, dans le chapitre-4 résumé dans la section 6.4, nous dérivons les équations pour les transferts d'énergie ainsi que les simulations utilisées. Nous analyserons alors les différents transferts pour deux cas différents : un pour des turbines en acier et l'autre pour des turbines en fer doux.

L'objectif de cette thèse est ainsi de découvrir quels sont les mécanismes dynamo responsables de la génération du champ magnétique dans les fluides conducteurs turbulents. Le volume de simulation que nous utilisons est montré en Fig. 1.1 où la rotation des turbines est implémentée grâce à une méthode de pénalisation (voir [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] p.5). La plupart des fonctions pour des calculs simples et de parallélisation sont déjà implémentées. Cependant, il reste à construire les outils pour calculer les nouveaux termes de transferts d'énergie.

. Améliorations du code SFEMaNS

. Approximations numériques dans SFEMaNS

Le code utilise une décomposition hybride pour les champs. L'utilisation des coordonnées cylindriques permet une décomposition en M modes de Fourier suivant θ permettant alors de réduire le problème 3D à des résolutions de problèmes 2D dans le plan (r, z). Ces problèmes sont résolus en utilisant le cadre des éléments finis et de la formulation faible des équations.

Les méthodes des éléments finis permettent d'approximer l'espace des fonctions par un sous-ensemble de fonctions appelées fonctions tests ϕ i . Ces fonctions sont définies sur une grille (voir Fig.

∂ t T + ∇ • (uT ) -K∆T = f,
A m,c ij = r,z r 3 2τ ψ j ψ i + K(∂ r ψ j ∂ r ψ i + ∂ z ψ j ∂ z ψ i + m 2 r 2 ψ j ψ i ) , (6.12
)

B m,c i = np j=1 r,z rF m,c j ψ j ψ i , (6.13)
où

F = 4T n -T n-1 2τ + f n+1 -∇ • ((2u n -u n-1 )(2T n -T n-1 )) .
Ce système linéaire est alors découplé en modes de Fourier permettant une parallélisation efficace. De plus, comme il s'agit d'une simple inversion de matrice, des outils efficaces d'algèbre linéaire peuvent être utilisés. Pour évaluer les intégrales dans A et B, une quadrature de Gauss est utilisée qui permet d'approximer l'intégrale en utilisant une somme pondérée de points bien choisis.

. Amélioration de la distribution du maillage de simulation

Une des principales limites des grilles de simulation dans SFEMaNS est leur utilisation de la mémoire. Le problème provient du fait que chaque processeur garde en mémoire l'intégralité de la grille de simulation et pas seulement la partie lui étant attribuée. Comme chaque processeur ne posséde que de 4GB à 8GB de mémoire RAM, il est impossible de simuler des grilles 2D avec plus de 500.000 points. Ceci limite sévèrement le nombre de Reynolds cinétique maximal pouvant être simulé et est très inefficace pour la mémoire rendant difficile le calcul de nouvelles quantités. Cette limite est par ailleurs d'autant plus facile à atteindre pour les équations de la MHD car elles nécessitent trois grilles (pour le champ de vitesses, de pression et magnétique).

Pour résoudre ce problème, nous avons implémenté une nouvelle méthode pour distribuer la grille de telle manière que chaque processeur ne conserve en mémoire que les cellules et points lui étant strictement nécessaires pour construire l'algèbre linéaire. Cette nouvelle implémentation permet d'obtenir des simulations contenant 300.000 points de simulation par processeur pour des problèmes de Navier-Stokes 2D. Ainsi, en utilisant de nombreux processeurs, il nous est possible de faire des simulations avec des centaines de millions de points sur nos maillages 2D.

Deux principales difficultés ont été rencontrées. Premièrement, il est nécessaire de respecter les conventions (figure 6.2) sur la numérotations des points, voisins et bords des cellules pour que les algèbres linéaires soient construites correctement. Deuxièmement, à cause de la parallélisation, certains processeurs ont besoin de l'information de cellules ne leur étant pas attribuées.

Par exemple, dans la figure 6.3 , le processeur 0 (bleu foncé) a besoin de l'information des cellules en bleu clair, rouge et jaune contenant un des gros points bleus. Ce problème est résolu en sauvegardant pour chaque processeur ces cellules voisines de manière à pouvoir reconstruire les interactions entre ces cellules. De plus, pour des raisons techniques, il est nécessaire de construire trois numérotations pour chacun des éléments de la grille (points, segments et cellules) : une numérotation globale, une au niveau du processeur et une au niveau de la cellule.

. Étude de méthodes avec préconditionnement pour résoudre Navier-Stokes

Cette section est un bref résumé de la section 2.3 composée d'un article soumis à Computer Methods in Applied Mechanics and Engineering (Octobre 2023).

Le but est d'implémenter et tester des méthodes de préconditionnements pour résoudre le problème de Stokes dans le cadre de SFEMaNS. L'objectif final est de voir si ces méthodes sont plus rapides que celle de prédictioncorrection pour l'instant utilisée dans SFEMaNS. Cela pourrait économiser beaucoup de temps de calcul mais la conclusion est que, pour l'instant, cela n'est pas le cas.

Pour des raisons de simplicité, l'aspect mathématique est grossièrement vulgarisé.

Contexte

Résoudre le problème de Stokes généralisé a une longue histoire commençant par Fortin et Glowinski [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF] et par Cahouet et Chabart [START_REF] Fortin | Augmented Lagrangian methods[END_REF]. Dans [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF], les auteurs ont introduit la méthode du Lagrangien augmenté permettant d'améliorer le conditionnement du complément de Schur en pression. Dans [START_REF] Fortin | Augmented Lagrangian methods[END_REF], les auteurs ont proposé un préconditionnement du problème de Stokes se comportant bien dans les deux limites µτ h -2 → 0 et µτ h -2 → ∞ avec µ la viscosité, τ le pas de temps et h la taille caractéristique des éléments du maillage de simulation. Ces méthodes sont par la suite combinées [START_REF] Cahouet | Some fast 3D finite element solvers for the generalized Stokes problem[END_REF][START_REF] Benzi | An augmented Lagrangian-based approach to the Oseen problem[END_REF][START_REF] Benzi | Modified augmented lagrangian preconditioners for the incompressible navier-stokes equations[END_REF] donnant alors un bon préconditionnement pour le problème de Schur généralisé. Ces différentes approches sont testées numériquement dans [START_REF] Benzi | Analysis of augmented lagrangian-based preconditioners for the steady incompressible navier-stokes equations[END_REF][START_REF] Farrell | An augmented lagrangian preconditioner for the 3d stationary incompressible navier-stokes equations at high reynolds number[END_REF][START_REF] Moulin | Augmented lagrangian preconditioner for large-scale hydrodynamic stability analysis[END_REF] et montrent en effet que le nombre de MINRES ou GMRES itérations pour obtenir une solution se comporte correctement en fonction de la viscosité et de la taille caractéristique de la grille de simulation.

Le problème de Stokes correspond à résoudre pour la pression p et la vitesse u : avec B correspondant à la matrice en éléments finis définissant le gradient ∇ et A := 1 τ M V + µE V où M V correspond à la matrice identité du maillage en vitesse en éléments finis (appelée matrice de masse) et E V la matrice associée à l'opérateur 1 4 (∇ + ∇ T ) : (∇ + ∇ T ). Ainsi il est possible de résoudre ce système par inversion de matrices.

∂ t u - µ 4 [(∇ + ∇ T ) : (∇ + ∇ T )]u + ∇p = f, ∇ • u = 0 (
Le complément de Schur de ce problème est alors défini par S := BA -1 B T et rend le système (6.15) équivalent à : Pour diminuer ce nombre de conditionnement, la méthode du Lagrangien augmenté a été développée dans [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF] et transforme le système (6.15) en un système équivalent (donc possédant les mêmes solutions) :

SP = -BA -1 F, AU = F + B T P. ( 6
(A + λµB T M -1 Q B)U -B T P = F, BU = 0, (6.17)
où λ est un paramètre adimensionné positif et M Q la matrice de masse pour la pression. Il est de nouveau possible de construire un complément de Schur pour ce système et ainsi de simplement résoudre :

S λµ := B(A + λµB T M -1 Q B) -1 B T , (6.18a
)

S λµ P = -B(A + λµB T M -1 Q B) -1 F, AU = F + B T P. ( 6
T M -1 Q B) -1 .

Préconditionneur

Une autre méthode pour diminuer le nombre d'itérations nécessaires pour inverser le problème Schur est de tout d'abord le multiplier par une matrice permettant de le rapprocher de l'identité. Pour le problème (6.16), dans [START_REF] Fortin | Augmented Lagrangian methods[END_REF], les auteurs proposent de multiplier l'équation en pression par une de ces deux matrices :

CaCh 1 := τ -1 (B M -1 V B T ) -1 + µM -1 Q , (6.19a) CaCh 2 := τ -1 (L Q ) -1 + µM -1 Q . (6.19b)
De notre côté, pour le problème avec le Lagrangien augmenté (6.18b), nous proposons :

C λµ := τ -1 (B(M V ) -1 B T ) -1 + µ(1 + λ)M -1 Q . (6.20) Dans la limite µτ /h 2 → 0, on obtient S λµ → τ B(M V ) -1 B T et donc S λµ de- vient τ -1 (B(M V ) -1 B T ) -1 . Ainsi, en multipliant C λµ S λµ ∼ I d dans cette lim- ite.
De même, dans la limite µτ /h 2 → ∞, on obtient

S λµ → µ -1 B(E V + λB T M -1 Q B) -1 B T . Comme E V et B T M -1 Q B possèdent des valeurs similaires, on peut approximer S λµ dans la limite µτ /h 2 → ∞ par µ -1 (1+λ) -1 B(B T M -1 Q B) -1 B T . Alors l'inverse de S λµ se comporte approximativement comme µ(1 + λ)M -1 Q ce qui correspond à l'inverse de C λµ dans cette même limite.
Ainsi, pour ces deux limites, nous nous attendons à un nombre très faible d'itérations nécessaires pour des méthodes de résolution itérative comme GMRES. Des arguments similaires marchent aussi pour les préconditionneurs 6.19a et 6.19b pour le problème sans Lagrangien augmenté 6.16.

Cadre numérique

Pour tester ces méthodes, nous utilisons des éléments finis continus de Lagrange en 2D sur des grilles non-structurées qui sont des raffinements successifs de la figure 6.4. Les éléments finis utilisés dans ce résumé sont les éléments de Taylor-Hood P 2 /P 1 . Pour chaque raffinement, le nombre de processeurs est multiplié par 4 pour conserver un nombre de points par processeur équivalent. Le maillage testé le plus fin posséde environ 120 millions de points pour le maillage en vitesse. Le tout est codé dans SFEMaNS et utilise donc PETSC pour la parallélisation en domaines sur la grille 2D.

En fonction des situations, nous utilisons soit un solveur linéaire direct MUMPS [START_REF] Golub | On solving block-structured indefinite linear systems[END_REF] ou algébrique multigrilles BoomerAMG [START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF] pour inverser les matrices dont il est possible de construire les coefficients. Pour des raisons de reproductibilité voici les paramètres utilisés pour BoomerAMG :

-pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_strong_threshold 0.7 or 0.1 -pc_hypre_boomeramg_coarsen_type Falgout -pc_hypre_boomeramg_relax_type_all Chebyshev

τ -1 M V + µL V + µ(1 + λ)D = A + λµD, (6.21)
τ -1 M V + µL V + µD = A, , (6.22)
τ -1 M V + µL V ., (6.23)
avec D les termes non diagonaux de la matrice Laplacien L V . Ces trois expressions semblent a priori de bonnes approximations de

A + λµB T M -1 Q B.
Pour comparer leur efficacité, nous mesurons combien de temps par degré de liberté est nécessaire pour résoudre le problème. Les résultats pour ces différents préconditionneurs sont reportés dans la figure 6.7. Les deux types de solveur MUMPS et BommerAMG sont testés pour inverser le préconditionneur (6.23). Clairement seul le préconditionneur (6.23) permet d'obtenir une bonne parallélisation pour des maillages de plus en plus fins (le temps par degré de liberté restant constant). Aussi la méthode directe MUMPS ne marche plus au-delà de 2 millions de points par manque de mémoire vive. A+λµD, BoomerAMG A, BoomerAMG

τ -1 M V +µL V , BoomerAMG τ -1 M V +µL V , MUMPS

Résultats et conclusion pour le problème de Stokes

Nous étudions maintenant les performances des préconditionneurs pour le problème de Stokes dans sa totalité. Pour les préconditionneurs CaCh 1 (6.19a) et CaCh 2 (6.19b), le problème sans Lagrangien augmenté (6.16) est résolu.

Et, pour le préconditionneur C λµ (6.20), le problème avec Lagrangien augmenté (6.18b) est résolu en utilisant le préconditionneur (6.23) pour le problème en vitesse.

Dans la figure 6.8, nous montrons le résidu de la méthode GMRES (c'està-dire la précision obtenue sur la résolution du problème) en fonction de l'itération GMRES atteinte. Pour atteindre la précision relative demandée (de 10 -10 ) lors de la résolution itérative du problème, nous observons bien une diminution du nombre d'itérations GMRES pour le préconditionneur C λµ (6.23). Il est par ailleurs plus faible en prenant λ = 10 plutôt que λ = 1. Par contre, le nombre d'itérations pour CaCh 1 (6.19a) diverge sur les maillages fins. Cependant, en regardant l'efficacité en temps de calcul (fig. 6.9), les méthodes avec le moins d'itérations se retrouvent avec des performances similaires de l'ordre de la milliseconde par degré de liberté. Ce même ordre de grandeur est obtenu dans la littérature [START_REF] Benzi | Modified augmented lagrangian preconditioners for the incompressible navier-stokes equations[END_REF][START_REF] Benzi | Analysis of augmented lagrangian-based preconditioners for the steady incompressible navier-stokes equations[END_REF][START_REF] Farrell | An augmented lagrangian preconditioner for the 3d stationary incompressible navier-stokes equations at high reynolds number[END_REF]. Les performances en parallèle sont globalement bonnes (peu de variation du temps de calcul par degré de liberté en raffinant le maillage).

τ -1 (B M -1 V B T ) -1 + µ(1 + λ)M -1 Q avec µ = 1 (bleu), µ = 10 -2 (
Pour conclure, l'efficacité en temps pour résoudre les problèmes de Stokes et de Navier-Stokes en utilisant les techniques de complément de Schur n'est pas, pour l'instant, meilleure que des méthodes de prédiction-correction (en vitesse ou pression). Ces méthodes de prédiction-correction [START_REF] Henson | BoomerAMG: a parallel algebraic multigrid solver and preconditioner[END_REF][START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires ii[END_REF] ayant une efficacité d'environ 10 -2 ms par degré de liberté, il semble difficile de gagner deux ordres de grandeur même en optimisant les paramètres des préconditionneurs.

. Outils de post-traitement

Le code SFEMaNS permet de base de simuler les équations de la MHD et donne donc l'accès aux champs de vitesse u, de pression p et au champ magnétique. Ainsi, il est nécessaire d'implémenter des outils pour calculer, étudier et visualiser les transferts d'énergie. Dans un premier temps, il est nécessaire de calculer ces quantités (voir équations 6.31 et 6.32) dans SFEMaNS. Ensuite,

(B M -1 V B T ) -1 + µ(1 + λ)M -1
Q par différents préconditionneurs. Axe des Y : temps de simulation×nb. processeurs nb. vel.+press. degrés de liberté en millisecondes ; axe des X : nombre de points pour le maillage de vitesse. µ = 1 ( ); µ = 10 -2 ( ), µ = 10 -4 ( ).

il faut transformer les données en sortie de SFEMaNS dans des formats de fichiers permettant l'analyse statistique ou la visualisation. Finalement, au vu des grandes quantités de données, il est important d'automatiser la visualisation et l'analyse des champs obtenus. Cette automatisation est faite grâce au couplage entre Python et Paraview. Pour facilement transférer ces outils d'une plateforme à l'autre et les utiliser, ils sont tous disponibles pour notre groupe de recherche dans un GitLab privé en tant que projets différents. Un fichier README.md est inclus dans chacun de ces projets contenant un tutoriel simple. L'objectif de cette section est de résumer les principaux paramètres de ces outils développés en plus de détails dans le chapitre 3.

. Méthodes pour le calcul des termes de transferts d'énergie

Dans un but de conserver une haute précision et une grande efficacité, les différents termes de transferts énergétiques sont calculés en utilisant le code SFEMaNS. Le projet correspondant se nomme 'DR_ MHD' dans GitLab. Toutes les options utiles pour faire fonctionner le code sont dans le fichier 'data_example' et il suffit de le copier-coller dans le fichier 'data' utilisé pour les simulations temporelles.

Pour pouvoir analyser ce qu'il se passe à différentes échelles, il est nécessaire de les séparer. Ceci est en particulier important pour trouver les évènements rares et intenses de la turbulence pouvant être des indices de quasisingularités ou des processus de dynamo rapide. Pour faire cette séparation, nous utilisons un filtrage se présentant sous la forme d'une résolution d'équation différentielle. Ce choix de filtre permet d'avoir une méthode qui est rapide et efficace sur les grilles non-structurées comme c'est le cas pour SFEMaNS.

Dans le cadre de SFEMaNS, résoudre des équations différentielles linéaires est facile. Ainsi, nous définissons un champ filtré g = u ℓ comme la solution de : gℓ 2 ∆g = u, (6.24) avec u un vecteur ou scalaire et ℓ ∈ R + la taille du filtre.

Le filtrage est alors transformé en un problème d'algèbre linéaire similaire à ce qui est fait dans la section 6.2.1. L'équation (6.24) devient alors, pour un mode de Fourier m et la partie cosinus d'un scalaire ou d'un vecteur :

A m,c ij g m,c j = u m,c i , (6.25)
avec

A m,c ij = r,z r(ψ j ψ i + ℓ 2 (∂ r ψ j ∂ r ψ i + ∂ z ψ j ∂ z ψ i + m 2 r 2 ψ j ψ i )). (6.26)
Comme la construction de la matrice A prend beaucoup de temps, une matrice est sauvegardée pour chaque échelle ℓ.

Une fois le filtre construit, il est alors possible de calculer chacun des termes de transferts d'énergie. Comme il s'agit simplement de successions de produits scalaires et vectoriels ainsi que de dérivées, la seule difficulté provient du fait que les champs soient définis sur des maillages de simulation différents. Cela nécessite de parfois projeter des champs existants sur le maillage du champ magnétique vers le maillage du champ de vitesse. Comme le coût en mémoire vive et en temps de calcul est important, une attention particulière a été prêtée à l'optimisation de cette partie du code. Ainsi, par exemple, les quantités indépendantes de l'échelle du filtre ℓ (comme le courant électrique j) ne sont calculées qu'une seule fois et gardées en mémoire vive pour d'autres échelles. Du côté optimisation mémoire, les tableaux de données utilisés pour stocker des champs scalaires et vectoriels dans le code sont réutilisés plusieurs fois pour éviter plusieurs allocations et désallocations de mémoire. Ceci est fait tout en minimisant le nombre de tableaux utilisés. Ainsi un fichier nommé 'memory_usage.ods' (traduit par utilisation de la mémoire) a été créé permettant de suivre quels termes sont présents dans quels tableaux à chaque instant du code.

En plus du calcul de tous les termes de transferts, d'autres quantités utiles sont calculées comme les différentes hélicités et courants. Finalement, les moyennes spatiales sur l'ensemble du volume de simulation de toutes ces quantités sont sauvegardées pour chaque pas de temps temporel dans un fichier nommé 'fort.75'.

De nombreux tests ont permis de vérifier que les termes sont calculés correctement. Par exemple, pour ℓ = 0, certains termes de transfert doivent être nuls. Cependant cela n'est pas le cas pour tous à cause du repliement de spectre de Fourier (aussi connu sous le nom d'aliasing) lors des produits entre champs. Pour éviter cela, SFEMaNS a déjà implémenté les stratégies d'anti-aliasing pour de simples produits à deux termes f × g mais cela n'est pas suffisant pour des produits à trois termes ou plus f × (g × h). Ainsi, il y a une erreur de quelques pourcents sur le calcul des énergies de transfert contenant des produits multiples.

. Transferts des données SFEMaNS vers des formats de fichiers utilisables

Les sorties SFEMaNS se trouvent être dans un format spécifique contenant de nombreuses informations superflues à l'étude des termes mais nécessaires pour relancer les simulations temporelles. De plus, les différents modes de Fourier sont sauvegardés de manière complexe en fonction de la structure du maillage et de sa répartition sur les processeurs. Pour pouvoir utiliser et visualiser ces données de manière efficace, un code simple permettant de changer le format de ces fichiers a été développé. Cela est fait nécessairement dans le code SFEMaNS et a pour nom de projet sur GitLab 'suites_to_bins'.

Ce code prend toutes les itérations temporelles des différents termes déjà calculés et donne en sortie soit des fichiers binaires pour le traitement statistique via Python ou des fichier pvd/vtu pour la visualisation dans Paraview. Encore un fois, toutes les options sont incluses dans un fichier 'data_example' qui ont simplement besoin d'être copiées dans le fichier 'data' de la simulation temporelle. Les différentes options permettent de sélectionner les champs à transférer mais aussi de faire des moyennes temporelles de ceux-ci. Il est aussi possible de les normaliser par la valeur moyenne de l'énergie magnétique au cours du temps pour pouvoir ainsi mieux mettre en évidence le mode propre de la croissance exponentielle du champ magnétique dans la phase de croissance de la dynamo.

. Automatisation de la visualisation avec Paraview

La visualisation est un point important pour l'analyse des phénomènes physiques car elle permet de mettre en évidence où et quand ils arrivent. De nombreuses options permettent de choisir différentes échelles de couleurs et de stipuler les iso-surfaces à représenter. Par défaut, elles sont en échelles bi-logarithmiques car les quantités étudiées varient grandement en intensité à cause des turbines.

Ces visualisations sont codées de manière à éviter de charger en mémoire plusieurs fois un même champ pour des visualisations différentes. Ainsi, pour parcourir tous les champs à visualiser, d'autres fonctions sont utilisées pour créer les tableaux Paraview à donner aux fonctions de visualisation. Cela permet par ailleurs de faire des calculs supplémentaires sur les champs avant de les visualiser grâce aux modules de calculs dans Paraview. Pour rendre les visualisations plus rapides, une parallélisation de celles-ci a été mise en place.

Comme la quantité d'images produites est importante, il est nécessaire d'avoir un outil permettant de les visualiser directement sur un super-calculateur sans les transférer en local. Par exemple, dans notre cas avec environ 20 champs scalaires sur 120 itérations temporelles et trois échelles de filtrage, cela donne des dizaines de milliers d'images ( 200 Go de mémoire) et ceci seulement pour les isosurfaces 3D. Ainsi un GUI simple a été codé en Python dans le projet GitLab 'image_viewer_sfemans'. Il suffit alors simplement de donner le répertoire contenant les images pour y avoir accès (voir figure 6.11). Par exemple, quelques outils d'analyse statistique en Python qui ont été implémentés en début de thèse ont ainsi pu être restructurés par Nikita Allaglo (stagiaire de L3, printemps 2023) et Rémi Bousquet (stagiaire de M2, printemps 2023). Malheureusement, ils n'ont pas pu être utilisés dans le cadre de cette thèse.

. Transferts énergétiques dans l'expérience de von Kármán Sodium

Cette section propose un résumé du chapitre 4 qui se présente sous la forme d'un article soumis à Physics of Plasmas (Août 2023). Elle contient les démonstrations analytiques des différents transferts d'énergie entre champs et entre échelles pour les équations de la magnétohydrodynamique avec variations de perméabilité magnétique et de conductivité électrique. Ces transferts sont ensuite appliqués à l'étude de deux simulations numériques ayant généré un champ magnétique mais avec des matériaux différents pour les turbines (acier et fer doux).

. Bilans énergétiques locaux

Pour évaluer les bilans d'énergie locaux permettant de décrire les transferts d'énergie entre champ magnétique et champ de vitesse, nous utilisons une méthode de filtrage. Cette stratégie s'inspire de travaux mathématiques pour dériver des dissipations anormales dues aux singularités dans le contexte des équations de Navier-Stokes et dans le contexte des équations MHD. Notre travail étend cette méthode en introduisant des variations possibles des propriétés des matériaux (µ, σ). Ceci est important pour bien prendre en compte ce qui se passe dans nos simulations ainsi que pour évaluer l'impact des sauts de conductivité électrique et de perméabilité magnétique.

Ainsi, en appliquant la méthode de filtrage (décrite en section 6.24) aux équations de la MHD (6.1) et (6.2), nous obtenons : En sommant le produit scalaire de ρu/2 avec l'équation (6.27) et le produit scalaire de ρu ℓ /2 avec l'équation (6.1) et en faisant des produits similaires pour b, les équations pour l'énergie cinétique filtrée E c = ρuu ℓ 2 et magnétique filtrée Elle est appelée dissipation anormale car, dans la limite ℓ → 0, elle peut ne pas disparaître [23]. De même, D b est un transfert d'énergie entre échelles mais cette fois dû aux interactions entre le champ de vitesse et le champ magnétique. Comme dans le cas cinétique, ce terme peut produire une contribution anormale si les champs sont suffisamment irréguliers. Si la contribution est positive, il s'agit alors d'une dissipation anormale. Cependant, il est aussi possible que sa contribution soit négative résultant alors en une amplification du champ magnétique. Cela correspondrait à une situation pour laquelle l'effet dynamo est produit par un champ de vitesse irrégulier [7], un mécanisme aussi connu sous le nom de dynamo rapide [8].

∂ t u ℓ + ∂ i u i u ℓ -ν∆u ℓ + ∇p ℓ ρ = - b × j ℓ ρ + f ℓ , (6.27)
∂ t b ℓ = ∇ × (u × b ℓ ) -∇ × ( j σ ) ℓ , (6.28)
∇ • u ℓ = 0, (6.29)
E m = bb ℓ
J M = 1 2µ { j σ × b ℓ + ( j σ ) ℓ × b + b × (u × b ℓ ) + b ℓ × (u × b)}, (6.34)
P = ρ 2 (u • f ℓ + u ℓ • f ), (6.35)
D ν = ρν∂ i u j ℓ ∂ i u j , (6.36)
D σ = j • j ℓ σ , D σ = ( j σ ) ℓ - (j ℓ ) σ • j 2 , (6.37)
D u = ρ 2 (uu i ∂ i u ℓ -uu i ∂ i u ℓ ), (6.38)
D b = 1 4 (u ℓ • (b × j) + u • (b × j ℓ ) -(u × b ℓ ) • j -(u × b) • j ℓ ), (6.39)
T u→b = 1 4 (u ℓ • (b × j) + u • (b × j ℓ ) + (u × b ℓ ) • j + (u × b) • j ℓ ), (6.40)
M µ = - b • b ℓ 2µ 2 ∂ t µ, (6.41)
D µσ = j 2σ • ∇ × (b ℓ ) µ -H ℓ , D µT = - (u × b) 2 • ∇ × (b ℓ ) µ -H ℓ . ( 6
Du côté droit, le terme (-T u→b ) correspond à la puissance de la force de Lorentz et reflète le transfert direct d'énergie entre les champs u et b. Il apparaît naturellement avec un signe opposé du côté droit du bilan d'énergie magnétique (6.32). Il est ainsi le candidat principal pour l'apparition et le maintien du champ magnétique et est donc un terme source pour l'effet dynamo.

Le dernier terme, P, est simplement l'injection d'énergie mécanique par le forçage (dans notre cas les turbines).

En plus des termes T u→b et D b , de nouveaux termes intéressants apparaissent dans l'équation (6.32). Du côté droit, J M est simplement le courant d'énergie décrivant comment l'énergie magnétique est transportée par le fluide. D σ est la dissipation Joule et donc une perte d'énergie pour le champ magnétique. D σ est la contribution aux transferts d'énergie entre échelles due aux sauts de conductivité électrique. De même, D µσ et D µT sont des termes de transferts d'énergie entre échelles dus aux irrégularités de la perméabilité magnétique. Finalement, du côté gauche, un autre terme source pour l'effet dynamo apparaît, M µ . Il est produit par des variations temporelles de perméabilité magnétique et est un transfert direct des turbines vers le champ magnétique.

. Simulations numériques

Nous étudions alors ces transferts d'énergie pour la géométrie de von Kármán Sodium représentée en figure 6.1. Nos simulations ont été faites en utilisant le code SFEMaNS présenté en section 6.2 en utilisant le supercalculateur Jean-Zay de l'IDRIS. Les simulations utilisées correspondent à deux séries décrites dans [START_REF] Giesecke | Electromagnetic induction in non-uniform domains[END_REF] (Table 3). Nous utilisons des unités adimensionnées de telle sorte que la longueur caractéristique du système L corresponde au rayon du cylindre interne R cyl . En dénotant σ 0 la conductivité électrique du sodium liquide, ρ sa densité et µ 0 la perméabilité magnétique du vide, le champ magnétique est adimensionné en utilisant l'échelle de Alfvén B = U √ ρµ 0 , avec U = ωR cyl où ω est la vitesse angulaire des turbines. Les deux paramètres de contrôle sont alors R m = µ 0 σ 0 R 2 cyl ω, le nombre de Reynolds magnétique et R e = R 2 cyl ω/ν, le nombre de Reynolds cinétique. Nous définissons aussi la conductivité électrique relative σ r et la perméabilité magnétique relative µ r . Ces quantités peuvent ne pas être uniformes pour des murs et turbines de matériaux différents. Dans notre cas, σ r = 1, µ r = 1 dans la région {(r, θ, z) ∈ Nous allons maintenant succinctement résumer les principaux résultats observés pour les deux phases (en croissance et en saturation). Pour cela, nous nous concentrerons seulement sur le bilan d'énergie pour le champ magnétique et ne discuterons que l'effet dynamo et non les diverses contributions anormales.

. Mécanismes dynamo

En phase de croissance, nous obtenons les bilans énergétiques pour le champ magnétique montrés en figure 6.12. Pour mieux souligner l'importance respective des différents termes, ils sont renormalisés par l'énergie magnétique moyenne ⟨E m ⟩(t). Cela permet d'éviter de seulement voir la croissance exponentielle de tous ces termes due à la croissance exponentielle du champ magnétique. Pour résumer, la puissance de la force de Lorentz T u→b est nécessaire au démarrage de la dynamo dans les deux cas, soulignant le rôle crucial du champ de vitesse pour initialiser les mécanismes dynamo. Cependant, le développement de la dynamo et son maintien diffèrent fortement pour les deux cas. Pour les turbines en fer doux, la rotation des turbines via M µ devient le seul acteur dynamo pendant que l'action de la force de Lorentz devient anti-dynamo. Pour les turbines en acier, le phénomène dynamo provient uniquement du champ de vitesse grâce à plusieurs mécanismes se passant au niveau des turbines et de la couche de mélange.

. Conclusion

Dans cette thèse, l'objectif principal a été de comprendre l'effet dynamo dans les fluides conducteurs et turbulents. Pour cela, la configuration de von Kármán Sodium a été simulée en utilisant le code SFEMaNS. Le choix de cette configuration est motivé par la possibilité de comparer les résultats numériques à des expériences faites au CEA de Cadarache. Pour mieux comprendre les mécanismes dynamo, nous avons décidé d'étudier les transferts d'énergie entre les champs de vitesse et magnétique ainsi qu'entre les différentes échelles de l'écoulement. Cette approche est une extension directe à la magnétohydrodynamique de ce qui a été fait pour les équations de Navier-Stokes dans le cadre de la thèse de H.Faller [START_REF] Guermond | SFEMaNS documentation[END_REF]. Ces choix ont mené à l'implémentation de nombreux outils de post-traitement ainsi qu'à quelques modifications dans le code SFEMaNS.

Une grande partie des améliorations que nous avons apportées à SFE-MaNS vient de la nouvelle distribution du maillage de simulation sur les processeurs. Cela a permis d'économiser de la place mémoire et d'atteindre 300.000 points de simulations par processeur pour des simulations 2D de Navier-Stokes. Ainsi, il est maintenant possible de lancer des simulations numériques pour de plus grands nombres de Reynolds cinétique à la fois pour l'hydrodynamique et la MHD. En plus de cela, nous avons testé une méthode de préconditionnement pour résoudre le problème de Stokes dans SFE-MaNS pour une éventuelle utilisation dans des méthodes de simulations temporelles implicites et/ou explicites. Les différentes techniques de préconditionnement testées obtiennent une performance similaire d'environ 10 -3 secondes de temps de calcul dépensé par degré de liberté avec une très bonne efficacité de la parallélisation. Cependant, elles ne sont pas suffisamment performantes comparées à la méthode de prédiction-correction utilisée actuellement dans SFEMaNS (ayant un temps de calcul par degré de liberté d'environ 10 -5 secondes).

De plus, pour calculer et analyser de manière efficace les différents transferts d'énergie, nous avons mis en place de nombreux outils de post-traitement. Une attention particulière a été portée pour économiser la mémoire et le temps de calcul durant l'évaluation de ces termes. Un outil général permettant de transférer les sorties de SFEMaNS en des fichiers simples d'utilisation pour un traitement statistique ou de la visualisation a été développé. Une automatisation de Paraview pour la visualisation des structures des champs avec de nombreuses fonctionnalités a été développée. Tous ces outils ont été nécessaires pour rendre l'analyse et la visualisation des transferts d'énergie efficaces au vu de la quantité de données à analyser. Ces outils sont par ailleurs indépendants les uns des autres et peuvent être donc utilisés pour toute simulation future utilisant SFEMaNS.

Enfin, la dérivation des bilans d'énergie pour la magnétohydrodynamique a amené à la découverte de nombreux termes de transferts d'énergie. En particulier, l'influence des variations temporelles et spatiales de la perméabilité magnétique sur les transferts d'énergie entre échelles est maintenant correctement décrite. Ainsi, dans nos deux cas d'étude, nous observons des mécanismes dynamo très dépendants de la perméabilité magnétique des turbines. Dans le cas des turbines en acier, la puissance de la force de Lorentz permet de compenser tous les termes de dissipation alors que, pour des turbines en fer doux, elle devient anti-dynamo. Pour des turbines en fer, dans le régime saturé, la seule contribution à l'effet dynamo provient alors de la rotation des turbines ayant une grande perméabilité magnétique. Dans le futur, une analyse plus complète des différents transferts d'énergie serait utile. En particulier, des analyses à plus grands nombres de Reynolds cinétique des transferts d'énergie responsables des dissipations anormales pourraient mener à des indices sur l'apparition de quasi-singularités et de phénomènes de dynamo rapide dans la MHD. De plus, des simulations pour des perméabilités magnétiques intermédiaires devraient apporter une meilleure compréhension de la bascule se passant entre la puissance de la force de Lorentz et l'effet de la rotation de la perméabilité magnétique. Dans toutes ces simulations, il serait par ailleurs intéressant de suivre les évènements rares et extrêmes et les autres structures cohérentes grâce aux outils de visualisation. Finalement, un traitement statistique de ces données serait nécessaire pour mieux souligner les corrélations entre les différents phénomènes et ainsi construire de nouveaux modèles permettant de mieux comprendre l'effet dynamo.

Titre:

  Comment produire du champ magnétique dans les fluides, les étoiles et les planètes ? Mots clés: Magnétohydrodynamique, analyse numérique, mécanique des fluides

Figure 1 . 1 :

 11 Figure 1.1: Representation of the different simulated volumes: fluid in movement (green), fluid kept at rest (yellow), other conductive solid walls (red), and turbines in movement (dark red).

Figure 2 . 1 :

 21 Figure 2.1: Refinement of a mesh. No refinement (left panel), one refinement (middle panel), two refinements (right panel).

3 v e r t e x 1 v e r t e x 2 Figure 2 . 2 :

 31222 Figure 2.2: Nodes, edges and neighbours convention for a cell.

Figure 2 . 3 :

 23 Figure 2.3: Numbering of points and edge orientations in one cell for P 1 (left panel), P 2 (middle panel) and P 3 (right panel) finite elements.

  that mesh: • mesh_glob%jj(:,m) gives the points ∈ [[1, N ]] of the cell m (there are 3 points per cell because we use P 1 finite elements), • mesh_glob%jce(:,m) gives the edges ∈ [[1, M e ]] of the cell m, • mesh_glob%neigh(:,m) gives the neighbours ∈ [[1, M ]] of the cell m (or 0 if the neighbour is on the external boundaries).

2 . 1

 21 are respected. This means that, for Fig. 2.6 (left), v 1 =mesh_glob%jj(1,m), v 2 =mesh_glob%jj(1,m) and v 3 =mesh_glob%jj(3,m). Similarly, e 1 =mesh_glob%jce(1,m), e 2 =mesh_glob%jce(1,m) and e 3 =mesh_glob%jce(3,m).

Figure 2 . 4 :o c e s s o r 4 Figure 2 . 5 :

 24425 Figure 2.4: Distribution of the points at the edges between different processors: processor 0 in dark blue, processor 1 in cyan, processor 2 in yellow, processor 3 in dark red.

  and whose processor number in p c is e l ∈ [[1, M e (p c )]], the relation reads: e g = e l + mesh_p1_loc%disedge(p c ) -1.

•

  for attributed cells, m l ∈ [[1, M (p c )]] -mesh_p1_loc%jj(:,m l ) which are the local indices of the cell's points, -mesh_p1_loc%neigh(:,m l ) which are the local indices of the neighbouring cells if attributed to p c , 0 if interface, -1 if not attributed, -mesh_p1_loc%jce(:,m l ) which are the global indices of the cell's edges, -mesh_p1_loc%i_d(m l ) which is the cell's global index, • for points, v l ∈ [[1, N (p c ) + N extra (p c )]] -mesh_p1_loc%loc_to_glob(v l ) which is the point global index, -mesh_p1_loc%rr(:,v l ) which are the (r, z) coordinates of the point, • for edges at the interface, m s ∈ [[1, M es (p c )]] -mesh_p1_loc%jjs(:,m s ) which are the global indices of the edge's points, -mesh_p1_loc%sides(m s ) which is the interface number, -mesh_p1_loc%neighs(m s ) which is the local index of the cell containing that edge, • for edges bordering a lower rank processor, edge e l ∈ [[1, M e extra (p c )]]

  (:,m l ) which are the global indices of the cell's points, -mesh_p1_loc%extra_jce(:,m l ) which are the global indices of the cell's edges, -mesh_p1_loc%extra_jcc(m l ) which is the cell's global index.

e 2 v n 2 e E 2 e I 2 e I 3 e I 1 m m 1 m 2 m 3 Figure 2 . 6 :

 22232326 Figure 2.6: Refinement of a cell from the left panel to the right panel. The numbering of cells, points and edges is shown, along with new edge orientations.

E i l == 2 *

 l2 M e (p c ) + e i l for the new inside edges, e r I i l M e (p c ) + m l -1 + i

2 . 5 (

 25 right panel), there are hence 5 extra cells for processor 1, 4 extra cells for processor 2, and 3 extra cells for processor 3.The requirement for the extra cells is to get the global index of each object. It is built using the local indices in their attributed processor and then applying the usual global shift done to get the global indices.

2 / 3 extra

 23 (p c)-N 1 extra (p c ). This is done so that, once again, for v l ∈ [[1, N 2/3 (p c )]], the point is attributed to the processor p c and, forv l ∈ [[N 2/3 (p c ) + 1, N 2/3 (p c ) + N 2/3 extra (p c )]], it is not.For the newly created nodes, as they are created on an edge not attributed to p c (local index e l ∈ [[1, M e extra (p c )]] in p c ), they will be attributed to the edge's processor.

Fig. 1 :

 1 Fig. 1: Preconditioning of A + λµB T M -1 Q B using either MUMPS or BoomerAMG for the preconditioner. Y-axis:wall-clock time×nb. processor nb. vel.+press. degrees of freedom in millisecond; X-axis: number of velocity grid points. µ = 1 ( ); µ = 10 -2 ( ), µ = 10 -4 ( ).

Fig. 2 :

 2 Fig. 2: GMRES residual vs. iteration count for the preconditioners CaCh 2 (left panel) and CaCh 1 (right panel) with µ = 1 (blue), µ = 10 -2 (red), and µ = 10 -4 (green). Number of grid points: 29857 ( ); 118785 ( ), 473857 ( ) ; 1892865 ( ).

Figure 3 . 1 :

 31 Figure 3.1: Quantities saved in the different Fortran scalar tables during computation. The order has been optimized for both computation cost and memory usage.

Figure 3 . 2 :

 32 Figure 3.2: 'means info' file which references which columns contain which quantities in the 'fort.75' file.

Figure 3 . 3 :

 33 Figure 3.3: Example of the error made on D b . Here D b should be equal to zero everywhere as ℓ = 0.

Figure 3 . 4 :

 34 Figure 3.4: Example of 'data_for_suites' file in order to select which fields to process for visualisation or statistical analysis.

Figure 3 . 5 :

 35 Figure 3.5: Example of the three types of visualisations with bi-logarithmic scale and isosurfaces.

Figure 3 . 6 :

 36 Figure 3.6: Example of how camera rescaling is used to easily compare quantities on two meshes with different dimensions.

  def make_animations ( file_names , data_directory , screen_directory , nb_frames , rot = True , tr = True , DDD = True , rot_cp = False , tr_cp = False , rot_frames = -1 , tr_frames = -1 , DDD_frames = -1 , rot_cp_frames = -1 , tr_cp_frames = -1 , opt_vect = False , ids_colormap = None , calculator_name = ' ' , opt_contour = ' ' , opt_p = None , opt_blade_file = ' ' , opt_tr_z = None )

  3.7).

Figure 3 . 7 :

 37 Figure 3.7: Visualisation window using the image viewer for some scalar fields.

  Tracking dynamo mechanisms from local energy transfers: application to the von Kármán sodium dynamo M. Creff, H. Faller, B. Dubrulle, J.-L. Guermond, C. Nore.

  Magnetic field produced by dynamo action

Figure 4 . 1 :

 41 Figure 4.1: (a) Schematic of the VKS experimental device of [5] in non-dimensional units. The impellers counter-rotate as indicated. (b) Magnetohydrodynamics simulation of the von Kármán Sodium dynamo at kinetic Reynolds number R e = 1500, magnetic Reynolds number R m = 150 and relative magnetic permeability of the impellers µ r = 50: visualization of the time-averaged magnetic field lines and section of the azimuthal magnetic field component. Adapted from [22].

  e.g. (b r , b z ) {Ω r },{Ω z } ------→ b θ {h θθ } ---→ (b r , b z ). These averaged Ω-terms are shown in figure 4.4. Their large structures are invariant under the R π symmetry (as expressed in equation (82)). They are similar for the two values of µ r and their maximum absolute non-dimensional value is around 10 2 .

Figure 4 . 2 :

 42 Figure 4.2: Time-averaged helicity tensor {h ij } in the saturated regime for µ r = 50.

Figure 4 . 3 :

 43 Figure 4.3: Time-averaged helicity tensor {h ij } in the saturated regime for µ r = 1.

Figure 4 . 4 :

 44 Figure 4.4: Time-averaged Ω effect in the saturated regime. Top row: µ r = 50, bottom row µ r = 1

Figure 5 . 1 :Figure 5 . 2 :

 5152 Figure 5.1: Time evolution of kinetic energy budget(19) in the growing phase. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 3 :

 53 Figure 5.3: Time evolution of renormalized magnetic budget in the growing phase. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 4 :

 54 Figure 5.4: Dynamo growth rate averaged on the 128 snapshots in the dynamo growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 5 :

 55 Figure 5.5: Kinetic energy averaged on the 128 snapshots in the dynamo growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 6 :

 56 Figure 5.6: Renormalized and time-averaged magnetic energy in the growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 7 :

 57 Figure 5.7: Time-averaged viscous dissipation in the growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 8 :

 58 Figure 5.8: Renormalized and time-averaged Joule dissipation in the growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 9 :

 59 Figure 5.9: Renormalized and time-averaged irregular Joule dissipation in the growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 .

 5 Figure 5.10: Time-averaged anomalous velocity dissipation in the growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 11 :

 511 Figure 5.11: Renormalized and time-averaged anomalous magnetic dissipation in the growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 12 :

 512 Figure 5.12: Renormalized and time-averaged energy transfer between fields in the growth regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 5 . 13 :

 513 Figure 5.13: Renormalized and time-averaged rotating magnet effect for the case R e = 1500, R m = 150, µ r = 50 in the growth regime at different filtering scales: (left) η = 0, (middle) η = 1, (right) η = 16.

Figure 5 . 14 :

 514 Figure 5.14: Renormalized and time-averaged Joule induced in-between scales magnetic transfers for the case R e = 1500, R m = 150, µ r = 50 in the growth regime.

Figure 5 . 15 :

 515 Figure 5.15: Renormalized and time-averaged Lorentz induced in-between scales magnetic transfers for the case R e = 1500, R m = 150, µ r = 50 in the growth regime.

Figure 6 . 1 :

 61 Figure 6.1: Time evolution of kinetic energy budget (19) in the saturated phase. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

figure 5 . 3 (Figure 6 . 2 :

 5362 Figure 6.2: Time evolution of magnetic energy budget (20) in the saturated phase. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 . 3 :

 63 Figure 6.3: Dynamo growth rate averaged on the 128 snapshots in the dynamo saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 . 4 :

 64 Figure 6.4: Kinetic energy averaged on the 128 snapshots in the dynamo saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 . 5 :

 65 Figure 6.5: Time-averaged magnetic energy in the saturated regime. Cases are defined in figure 6.4. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 . 6 :

 66 Figure 6.6: Time-averaged viscous dissipation in the saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 . 7 :

 67 Figure 6.7: Time-averaged Joule dissipation in the saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 . 8 :

 68 Figure 6.8: Time-averaged irregular Joule dissipation in the saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 .

 6 Figure 6.8 shows the Joule dissipation due to the irregularity of the electrical conductivity. As expected, it is located on the lateral border {r = 1.4, 0 ≤ θ < 2π, -1 ≤ z ≤ 1} where the conductivity has a jump from the stagnant liquid sodium layer to the copper vessel. For all filter scales, it remains low (maximum absolute values around 10 -4 ) and is therefore not really relevant for the dynamo action.

Figure 6 . 9 :

 69 Figure 6.9: Time-averaged anomalous velocity dissipation in the saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 .

 6 Figure 6.10: Time-averaged anomalous magnetic dissipation in the saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 .

 6 Figure 6.11: Time-averaged energy transfer between fields in the saturated regime. Panels correspond to η = 0 for the first row, η = 1 for the second row, η = 16 for the third row, and the case R e = 1500, R m = 150, µ r = 50 on the first column and the case R e = 1500, R m = 300, µ r = 1 on the second column.

Figure 6 .

 6 Figure 6.12: Time-averaged rotating magnet effect for the case R e = 1500, R m = 150, µ r = 50 in the saturated regime at different filtering scales: (left) η = 0, (middle) η = 1, (right) η = 16.

Figure 6 .

 6 Figure 6.13: Time-averaged Joule induced in-between scales magnetic transfers for the case R e = 1500, R m = 150, µ r = 50 in the saturated regime.

Figure 6 .

 6 Figure 6.14: Time-averaged Lorentz induced in-between scales magnetic transfers for the case R e = 1500, R m = 150, µ r = 50 in the saturated regime.

4 )

 4 Les différentes quantités présentes dans ces équations sont : u le champ de vitesse, p le champ de pression, b la densité de flux magnétique (ou induction magnétique) , j le courant électrique, ρ la densité du fluide, ν la viscosité du fluide, σ la conductivité électrique du milieu et µ la perméabilité magnétique du milieu. La première équation (6.1) correspond à l'équation de Navier-Stokes pour les fluides incompressibles avec l'ajout de la force de Lorentz (j × b)/ρ. La seconde (6.2) est obtenue en utilisant les équations de Maxwell dans la limite quasi-statique avec la loi d'Ohm écrite dans le référentiel du laboratoire j = σ(E + u × b) permettant alors d'éliminer le champ électrique E des équations. La troisième (6.3) traduit l'incompressibilité du fluide : l'écoulement est de divergence nulle. De même, la quatrième (6.4) provient de la conservation des lignes de champ magnétique et l'inexistence de monopoles magnétiques.Les autres relations importantes pour le champ magnétique sont :∇ × H = j,H est le champ magnétique. Ces équations proviennent des équations de Maxwell et permettent de lier le champ magnétique au courant électrique. Pour des raisons de simplicité, nous appellerons par la suite b le champ magnétique car nous ne ferons plus référence à H. L'effet dynamo est alors l'amplification et le maintien du champ magnétique b par le mouvement du fluide conducteur u. Dans notre cas, cela provient d'une instabilité. En effet, b = 0 est toujours solution de l'équation (6.2). Ainsi, pour qu'un champ magnétique soit généré et maintenu, il est nécessaire que le terme non-linéaire (∇ × (u × b)) dépasse la dissipation par effet Joule (-∇ × j σ ). Alors, pour un système de taille caractéristique L et vitesse caractéristique V , il est à minima nécessaire d'avoir :

(6. 7 )

 7 R m est connu sous le nom de nombre de Reynolds magnétique et donne une borne inférieure pour la génération d'un champ magnétique. En effet, pour un champ de vitesse fixe, l'équation (6.2) devient linéaire en b et donc R m devient le seul paramètre de contrôle. Dans ce cadre de travail avec

  6.1) est une approximation de l'expérience. Le volume interne pour r ∈ [0, 1] contient le sodium liquide alors que la partie externe r ∈ [1, 1.6] est composée de deux couches de solides conducteurs : une couche avec la même conductivité que le liquide pour r ∈ [1, 1.4] (modélisant une couche de sodium immobile) et une autre avec une conductivité 4.5 fois supérieure pour r ∈ [1.4, 1.6] (modélisant la paroi en cuivre). Les turbines ont une perméabilité différente en fonction du matériau utilisé.

Figure 6 . 1 :

 61 Figure 6.1: Représentation des différents volumes simulés : fluide en mouvement (vert), fluide au repos (jaune), autres murs conducteurs (rouge) et turbines en mouvement (rouge foncé).

( 6 . 10 )

 610 où u est le champ de vitesse, f une source thermique et K la diffusivité thermique du fluide, cela revient à résoudre pour chaque mode de Fourier m et fonction test ϕ i :A m,c ij T n+1,m,c j = B m,c i ,(6.11)137 pour la partie cosinus et un pas de temps temporel n, avec:

Figure 6 . 2 :

 62 Figure 6.2: Convention pour les points, bord et voisins d'une cellule.

Figure 6 . 3 :Figure 6 . 4 : 3 Figure 6 . 5 :

 6364365 Figure 6.3: Distribution des points sur différents processeurs : processeur 0 en bleu foncé, processeur 1 en bleu clair, processeur 2 en jaune, processeur 3 en rouge.

Figure 6 . 6 :

 66 Figure 6.6: Numérotations des points et orientation des bords pour une cellule pour les éléments finis P 1 (panneau gauche), P 2 (panneau central) et P 3 (panneau droit).

6 . 14 )terme µ 4 (

 6144 Le ∇ + ∇ T ) : (∇ + ∇ T ) correspond au Laplacien plus les termes extra-diagonaux en ∂ x ∂ y pour le cas 2D et f une force externe. Dans sa forme discrétisée, le problème de Stokes s'écrit (U, P et F étant respectivement la vitesse, la pression et la force externe, discrétisées) : AU -B T P = F, BU = 0, (6.15)

Figure 6 . 7 :

 67 Figure 6.7: Preconditionement de A + λµB T M -1 Q B pour différents préconditionneurs et solveurs. Axe des Y : temps de simulation×nb. processeurs nb. vel.+press. degrés de libertés en millisecondes ; axe des X : nombre de points pour le maillage de vitesse. µ = 1 ( ); µ = 10 -2 ( ), µ = 10 -4 ( ).

Figure 6 . 8 :

 68 Figure 6.8: Résidu GMRES vs. nombre d'itérations pour le préconditionneur

  rouge), et µ = 10 -4 (vert). Maillage 0.1M pts. vit. ( ); Maillage 0.5M pts. vit. ( ), Maillage 1.9M pts. vit. ( ); Maillage 7.6M pts. vit. ( ). Préconditionneur de gauche à droite: Cahouet&Chabard CaCh 1 (λ = 0), Cahouet&Chabard CaCh 2 (λ = 0), Lagrangien augmenté C λµ avec λ = 1, Lagrangien augmenté C λµ avec λ = 10.

Figure 6 . 9 :

 69 Figure 6.9: Préconditionnement de τ -1(B M -1 V B T ) -1 + µ(1 + λ)M -1

  Cependant, au vu des grandes quantités de données, il est important d'automatiser ce processus. Une fois les champs transférés dans les fichiers pvd/vtu, il est possible de les visualiser grâce à Paraview et d'automatiser le tout via Python. Le projet GitLab se nomme 'paraview_python' et la version Paraview utilisée est 5.8.0. Il peut y avoir quelques problèmes si des versions Paraview différentes sont utilisées car il n'y a pas de garantie que les noms des fonctions Python soient les mêmes. Comme SFEMaNS est en coordonnées cylindriques, les différentes visualisations implémentées pour des quantités scalaires sont (voir figure 6.10) : des sections planes (θ, r) pour différentes hauteurs z, des sections planes (r, z) pour différents angles θ et des isosurfaces 3D.

Figure 6 . 10 :

 610 Figure 6.10: Exemple des trois types de visualisation avec échelles et contours bi-logarithmiques.

Figure 6 . 11 :

 611 Figure 6.11: Fenêtre de visualisation pour accès à distance des différents champs scalaires.

30 )

 30 Comme la densité ρ et la viscosité du fluide ν sont constantes dans nos simulations, elles n'impactent pas le filtrage et peuvent donc en être extraites.

  2µ sont :∂ t E c + ∇ • J N S + D ν + D u + D b = -T u→b + P, (6.31) ∂ t E m + ∇ • J M + D σ + D σ + D b + D µσ + D µT = T u→b + M µ , ℓ )u -ν∇(uu ℓ )) + 1 2 (u ℓ p + up ℓ ),(6.33)

. 42 )

 42 Tous les termes mesurent les transferts d'énergie pour toutes les échelles plus grandes que ℓ. Pour commencer, nous discutons les termes de l'équation (6.31). Du côté gauche de l'équation, J N S représente le courant d'énergie responsable du mouvement spatial de l'énergie et ainsi ne participe ni aux transferts d'énergie entre échelles et champs ni aux dissipations. D ν est la dissipation due à la viscosité du fluide : c'est globalement une perte d'énergie pour le champ de vitesse. D u est la dissipation anormale dite de Duchon-Robert. Ce terme correspond au transfert local d'énergie vers les échelles plus petites que ℓ et ainsi permet de caractériser où se forment les irrégularités du champ de vitesse.

  [1, 1.4] × [0, 2π) × [-1, 1]} pour modéliser le sodium , et σ r = 4.5, µ r = 1 in {(r, θ, z) ∈ [1.4, 1.6] × [0, 2π) × [-1, 1]} pour représenter les murs en cuivre. Les deux simulations que nous utilisons sont réalisées à un nombre de Reynolds modéré R e = 1500 mais avec différentes perméabilités magnétiques pour les turbines : µ r = 1 pour celles en acier et µ r = 50 pour celles en fer doux. Pour des conditions aux limites de pseudo-vide, les nombres de Reynolds magnétique pour avoir une dynamo sont R c m (µ r = 1) = 190 ± 10 et R c m (µ r = 50) = 90 ± 5. Ainsi, nous utilisons des R m au-dessus de ces seuils, respectivement R m (µ r = 1) = 300 et R m (µ r = 50) = 150, pour pouvoir observer les deux phases de croissance et saturation de l'effet dynamo. Chaque régime est étudié pour 128 pas de temps (1 tous les huitièmes de tour pour stroboscoper le mouvement des turbines). Nous utiliserons le symbole {.} pour la moyenne faite sur tous ces pas de temps. La moyenne sur tout le volume aura le symbole langle.⟩. L'échelle de Kolmogorov pour ces simulations est η = 1.16 × 10 -2 .

Figure 6 . 12 :

 612 Figure 6.12: Évolution temporelle du bilan d'énergie magnétique (6.32) sur l'ensemble du volume renormalisée en phase de croissance. Les panneaux correspondent à ℓ η = 0 pour la première ligne, ℓ η = 1 pour la deuxième ligne

Figure 6 . 13 :

 613 Figure 6.13: Transferts d'énergie entre les champs de vitesse et magnétique renormalisés et moyennés en temps en phase de croissance. Les panneaux correspondent à ℓ η = 0 pour la première ligne, ℓ η = 1 pour la deuxième ligne

Figure 6 . 14 :-Figure 6 . 15 :

 614615 Figure 6.14: Injections d'énergie par les turbines dans le champ magnétique, renormalisées et moyennées en temps en phase de croissance pour le cas R e = 1500, R m = 150, µ r = 50 aux échelles: (gauche) ℓ η = 0, (centre) ℓ η = 1, (droit) ℓ η = 16.

  This mesh is called mesh_glob. Its M cells are numbered from [[1, M ]], its N points are numbered from [[1, N ]] and its M e edges are numbered from [[1, M e ]]. Mesh structure information is stored in a number of tables and integers.

  with associated edges with local indices {e i l ∈ [[1, M e (p c )]]|i ∈ {1, 2, 3}} and with associated points with local indices {v i l ∈ [[1, N (p c )]]|i ∈ {1, 2, 3}} with all of them being attributed to the same processor p c ∈ [[1, N p ]].

Table 2 :

 2 3.6)). Lagrange elements for the velocity), using BoomerAMG, strong threshold 0.7.

		µ	Nb. points Nb. Proc.	118837 2	474065 8	1893697 32	7569665 128	30268417 512
			GMRES iter.	4	3	2	2	2
	A + λµD	1 10 -2	L2 vel. err. Times (s) GMRES iter. L2 vel. err. Times (s)	0.188E-06 12.3 4 0.178E-06 8.10	0.156E-07 61.5 3 0.151E-07 37.8	0.136E-08 176 2 0.134E-08 126	0.154E-08 529 2 0.118E-09 265	0.531E-08 1221 2 0.101E-08 745
		10 -4	GMRES iter. L2 vel. err.	4 0.172E-06	3 0.149E-07	2 0.132E-08	2 0.117E-09	2 0.104E-10
			Times (s)	3.50	7.85	21.6	36.4	71.0
		µ	Nb. points Nb. Proc.	118837 2	474065 8	1893697 32	7569665 128	30268417 512
			GMRES iter.	4	3	2	3	4
		1	L2 vel. err.	0.188E-06	0.156E-07	0.133E-08	0.123E-09	0.127E-10
			Times (s)	4.76	8.82	25.4	103	338
	A	10 -2	GMRES iter. L2 vel. err.	3 0.178E-06	3 0.151E-07	2 0.131E-08	3 0.115E-09	3 0.282E-10
			Times (s)	3.21	6.75	16.2	51.7	116
		10 -4	GMRES iter. L2 vel. err.	3 0.172E-06	2 0.149E-07	2 0.130E-08	2 0.115E-09	2 0.310E-10
			Times (s)	2.34	2.58	5.20	8.11	15.7

Preconditioning of A + λµB T M -1 Q B by A + λµD, (with τ = 1, λ = 1

, and P 2

Table 3 :

 3 Lagrange elements for the velocity), using BoomerAMG, strong threshold 0.7.

		µ	Nb. points Nb. Proc.	118837 2	474065 8	1893697 32	7569665 128	30268417 121053185 512 700
	MV + µLV τ -1	1 10 -2 10 -4	GMRES iter. L2 vel. err. Times (s) GMRES iter. L2 vel. err. Times (s) GMRES iter. L2 vel. err.	7 0.188E-06 0.156E-07 0.133E-08 0.107E-09 0.321E-10 0.325E-10 5 4 2 1 2 3.76 3.41 5.01 3.24 2.37 12.4 7 5 4 2 2 2 0.178E-06 0.151E-07 0.131E-08 0.103E-09 0.246E-10 0.335E-10 3.64 3.83 4.97 3.18 3.50 11.3 7 5 4 2 2 1 0.172E-06 0.149E-07 0.130E-08 0.103E-09 0.350E-10 0.361E-10
			Times (s)	3.11	3.07	4.67	2.85	3.16	7.20

Preconditioning of A + λµB T M -1 Q B by A, (with τ = 1, λ = 1, and P 2

Table 4 :

 4 Preconditioning

Table 5 :

 5 Preconditioning

		µ	Nb. points Nb. Proc.	118837 2	474065 8	1893697 32	7569665 128	30268417 512
			GMRES iter.	7	5	4	2	X
	MV + µLV τ -1	1 10 -2	L2 vel. err. Sol. times (s) LU fact. times (s) GMRES iter. L2 vel. err. Sol. times (s) LU fact. times (s) GMRES iter.	0.188E-06 1.25 5.42 7 0.178E-06 1.23 5.51 7	0.156E-07 1.53 13.1 5 0.151E-07 1.31 13.0 5	0.133E-08 4.16 68.8 4 0.131E-08 3.55 51.3 4	0.106E-09 4.77 209 2 0.105E-09 4.36 160 2	X X X X X X X X
		10 -4	L2 vel. err. Sol. times (s)	0.172E-06 1.20	0.149E-07 1.26	0.130E-08 3.58	0.105E-09 4.54	X X
			LU fact. times (s)	5.01	11.3	51.9	206	X

  Magnetohydrodynamics equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Mean field theory: α and Ω effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Filtered equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Local MHD coarse-grained energy budget at finite scale . . . . . . . . . . . . . . . . 3.3 Exact local MHD energy budget . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Definition of a local dynamo growth rate . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Connection with previous works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Application to the VKS dynamo 4.1 Numerical setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 SFEMaNS code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Description of the numerical runs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Scale filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Analysis using the mean-field framework . . . . . . . . . . . . . . . . . . . . . . . . . . Using various conductive materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents
	1 Introduction
	2 Background material and general formalism
	2.1 3 Local Energy budget
	3.1 7 Mechanisms for dynamo action
	8 Conclusion
	9 Acknowledgement
	A Navier-Stokes calculations
	B MHD calculations
	B.1

5 Dynamo growth mechanisms from local energy budgets 5.1 Time evolution of energy terms in the dynamo growth phase . . . . . . . . . . . . . . . 5.2 Spatial distribution of linear growth rate . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Kinetic and magnetic energy distributions in the growth regime . . . . . . . . . . . . . 5.4 Energy dissipation in the growth regime . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Local energy transfers between scales in the growth regime . . . . . . . . . . . . . . . . 5.6 Conversion of energy in the growth regime . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Dissipation and energy transfer through variation of magnetic permeability . . . . . . 6 Saturation of the dynamo from local energy budgets 6.1 Time evolution of energy terms in the dynamo saturated phase . . . . . . . . . . . . . 6.2 Spatial distribution of the local growth rate in the saturated regime . . . . . . . . . . . 6.3 Kinetic and magnetic energy repartition . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Energy dissipations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Local energy transfers between scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Conversion of energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Dissipation and energy transfer through variation of magnetic permeability . . . . . . B.2 Using various permeable materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C Linking kinetic energy transfers and velocity increments D Linking magnetic energy transfers and velocity increments 1 Introduction

  Pour simuler les équations de la magnétohydrodynamique et calculer les différents transferts d'énergie, nous utilisons le code SFEMaNS. Cette section est un résumé du chapitre 2 de ce manuscrit donnant une brève description du code SFEMaNS ainsi que quelques améliorations que nous y avons apportées. Après une description des approximations numériques utilisées, nous présentons une nouvelle méthode pour répartir la grille de simulation sur les processeurs. Finalement, nous étudions une autre méthode pour simuler les équations. Tout ce travail a été fait dans le cadre d'une collaboration avec J.-L. Guermond qui m'a invité à Texas A&M University du 5 Janvier au 3 Février 2022 et du 15 Juillet au 3 Décembre 2022. La code SFEMaNS [59] est écrit en Fortran-90 et est dédié aux systèmes axi-symétriques. Il utilise les coordonnées cylindriques et décompose alors les champs sur une base d'éléments finis dans le plan (r, z) et en mode de Fourier suivant θ. Cela permet une grande parallélisation grâce à l'utilisation du module MPI pour les modes de Fourier et de PETSC pour la parallélisation en domaines sur les éléments finis.

  6.4) composée de triangles 2D telle que, quand la taille caractéristique des cellules tend vers zéro, les fonctions tests génèrent l'espace entier des fonctions. Les grilles de simulations pour les champs magnétique et de vitesse sont en éléments P 2 et pour la pression P 1 . Cela est dû au fait que seules certaines paires d'éléments permettent d'assurer la convergence vers des solutions propres de Navier-Stokes. Ces combinaisons permettent alors d'écrire toute équation sous la forme d'un problème simple d'algèbre linéaire. Par exemple, pour l'évolution temporelle du champ de température pour un fluide :

  .16) Il est alors possible de résoudre le problème en pression indépendamment de la vitesse. La difficulté pour inverser S (et donc résoudre le problème) provient du fait que cette matrice est très loin d'être diagonale et prend donc beaucoup d'itérations pour des solveurs itératifs comme GMRES et MINRES.De plus, comme il n'est pas possible de construire les éléments de la matrice S pour des raisons de mémoire, on ne peut pas utiliser des méthodes de résolution plus directes. Mathématiquement, la quantité d'intérêt est le nombre de conditionnement ℓ 2 de S étudié dans[START_REF] Shih | Robust multigrid techniques for augmented Lagrangian preconditioning of incompressible Stokes equations with extreme viscosity variations[END_REF] Prop. 5.24] et dans[START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF] Prop. 50.14].

Table 6 .

 6 1: Options BoomerAMG utilisées pour tester les préconditionneurs.

	Choix

du préconditionneur pour le problème en vitesse

  λµB T M -1 Q B)X = Y .Pour accélérer cette inversion, un préconditionnement est nécessaire et nous étudions les préconditionneurs suivants :

	Pour résoudre le problème en pression, il est nécessaire d'inverser la ma-
	trice A + λµB T M -1 Q B. Cependant, il n'est pas possible de construire les co-

efficients de cette matrice et il est donc nécessaire d'utiliser une méthode itérative. Il s'agit alors de résoudre un problème en vitesse pour X : (A +

In the following, we will abusively call b the magnetic field.

This will be the case for example in the von Kármán setup if the impellers are made of permeable material such as, for example, soft iron[5].

Remerciements

This material is based upon work supported in part by the DMS-2110868 (JLG), by the Air Force Office of Scientific Research, USAF, under grant/contract number FA9550-23-1-0007 (JLG), the Army Research Office, under grant number W911NF-19-1-0431 (JLG), and the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contracts B640889, B641173 (JLG). The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. http://www.tacc.utexas.

By default the color scale is set as bi-logarithmic as the values of the studied fields greatly vary between the bulk and the turbines. Its maximal and minimal values are set by the parameter 'p' which sets them to ±10 p . Other types of color scales are available and can be set using 'id_colormap'. Creating a new one requires coding its values in 'color_and_opacity_maps.py'. For example, the bi-logarithmic is coded as:

[ -10 ** p , 0 .0 , 0 .0 , 0 .4 , -10 ** ( p -1 ) , 0 . , .5 , 1 . , -10 ** ( p -2 ) , 0 .0 , .9 , 1 .0 , -10 ** ( p -3 ) , 0 .0 , 1 . , .9 , 0 .0 , 1 .0 , 1 .0 , 1 .0 , 10 ** ( p -3 ) , 1 .0 , 1 .0 , 0 .0 , 10 ** ( p -2 ) , 1 .0 , 0 .5 , 0 .0 , 10 ** ( p -1 ) , 1 . , 0 .0 , 0 .0 , 10 ** p , 0 .4 , 0 .0 , 0 . 0 ] where each line corresponds to a new color being set at a specific value of the field. The same goes for opacity scales which are used for 3D isosurfaces.

Isosurfaces are created for all of them. By default, their values depend on the selected color scale (a logarithmic color scale will make logarithmic isosurfaces). Others can be set using 'opt_contour' and new ones created in 'contour_definitions.py'. Here are all the ones currently present:

' positive_linear_tenth ': [ 1 . , .9 , .8 , .7 , .6 , .5 , .4 , .3 , .2 , . 1 ] , ' linear_tenth ': [ 1 . , .9 , .8 , .7 , .6 , .5 , .4 , .3 , .2 , .

1 , 0 , -.1 , -.2 , -.3 , -.3 , -.4 , -.5 , -.6 , -.7 , -.8 , -.9 , -1 . ] , ' positive_linear_ddd ': [ 1 . , .9 , . 75 , . 5 ] , ' linear_ddd ': [ 1 . , .9 , . 75 , .5 , -.5 , -. 75 , -.9 , -1 . ] , ' positive_log_10 ': [1 , .1 , . 01 , . 001 ] , ' log_10 ': [1 , .1 , . 01 , . 001 , 0 , -1 , -.1 , -. 01 , -. 001 ] , ' positive_log_10_ddd ': [1 , .1 , . 01 ] , ' log_10_ddd ': [1 , .1 , . 01 , -1 , -.1 , -. 01 ] , ' 3waygreater ': [ 0 .9 , -0 . 9 ] }

The ones ending with 'ddd' are used for 3D visualisations as too many 3D isosurfaces make the pictures too bloated.

One last important option, specific to the translation animation, is 'opt_z' which sets the various heights at which the visualisation will be done. It is very practical in order to study what takes place just over and under the blades and impellers.

Specific animations have been developed for vector fields ('animation_translation_compact' and 'animation_rotation_compact'). They do the same thing as the previous

A Navier-Stokes calculations

Using the scalar product of ρu/2 with equation ( 15) and the scalar product of ρu /2 with equation (??) leads to the equation for the transferred kinetic energy from all scales larger than to the scale :

After some calculations, the left hand side yields a simpler form (using multiple times ∇ • u = 0 and ∇ • u = 0):

And the first part of the right hand side can then be written as:

B MHD calculations

Using the scalar product of b/2µ with equation ( 16) and the scalar product of b /2µ with equation ( 2) leads to the equation for the transferred magnetic energy from all scales larger than to the scale :

The variation of the physical properties σ and µ leads to extra terms that we develop in the following.

B.1 Using various conductive materials

The sharp variation in the electrical conductivity at {r = 1.4R cyl , -H/2 ≤ z ≤ H/2} between the liquid sodium lateral layer (σ = σ 0 ) and the copper wall of the container (σ = 4.5σ 0 ) needs to be taken into account:

Using:

Then:

Hence:

B.2 Using various permeable materials

The sharp variation in the magnetic permeability between the liquid sodium and the moving impellers (when µ = µ 0 ) needs to be taken into account:

Hence :

C Linking kinetic energy transfers and velocity increments

Following [26], the filtering process discussed in section § 4.4 can be interpreted as a smoothing operation at a resolution achieved through the function ϕ (smooth non-negative function with unit integral) such that the filtered velocity field reads:

Using increments of the velocity fields, the third-order structure function can be written as (see equation (2.7) in [27]):

with J S N S = -uu 2 + uu 2 . The term S N S was denoted D I in [26]. Note that it differs from the anomalous dissipation term D u by an additional flux term (in divergence). Integrated throughout

D Linking magnetic energy transfers and velocity increments

In the MHD case with µ = const., we introduce:

In the special case µ = const., this formula simplifies like in [24]:

Treating both parts of the integral independently:

Grouping two by two similar terms:

using that for fields of null divergence: