Relaxation to magnetohydrodynamics equilibria via collision brackets
Metriplectic dynamics is applied to compute equilibria of fluid dynamical systems. The result is a relaxation method in which Hamiltonian dynamics (symplectic structure) is combined with dissipative mechanisms (metric structure) that relaxes the system to the desired equilibrium point. The specific metric operator, which is considered in this work, is formally analogous to the Landau collision operator. These ideas are illustrated by means of case studies. The considered physical models are the Euler equations in vorticity form, the Grad-Shafranov equation, and force-free MHD equilibria.
The computation of general 3D magnetohydrodynamic (MHD) equilibria plays a fundamental role in simulations of stellarators and it is important for tokamaks as well, due to deviations from axisymmetry (magnetic islands, ripples, and resonant magnetic perturbations).
This problem has always attracted interest in the plasma physics community, leading to different numerical approaches [1, 2, 3, 5, 4]. Nonetheless, the efficient computation of three-dimensional MHD equilibrium is still an open issue.
In the geometric mechanics community, on the other hand, there has been a significant effort directed to the study of the appropriate geometric structures for the description of dissipative systems and irreversible dynamics. Such a structure has been proposed by Morrison [12, 13] and it is referred to as metriplectic dynamics since it combines the symplectic structure of Hamilton’s equations with the metric structure of gradient flows. Several physical systems can be cast in metriplectic form, e.g., free rigid body with suitable chosen torque , resistive MHD  and the Lindblad equation for open quantum systems . However, the geometric properties of metriplectic flows can also be exploited to design artificial dynamical systems that relax to an equilibrium of the considered physical system. The advantages of such methods come from properties directly implied by the geometric structure and the energy-Casimir principle [14, 10, 11].
A related but different approach has been proposed by Flierl and Morrison  and developed further by Chikasue, Furukawa and Morrison [6, 7, 8, 9]. In such an approach, the relaxation mechanism is constructed on the basis of the symplectic structure only, essentially by squaring the Poisson brackets. This method has the properties of minimizing the energy functional of the system, while preserving all the other Casimir invariants; for the case of ideal MHD this implies the preservation of the topology determined by the initial conditions.
The present work is set in the framework of metriplectic dynamics. A specific metric operator is constructed on the lines of the Landau operator for Coulomb collisions, which has a metriplectic structure already discovered by Morrison [12, 13]. The basic idea is developed for three case studies in order to explore the advantages and disadvantages of the proposed relaxation method.
For a class of dynamical systems arising in fluid and kinetic theory of plasmas, equilibrium states can be characterised by a variational principle.
Typically equilibrium states are extrema of an entropy functional under the constraints imposed by the first integrals of the systems, such as mass, energy or momentum. For instance, in the Boltzmann equation (with collisions) equilibria are obtained by extremizing the entropy at constant energy, momentum, and particle number . Moreover, the Boltzmann equation has the additional property that a solution of the initial value problem relaxes, as time goes to infinity, to an equilibrium because of the celebrated H theorem ; equilibria can therefore be identified by time-evolution of properly chosen initial conditions. In general this is not the case: ideal systems with no dissipation mechanisms will not relax to an equilibrium. Therefore, in order to design a relaxation method for the computation of equilibria, some dissipation mechanism has to be introduced.
The idea proposed by Morrison  shows the possibility to define a dissipative dynamics that relaxes to a solution of the variational problem for the equilibrium. These concepts will be explained in more detail here with the help of specific physical models.
2.1 Physical Models
Three specific case studies are presented in order to illustrate the proposed idea.
The first example uses the vorticity form of the 2D Euler equations,
where is time and are the Cartesian coordinates in the two-dimensional space. The dynamical variable is the vorticity of an incompressible flow in two-dimensions and is the stream function. The two-dimensional Laplacian is and , for any pair of functions .
The equilibrium states of the Euler system are reached when . Then from , the vorticity must be proportional to a function of the stream function, . Substituting this expression into the Poisson equation of (1) leads to the non-linear eigenvalue problem
for the pair ; given a solution, the corresponding vorticity field is determined by . Complemented with boundary conditions, the nonlinear eigenvalue problem (2) determines a whole class of equilibrium states : for every choice of , each solution corresponds to an equilibrium.
where and are the radial and axial coordinates of a cylindrical reference system ; is the speed of light in vacuum (c.g.s. units). Physically the unknown is a flux function and the right-hand side is with the -component of the current density.
The Grad-Shafranov equation is formally analogous to equation (2), if the Laplace operator and the stream function are replaced by and the flux function , respectively. The same considerations about the equilibrium states apply.
As a last example, force-free MHD equilibria (also known as Beltrami fields or Taylor-relaxed states) are considered. The magnetic field , where are the Cartesian coordinates in three-dimensional space, satisfies the force-free equilibrium condition if
where is in general a function. If is constant, the Beltrami equation reduces to the eigenvalue problem for the curl operator.
2.2 Variational Principle
The equilibria of the considered systems can be characterised as the extrema of an entropy functional with the constraint of constant Hamiltonian.
We consider either the case of a scalar field or a multicomponent field , defined over a spatial domain . Let and be the entropy and Hamiltonian functionals, respectively. The problem of finding the extrema of at constant is written as
where is the Lagrange multiplier.
In the Euler example, the dynamical variable is the scalar vorticity . The energy functional is the kinetic energy of the fluid (per unit mass) and it is written as
We restrict the entropy functional to be of the form
where is smooth and with derivative being a monotonic function. The functional derivatives are readily computed,
and equation (5) becomes
One can now compare this result with equation (2), and deduce a relationship between the choice of the entropy functional and a particular physical equilibrium described by the function , namely,
The inverse exists since is monotonic. This also implies that only equilibria with a monotonic can be described in this way.
In the case of the Grad-Shafranov equation, the natural variational principle  seeks the extrema of an action functional written in terns of a Lagrangian density. However, this is not in the form (5). In order to obtain a variational principle in the form (5), let us introduce the variable . The flux function is determined from by solving the linear elliptic problem
equipped with the desired boundary conditions. Then we define the energy functional
and we consider functionals of the form
As an example, the entropy
where and are positive constants, leads to the Herrnegger-Maschke solutions of the Grad-Shafranov equation ( and reference therein). Other choices, which will lead to different physical equilibria, can be made. The functional derivatives are
where . Equation (5) becomes
Using the entropy (11) for the sake of illustration, one obtains
and since , one obtains the weighted linear eigenvalue problem
that characterizes the Herrnegger-Maschke solutions.
As a last example, let us address Beltrami fields, that is the minimizers of the magnetic energy at constant magnetic helicity [30, 31]. Thus the Hamiltonian functional should be the magnetic helicity,
where is the vector potential. We fix the Coulomb gauge, so that
plus suitable boundary conditions. This establishes a one-to-one relationship between and . Correspondingly, the dissipated entropy is actually the physical energy of the magnetic field
We have now two equivalent choices. The standard choice consists in setting and computing
and condition (5) gives the Beltrami equation directly with a constant . Alternatively, one can set , so that
and condition (5) reduces to
which should be solved together with (15). Since is a constant, this formulation is equivalent to the Beltrami equation. We shall see that the latter choice is more convenient in terms of computational cost.
2.3 Metriplectic dynamics
The metriplectic formulation of the dynamics of a (possibly multi-componet) field reads: for every functional , the function must satisfy
where is a Poisson bracket, that is, an anti-symmetric, bilinear operation on the functionals, satisfying the Leibniz and Jacobi identities , whereas is a metric bracket, that is, a symmetric, bilinear operation with a definite sign. In the following we shall assume that the metric brackets are negative semi-definite, but this is just a convention. The functional plays the same role as the test-function in a weak formulation.
The Hamiltonian and entropy functionals must satisfy the conditions
for all . Such combatibility conditions imply
that is the entropy is dissipated at constant Hamiltonian. The qualitative idea is that, if is bounded from below, the system will evolve on the manifold , where is the initial condition, toward a state that satisfies . If the metric brackets vanish only in the direction of the Hamiltonian functional, i.e.,
then the relaxed state is a solution of the variational problem (5). This is not always the case: some metric brackets have a larger “null space” so that the set of relaxed states strictly contains the solutions of (5). If this happens we say that the operator does not completely control the relaxation process.
2.4 The metric operator and its applications to equilibria
In this work we shall focus on the metric part of the dynamics and consider the class of metric brackets introduced by Morrison  as a generalisation of the Landau collision operator. Such operators will be referred to as integral collision-like operator.
This general form for two arbitrary functionals and can be written as 
where and is a matrix with either scalar- or matrix-valued entries, depending on whether is a scalar or a multi-component field, respectively. Symmetry requires .
In order to ensure the conservation of a given Hamiltonian we choose the kernel of the metric brackets according to
where . A rigorous proof of (18) for this class of operators is still not available. We shall however show in numerical experimets that the corresponding dynamics relaxes to a solution to (5) as desired.
Nonetheless this choice of the metric brackets leads to integral operators that are as challenging as the full Landau collision operator. Even though structure-preserving methods for the discretization of such operators are now available [17, 18], we have introduced a simplified class of brackets leading to diffusion-like operators. Specifically we define
where is an effective diffusion coefficient with .
3 Computational Aspects
For the time discretization the Crank-Nicolson scheme  has been chosen, in order to guarantee discrete energy conservation, at least for quadratic energy functionals. The discrete entropy is proven to be dissipated. A finite element discretization has been chosen for the spatial operators. Both the integro-differential operator and its local version have been implemented in FEniCS, a computing platform for solving partial differential equations [19, 20].
4 Numerical Experiments
A gallery of different numerical experiments, for the models of section 2, is presented here.
4.1 2D Euler Equation
The setup for the simulations is as follows: an anisotropic Gaussian is chosen as initial condition and Dirichlet homogeneous boundary conditions are applied for the domain . The entropy functional is quadratic in the dynamical variable , namely , from which the variational principle predicts a linear relation between and the stream function , . A local version of the collision-like metric operator, derived from equation (19), has been used to evolve the system. The energy functional is verified to be preserved with a relative precision of .
Figure 1 shows the color map of the vorticity field with the contours (represented with white solid curves) of the stream function at equilibrium, that is, at the end of the simulation.
In order to verify whether the equilibrium condition selected by the choice of the entropy functional has been reached or not, it is useful to consider specific diagnostics. First of all, the time evolution of the entropy functional gives information about the relaxation process. When no appreciable variation of the entropy functional occurs, the system has reached an equilibrium. Moreover, it is interesting to make use of another diagnostic which in the following will be called a “scatter plot”. It is constructed by plotting for every grid node the corresponding discrete values and : when the system is far from the equilibrium, these points are uniformly distributed on the cartesian plane. On the other hand, as the system relaxes to the equilibrium, they show a functional relation, which can then be compared with what is theoretically expected from the variational principle. The evolution of the entropy functional and the scatter plot of the system is a key diagnostics in order to understand whether an equilibrium condition has been reached or not.
Figure 2 shows the scatter plot at the beginning and at the end of the system temporal evolution, together with the time evolution of the entropy functional.
A fit of the functional relationship is performed, confirming that a linear functional relationship between and is found. In the case of the choice of a quadratic entropy functional it is also possible to have a comparison with analytical results, thus providing a good verification for the method. In fact equation (2), which describes the physical equilibria, reduces to a linear eigenvalue problem , which can be solved analytically. The analytical eigenvalues can be compared with the result of the fit in the numerical simulation.
The eigenvaues in a domain with Dirichlet boundary conditions are
In our example , and the eigenvalue corresponding to the fundamental state is given by . The linear fit of figure 2 is , in good agreement with the analytical result.
The numerical experiments conducted confirm that the equilibrium reached is independent of the initialisation chosen for the simulation, being driven by the choice of the entropy functional only. All the numerical results are found in agreement with the analytical solution.
The same test case has been simulated in a complicated mesh, with the domain constructed from a unitary circle mapped by the Czarny mapping .
Figure 3 shows the color map of the dynamical variable and the solid white lines representing the contours of the streaming function at equilibrium.
In figure 4 the scatter plot at the beginning and equilibrium state of the simulation and the time evolution of the entropy functional are shown. Again it is to be noted that as the system reaches the equilibrium state, i.e. as the entropy functional is minimized, the expected linear functional relationship between and appears.
The more complicated domain does not influence the selection of the physical equilibrium, which is due to the choice of the entropy functional only. The method can thus be applied in complicated domains.
4.2 Grad Shafranov
Numerical simulations for the Grad Shafranov model have also been performed. The common simulation setup used here is the following: anisotropic Gaussian as initial condition, homogeneous Dirichlet boundary conditions on a rectangular domain . The entropy functional has been chosen according to the equilibrium to be simulated. As in the simulations for the Euler case, the energy functional is preserved with up to a relative precision of .
In figure 5 the color map of the dynamical variable is shown together with the contours of the flux function . The profiles of the two fields are aligned, as the simulation is at equilibrium.
A more quantitative analysis of the equilibrium state is shown in figure 6. Here the functional relationship between the dynamical variable and the flux function is presented at two different simulation times: at the beginning and at convergence. An inset shows the time evolution of the entropy functional, which is minimized. The functional relationship at convergence is the one expected from the variational principle.
An entropy functional quadratic in the dynamical variable , , has been chosen in order to run the same simulation on a Czarny mapped domain . The variational principle still describes a linear functional relation between and the flux function , .
Figure 7 shows the color map of with the white solid lines representing the contours of .
Figure 8 shows the scatter plot diagnostics at the beginning and at convergence. As the entropy functional is minimized, the functional relation collapses to a linear function of the two variables and , as expected from the variational principle.
We acknowledge useful discussions with Yaman Güçlü and Edoardo Zoni regarding the Euler equilibrium domain mapping.
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
-  Hirschman S P and Whitson J C 1983, Physics of Fluids, 26, 3553
-  Hirschman S P and Whitson J C 1986, Computer Physics Communications, 43, 143
-  Reiman A H and Greenside H S 1986, Computer Physics Communications, 43, 157
-  Hudson S R, Dewar R L \etal. 2012, Plasma Physics and Controlled Fusion, 54, 014005
-  Hirshman \etal. 2011, Phys. Plasmas, 18, 062504
-  Chikasue Y and Furukawa M 2015 Physics of Plasmas, 22, 02251
-  Chikasue Y and Furukawa M 2015 Journal of Fluid Mechanics, 744, 443-59
-  Furukawa M and Morrison P J 2017 Plasma Physics and Controlled Fusion, 59, 054001
-  Furukawa M, Morrison P J, Watanabe T and Ichiguchi T 2018, Physics of Plasmas, 25 082506
-  Gay-Balmaz F and Holm D 2013, Nonlinearity, 26, 495-524
-  Gay-Balmaz F and Holm D 2014, Nonlinearity, 27, 1747-1777
-  Morrison P J 1984, Physics Letters A, 100, 423-7
-  Morrison P J 1986, Physica D, 18, 410-9
-  Morrison P J 1998, Reviews of Modern Physics, 70, 2
-  Flierl G R and Morrison P J 2011, Physica D, 240
-  Lenard A 1900, Annals of Physics, 3, 390-400
-  Hirvijoki E and Adams M F 2017, Physics of Plasmas, 24, 032121
-  Kraus M and Hirvijoki E 2017, Physics of Plasmas, 24, 102311
-  Martin Alnæs S \etal. 2015. Archive of Numerical Software 3, 100, 9-23
-  Logg A, Mardal K A, Wells G N \etal. 2012, Springer
-  Crank J, Nicolson P 1947, Proc. Camb. Phil. Soc., 43 (1), 50-67
-  Bloch A M, Morrison P J and ratiu T S 2013, Springer Proc. in Mathematics & Statistics, 35, pp. 371–415
-  Materassi M and Tassi E, Physica D, 241, 6
-  Mittnenzweig M and Mielke A 2017, Journal of Statistical Physics, 167, 205-233
-  Materassi M and Morrison P J 2018, Journal Cybernetics and Physics , Accepted
-  Mc Carthy P J 1994, Physics of Plasma, 6, 3554
-  Shafranov V D 1958, Sov. Phys. JETP, 6, 545
-  Grad H and Rubin H 1958, Proc. of the Second United Nations Conference on the Peaceful uses of Atomic Energy United Nations, Geneva vol. 21, p.190
-  Lao L L, Hirschman S P, Wieland R M 1981, American Institute of Physics, 24, 1431
-  Woltjer L 1958, Proceedings of the National Academy of Sciences, 44, 6
-  Taylor J B 1974, Phys. Review Lett., 33, 19, 1139-1141
-  Czarny O and Huysman G 2008, Journal of Comp. Phys., 16, 227, 7423-7445