Recently, Jouve et al(A&A, 2008) published the paper that presents the numerical benchmark for the solar dynamo models. Here, I would like to show a way how to get it with help of computer algebra system Maxima. This way was used in our paper (Pipin & Seehafer, A&A 2008, in print) to test some new ideas in the large-scale stellar dynamos. In the present paper I complement the dynamo benchmark with the standard test that address the problem of the free-decay modes in the sphere which is submerged in vacuum.
Deep Dive into Benchmarking the solar dynamo with Maxima.
Recently, Jouve et al(A&A, 2008) published the paper that presents the numerical benchmark for the solar dynamo models. Here, I would like to show a way how to get it with help of computer algebra system Maxima. This way was used in our paper (Pipin & Seehafer, A&A 2008, in print) to test some new ideas in the large-scale stellar dynamos. In the present paper I complement the dynamo benchmark with the standard test that address the problem of the free-decay modes in the sphere which is submerged in vacuum.
In equations above, the turbulent contribution is expressed through the components of the mean electromotive force (MEMF) E = u × b , where u , b are the small-scale fluctuated velocity and magnetic field respectively, Ω = Ω (r, θ) -the given angular velocity distribution.
For the sake of simplicity we restrict consideration to the case of αΩ dynamo with isotropic turbulent diffusion. Hence, we have
where η T -turbulent diffusion, α (r, θ) -dimensionless function to model the α effect, G = ∇ log ρ -stratification parameter. We adopt the model parameters given in [2]. The boundary conditions are ∂rB ∂r = 0, A = 0 -at the bottom and vacuum conditions -at the top of convection zone. For the computation all the equations are written in dimensionless form with new radial coordinate x = r/R . Moreover, to project equations on the basis function we use the coordinate transformation to the interval where basis is orthogonal.
As the first step we consider the solutions for the largest free-decay modes. This intends to test the accuracy of the boundary conditions implementation procedure and the speed of convergence. We neglect all the dynamo effects in (1,2) and return to equations:
∂ ∂θ
∂ ∂θ
where for the sake of simplicity we have assumed that magnetic diffusivity is constant over the depths and equations are written in dimensionless form. Lets consider the integration domain on the radial coordinate to be x ∈ [0, 1] and for latitude -from pole to pole. Having the regular conditions at the origin and vacuum boundary conditions at the top, we can get the analytical solutions of (6, 7). They can be expressed via the spherical Bessel functions. For the sake of simplicity we restrict ourselves to solutions for the largest modes. We decompose the eigenmodes as follows
A (x, θ, t) = e λt n m ânm sin (θ) S (A) nm (x) P 1 m (cos (θ)) ,
where P 1 m are the associated Legendre polynomials and other functions were found via basis recombination of the Legendre polynomials, namely,
Here, each element of basis satisfies the boundary conditions individually. For the largest decay modes we have
Similar to [3], the eigenvectors are scaled as S (B) 1
The results of benchmark are given in the Table 1.
Next benchmark is related to one given in the paper [2]. The detail of the model can be found in that paper. We consider the results for test B, which is for αΩ dynamo with external vacuum boundary conditions and jump othe magnetic diffusivity below the botom of convection zone. For the magnetic fields we have used the following set of the modes:
3.83e-5 6.849e-4 4 9.984e-8 4.365e-9 5 1.207e-10 2.497e-12 6 8.707e-14 9.068e-16 7 1.77e-14 4.80e-19 8 5.32e-15 6.263e-23
3.98e-9 5.651e-11 1.255e-12 9.68e-14 2.142e-16 7.72e-14 2.136e-20 1.38e-14 1.325e-24 1.90e-14 4.427e-29 where factor f = 2 1 -x i appears due the coordinate transformation [x i , 1] → [-1, 1] and
x i = 0.65 is bottom of integration domain.
The results are shown at the table 2 and Figure 1. Note, that our results are quite close to those in [2]. The small difference can be attributed to the difference in the method of solution of the problem.
Here, I give some details of the above computation within Maxima. The strongest side of the computer algebra system is that most of the computational work is done exactly. For example, all the derivatives in (6, 7) are computed directly by substitution of (8, 9) to (6, 7). The integration over space domain has to be done to proceed with the Galerkin method. This is carried out with help of the Gauss-Legendre procedure. The procedure is based on the so-called “collocation points and weights” method. The collocation points are the set of zeros of the Legendre polynomial of the order > N where N is the number of the greatest mode among n = 1, . . . N in (8, 9). The collocation points and weights can be found via subprogram “pseudp.mac”. It is given in Appendix. The maxima session is started with the general definitions: kill(all)$ batchload("pseudp.mac")$ (xb:0.,xe:1.)$ bftorat:true$ float2bf:true$ ratprint:false$ fpprec:64$ ratepsilon:1.e-32$ Note, that only a half of [-1, 1] intervals is used in computations. This leads to the computational economy and to increase the spectral resolution of the code. The following are the part of the code which defines the basis functions on radius (S As you see, we use odd or even associated Legendre polynomial in respective of the parity choice. To accelerate the integration procedure we calculate the matrices of basis functions over the sets of collocation points, -328.0254017082749, -373.0041226108443, -410.8018157235754, -505.412713312208, -577.4557504236772, -665.4381736491396, -809.3382250971157, -851.3024977594987, -882.9750962899574, -932.6139588501989, -970.2769707279747, -1049.716092377411, -1259.650734653339, -1830.07144488619, -2585.672953057317] The relevant spherical Bessel functions of the problem are jn(n,x):=spherical_bessel_j(n,x)$ jnR(n,x):=x*spherical_bessel_j(n,x)$
…(Full text truncated)…
This content is AI-processed based on ArXiv data.
