Symbolic Computations in Simulations of Hydromagnetic Dynamo

The compilation of spectral models of geophysical fluid dynamics and hydromagnetic dynamo involves the calculation of a large number of volume integrals from complex combinations of basis fields. In this paper we describe the automation of this computation with the help of systems of symbolic computations.


Introduction
Many problems of geophysical fluid dynamics reduce to the solution of the Navier-Stokes equations in a rotating coordinate system in an area having spherical symmetry [1]. This area is usually a sphere or a spherical shell. In the most general case we can speak of a spherical shell (inner radius is zero in sphere case). The research of the planetary or stellar dynamo suggests the study of an even more complex system describing magneto-convection. For these systems, exact solutions are unknown and various numerical methods are being intensively used.
The spectral approaches are very popular, when solutions of the non-stationary problem are represented as a linear combination of a finite number of stationary basis fields (modes) with timedependent amplitudes. For these amplitudes, using one of the weighted residuals methods, a dynamic system is composed. This system together with a set of modes, is simplified model of the initial spatially distributed hydrodynamic or magnetohydrodynamics system [2]. The implementation of this technique involves first construction the dynamic system, and then its numerical study. This system will be quadratically nonlinear with constant coefficients, which are determined by the values of the control parameters of the model and the selected modes.
It should be noted that the construction of the system may require much greater computational effort than its subsequent solution. This is due to the fact that the calculation of the coefficients of the model reduces to the calculation of volume integrals over the spherical shell from various multiplicative combinations of the modes themselves and/or the operator ∇. The number of these coefficients increases in proportion to the third power of the number of modes used. Even for very rough models, including ∼ 10 modes, the coefficient count goes to thousands.
However, most of the coefficients of system are zero. The reason is that the spherical symmetry of the region determines the representation of the modes in latitudinal-longitude direction with the e-mail: gvodinchar@ikir.ru help of spherical harmonics. The orthogonality of the spherical harmonics and some related functions on the surface of the sphere determines the vanishing of a large number of coefficients. Sometimes this orthogonality can be shown in a general form, sometimes it is veiled and manifested only in the process of integration, when the integral over the surface of the sphere is first computed, and then along the radius.
The value of every coefficient is the measure of interaction of the corresponding modes in the process, described by term of equation. For this reason exact equality to zero some coefficients, established by analytic integration, can carry considerable physical information. To automate the process of compiling the described dynamic systems, the systems of analytical computations are convenient. In any such system, we can programmatically implement the following calculation scheme: 1. Analytical expressions are construct for each mode.
3. Analytical integration over the sphere surface is performed. This integration gives zero or some expression depending on the radial variable; 4. If the integral over the sphere surface is non-zero, analytic or numerical integration along the radius is performed.
It is clear that the third and fourth point can be performed by standard means of any system of analytical calculations, because any such system allows you to do analytical and numerical integration. To implement the first and second items, it is necessary to calculate the operations of vector analysis in analytical form in spherical coordinates. For this, it is convenient to use systems of analytical computation, containing vector analysis command packages. For example, the Maple system contains powerful package of such calculations VectorCalculus [3,4]. However, the necessary minimum set of vector analysis operations is easy to programmed in any system. Now, we consider the described range of problems and their solution in more detail

The original equations and basic systems
The motion of an incompressible fluid in a rotating spherical shell Ω, typical of geophysical fluid dynamics, is described by the Navier-Stokes equations in a dimensionless form [5]: where v is the velocity filed, p is the effective pressure (including the centripetal term), f is the mass density of the external forces (the specified field), e z is the unit vector of the rotation axis, Re is the Reynolds number, and Ro is the Rossby number.
In the future, we always use the standard connected Cartesian (x, y, z) and spherical (r, θ, ϕ) coordinate systems, whose origin O is placed in the center of the shell Ω, and the z-axis is directed along the rotation axis.
The system (1) is supplemented by the boundary conditions for the fields v and p on the inner r = r i and outer r = 1 boundaries of the shell. In case r i = 0, the internal conditions assume the requirements for the finiteness of the fields at the center. The concrete form of the boundary conditions does not yet play a role, we shall only consider them linear and homogeneous.
In the problems of planetary and stellar dynamo, the system (1) is supplemented by equations for the magnetic field [6]: where Rm is the magnetic Reynolds number, and on the right side of the first of equations (1) there is also a Lorentz term (∇ × B) × B. For magnetic field homogeneous linear boundary conditions are also formulated.
Since the velocity and magnetic induction fields has zero divergence, they can be represented as a sum of the toroidal and poloidal components [1], and the inhomogeneity of the fields in the (θ, ϕ)direction can be expanded on a spherical harmonics. Then the toroidal basis modes for velocity and magnetic induction will have the form: and the poloidal basis modes will have the form: In these formulas Y m n (θ, ϕ) are the spherical harmonics, r is the radius-vector, the index k corresponds to the quantization of the spatial spectrum in the r -direction, and the specific form of the radial functions u T knm , b T knm , u P knm , b P knm can be different. For example, they can be chosen so that the modes (3, 4) are modes of free dissipation of the fluid and the magnetic field in the shell. Also, the form of the radial functions depends on specific boundary conditions. We shall further assume that these functions are known.
The solutions of problem (1,2) are construct in the form of linear combinations of modes (3, 4) with time-dependent amplitudes. The solenoidal nature of the velocity and magnetic induction fields, as well as the boundary conditions, will be satisfied automatically. We shall further exclude the pressure field from the solution. Also we will assume that in the field of external forces the temporal and spatial variables are separated, i.e. f (r, t) = ξ(t)f 0 (r).
So, we will use decompositions where v i (r) and B i (r) are any of the fields (3, 4), and the multi-index i = (ε, k, n, m), ε ∈ {T, P}.

The equations of the model
We supplement the magnetohydrodynamics system (1, 2) with the Lorentz term and used decompositions (5). To obtain the model equations in accordance with the methods of weighted residuals, it is necessary to multiply these equations by some weight functions and integrate over the Ω [2]. Let us consider this using the example of the Galerkin method, when the basis modes are used as weight functions.
First we consider the question of excluding the pressure field. In the most general case, it is sufficient to simply apply the rotor operation to both parts of the Navier-Stokes equation. Clearly, this will give not just a transition to the equation-consequence, but it will be an equivalent transformation. However, in the case using of zero boundary conditions for the velocity field, the pressure disappear as a result of the Galerkin procedure itself. Indeed, multiplying the Navier-Stokes equation by mode v l (r) and integrating, we obtain for the term with the pressure (taking into account the solenoidal nature of velocity and the divergence theorem) Therefore, we assume that a term with pressure is excluded. These two cases do not differ in principle, the difference is only in the technical details of the subsequent calculations, so we confine ourselves to the second case. The implementation of the Galerkin procedure will lead to the system: In this system the coefficients are calculated by the following formulas: The main problem is precisely in the computation of these integrals. It should be noted that the calculations are carried out in a spherical coordinate system, where the expressions are very cumbersome even for the basic differential operators like gradient, curl, and divergence. It is very cumbersome, even in view of the solenoidal nature of velocity modes, the invariant representation [7]: The transition from such invariant representations through the gradient and curl operators to the tensor formalism in spherical coordinates does not simplify the situation. If we take into account that the modes of speed and induction themselves have a very complicated form (3,4), it becomes clear that analytical calculations without the use of computer systems are unreal. Even taking into account the obvious symmetries A li = A il , P li = −P il , H si = H is , C li = C il , D si = D is situation varies little. The last two symmetries are a consequence of the self-adjointness of the Laplace operator.
The situation will be somewhat simplified if the normalized eigenmodes of the free dissipation of velocity and magnetic induction are used as the basis modes. Then A li = H li = δ li and matrices with elements C li and D si become diagonal. However, the most difficult to calculate are the coefficients of the nonlinear terms of the system (7).
Consider, for example, the integrand in the coefficient W si j for the case when all three modes are poloidal: It is clear that after substituting expressions for spherical functions and radial functions, we obtain a construction for which even a simple calculation of all curles and vector products is not possible manually. But we still need to take the integral.
From all of the above, it follows that there is simply no alternative to the use of systems of analytical computation for compiling models of type (7). Here it is also necessary to take into account that after selecting the modes and calculating the model it may turn out that it does not describe all the necessary processes and its interactions inherent in the dynamo-system (1,2). For example, if all the coefficients are zero, we find that the model (7) does not contain information about the rotation of the shell. In these case, we have to change the set of modes and, at least in part, recalculate the coefficients.

The use of symbolic computations for the compilation of model equations
In the Introduction, we mentioned that it is convenient to compile models (7), if the system of symbolic computations used contains a library of functions of vector analysis, and supports work in spherical coordinates. But even if there is no such library, it's easy to describe a minimal set of user-defined structures and functions, sufficient to solve the described problems. For convenience, we will use names of data types and functions from the VectorCalculus library of the Maple package [3,4]. First of all, we need to define the data structure VectorField, representated of vector fields. This structure must obviously be described as a structural type with three components corresponding to three projections of the field onto the vectors of the local spherical basis. The type of each component is a character expression. We deliberately avoid using the term «field» to refer to components of structural type, so that there is no confusion with scalar and vector fields. If the system used supports an object-oriented paradigm, then VectorField is defined as a class. Each vector field will be a variable of the VectorField type or an object of the VectorField class.
Next, we describe functions or methods that implement a set of necessary operations: 1. addition and subtraction of vector fields (using the object approach it is better to overload the addition and subtraction operations for class objects); 2. product of a vector field and a scalar field -ScalarVectorProduct; 3. scalar product of two vector fields -DotProduct; 4. vector product of two vector fields -CrossProduct; 5. divergence of vector field -Divergence; 6. curl of vector field -Curl; 7. gradient of scalar field -Gradient.
Formulas for the symbol representation of these operations (formulas of vector analysis in spherical coordinates) can be found in many sources, for example, in the popular handbook [7].
Then we need to specify analytic expressions for spherical functions, unless, of course, they are standard in the system used. To determine them we use known recurrence relations for Legendre polynomials and associated Legendre functions [7]. We can simply define the spherical harmonics Y m n (θ, ϕ) in the system as a function of arguments (n, m, θ, ϕ) by recursion over integer variables n and m, but then each of its calls will take a long time. Therefore, it is better to simply define a twodimensional array Y[n,m] of symbolic expressions with variables (θ, ϕ), having previously determined the largest expected value of index n.
Similarly, having determined the maximum assumed value of the radial index k, we form the expression arrays for the radial functions u T knm , b T knm , u P knm , b P knm . Then, we need to determine the specific type of basis modes. For example, if the modes of free dissipation are used for them, then all radial functions do not depend on index m, and are some linear combinations of power functions and spherical Bessel functions [8]. Expressions for such functions are also easily obtained from known recurrence relations [7].
So, let the expressions for all radial functions and spherical harmonics be constructed and written into the corresponding arrays.
We now define two arrays of multi-indexes for the velocity and magnetic field modes. Their elements are the main input data of the described program of symbolic computations. Based on the values of the elements of these two arrays, we compose arrays of VectorField type, calculating their elements by the values of the multi-index and formulas (3, 4), using certain or standard vector analysis operations. But, at the same time, it is still better not to use specific expressions for radial functions, simply defining them as arbitrary functions of the argument r. This will make subsequent calculations more efficient.
After filling the array of fields, we can proceed directly to the calculation of the arrays of coefficients (8) in cycles.
At each step of the cycle an integrand is constructed, which is multiplied by the product of the scale factors of the spherical coordinates. This expression is first analytically integrated over ϕ in the range from −π to π. If the integral is zero, then the corresponding coefficient is assigned a value of zero. Otherwise, the obtained expression is integrated over θ in the range from 0 to π. If the integral is zero, then the corresponding coefficient is assigned a value of zero.
With this calculation technique, coefficients that are exactly zeroed out at the stage of integration over the sphere are determined. And specific expressions for radial functions are not needed. The system only needs to know, that these are some functions of the argument r. Due to this, the system does not waste computer time and RAM for processing expressions for radial functions. And only if the integral over the sphere turns out to be nonzero, it is necessary to substitute specific expressions for radial function in place of abstract functions and carry out integration over the r. And this last integration can already be performed numerically, because, as a rule, expressions for radial functions contain parameters that are determined numerically. In such a situation, analytical integration is inexpedient.
An exception is the work in the sphere, since at the zero lower limit part of the integrals along the radius become improper. In this case, it is necessary to first set the expressions for the radial functions without specifying the numerical coefficients, to try to perform analytical integration, and, if successful, substitute the values of the coefficients already in the final expression.

Conclusion
In this paper, the technology of using the systems of symbolic computation in some problems of geophysical hydrodynamics is described. It was implemented by the authors in the Maple12 package and was successfully used to develop and study some models of geodynamo, for example [8][9][10][11]. The technology allows them to automate the routine operations for calculating the coefficients of Galerkin approximations in hydrodynamic problems and significantly increases the reliability and speed of computations.