Control namelist

This namelist defines the numerical parameters of the problem plus the variables that control and organize the run.

• mode (default mode=0) is an integer which controls the type of calculation performed.

  mode=0	Self-consistent dynamo
  mode=1	Convection
  mode=2	Kinematic dynamo
  mode=3	Magnetic decay modes
  mode=4	Magneto convection
  mode=5	Linear onset of convection
  mode=6	Self-consistent dynamo, but with no Lorentz force
  mode=7	Super-rotating inner core or mantle, no convection and no magnetic field
  mode=8	Super-rotating inner core or mantle, no convection
  mode=9	Super-rotating inner core or mantle, no convection and no Lorentz force
  mode=10	Super-rotating inner core or mantle, no convection, no magnetic field, no Lorentz force and no advection

• tag (default tag="default") is a character string, used as an extension for all output files.

• n_time_steps (default n_time_steps=100) is an integer, the number of time steps to be performed.

• tEND (default tEND=0.0) is a real, which can be used to force the code to stop when :math:t=tEND. This is only used when t/=tEND.

• alpha (default alpha=0.5) is a real. This is the weight used for current time step in implicit time step.

Default scales

• n_tScale (default n_tScale=0) is an integer, which determines the time scaling

  n_tScale=0	Use viscous time scale.	\(d^2/\nu\)
  n_tScale=1	Use magnetic time scale.	\(d^2/\eta\)
  n_tScale=2	Use thermal time scale.	\(d^2/\kappa\)
  n_tScale=3	Use rotational time scale.	\(\Omega^{-1}\)

• n_lScale (default n_lScale=0) is an integer which determines the reference length scale.

  n_lScale=0	Use outer core.
  n_lScale=1	Use total core.

• enscale (default enscale=1.0) is a real. This is the scaling for energies.

Update control

• l_update_v (default l_update_v=.true.) is a logical that specifies whether the velocity field should be time-stepped or not.

• l_update_b (default l_update_b=.true.) is a logical that specifies whether the magnetic field should be time-stepped or not.

• l_update_s (default l_update_s=.true.) is a logical that specifies whether the entropy/temperature should be time-stepped or not.

• l_update_xi (default l_update_xi=.true.) is a logical that specifies whether the chemical composition should be time-stepped or not.

• l_update_phi (default l_update_phi=.true.) is a logical that specifies whether the phase field should be time-stepped or not.

Time step control

A modified Courant criterion including a modified Alfven-velocity is used to account for the magnetic field. The relative and absolute importance of flow and Alfven-velocity can be controled by courfac and alffac respectively. The parameter l_cour_alf_damp allows to choose whether the actual Alven speed is used to estimate the Courant condition or if damping is included. Practically, the timestep size is controlled as follows

\[\delta t < \min_{V}\left( c_I\,E,\, \dfrac{\delta r}{|u_r|},\, \dfrac{\delta h}{u_h} \right)\]

where \(u_h=(u_\theta^2+u_\phi^2)^{1/2}\), \(\delta h = \dfrac{r}{\sqrt{\ell(\ell+1)}}\), and \(\delta r\) is the radial grid interval. The first term in the left hand side accounts for the explicit treatment of the Coriolis term.

\[{|u_r|}=c_F{|u_{F,r}|}+c_A\dfrac{u_{A,r}^2}{\left[u_{A,r}^2+\left(\frac{1+Pm^{-1}}{2\delta r}\right)^2\right]^{1/2}}\,,\]

where \(u_{F,r}\) is the radial component of the fluid velocity and \(u_{A,r}=Br/\sqrt{E\,Pm}\) is the radial Alven velocity. The denominator of the rightmost term accounts for the damping of the Alven waves.

• dtMax (default dtMax=1e-4) is a real. This is the maximum allowed time step \(\delta t\). If \(\delta t > \hbox{dtmax}\), the time step is decreased to at least dtMax (See routine dt_courant). Run is stopped if \(\delta t < \hbox{dtmin}\) and \(\hbox{dtmin}=10^{-6}\,\hbox{dtmax}\).

• courfac (default courfac=2.5) is a real used to scale velocity in Courant criteria. This parameter corresponds to \(c_F\) in the above equation.

• alffac (default alffac=1.0) is a real, used to scale Alfven-velocity in Courant criteria. This parameter corresponds to \(c_A\) in the above equation.

• intfac (default intfac=0.15) is a real, used to scale Coriolis factor in Courant criteria. This parameter corresponds to \(c_I\) in the above equation.

• l_cour_alf_damp (default l_cour_alf_damp=.true.) is a logical. This is used to decide whether the damping of the Alven waves is taken into account when estimating the Courant condition (see Christensen et al., GJI, 1999). At low Ekman numbers, this criterion might actually lead to spurious oscillations/instabilities of the code. When turn to False, \({|u_r|}=c_F{|u_{F,r}|}+c_A{|u_{A,r}|}\).

• time_scheme (default time_scheme='CNAB2') is a character string. This is used to choose the time step integrator used in the code among the following implicit-explicit time schemes:

  time_scheme=’CNAB2’	Crank-Nicolson and 2nd order Adams-Bashforth scheme
  time_scheme=’CNLF’	Crank-Nicolson and Leap-Frog scheme
  time_scheme=’MODCNAB’	Modified CN/AB2
  time_scheme=’SBDF2’	Semi-implicit backward difference scheme of 2nd order
  time_scheme=’SBDF3’	Semi-implicit backward difference scheme of 3rd order
  time_scheme=’SBDF4’	Semi-implicit backward difference scheme of 4th order
  time_scheme=’ARS222’	Semi-implicit S-DIRK of 2nd order
  time_scheme=’ARS232’	Semi-implicit S-DIRK of 2nd order
  time_scheme=’CK232’	Semi-implicit S-DIRK of 2nd order
  time_scheme=’LZ232’	Semi-implicit S-DIRK of 2nd order
  time_scheme=’PC2’	Semi-implicit S-DIRK of 2nd order
  time_scheme=’CB3’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’ARS343’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’MARS343’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’ARS443’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’BPR353’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’BHR553’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’DBM453’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’LZ453’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’KC343’	Semi-implicit S-DIRK of 3rd order
  time_scheme=’KC564’	Semi-implicit S-DIRK of 4th order
  time_scheme=’KC674’	Semi-implicit S-DIRK of 4th order
  time_scheme=’KC785’	Semi-implicit S-DIRK of 5th order

Run time

The total desired runtime (in human units and not in CPU units) can be specified with the three variables runHours, runMinutes and runSeconds.

• runHours (default runHours=0) is an integer that controls the number of run hours.

• runMinutes (default runMinutes=0) is an integer that controls the .

• runSeconds (default runSeconds=0) is an integer that controls the number of run hours.

Here is an example for a run of 23h30:

runHours   = 23,
runMinutes = 30,

Hyperdiffusivity

Hyperdiffusion can be applied by multiplying the diffusion operators by a factor of the form

\[d(\ell)=1+D\left[\frac{\ell+1-\ell_{hd}}{\ell_{max}+1-\ell_{hd}} \right]^{\beta}\]

for the spherical harmonic degrees \(\ell \geq \ell_{hd}\).

• difnu (default difnu=0.0) is a real. This is the amplitude \(D\) of the viscous hyperdiffusion.

• difkappa (default difkappa=0.0) is a real. This is the amplitude \(D\) of the thermal hyperdiffusion.

• difchem (default difchem=0.0) is a real. This is the amplitude \(D\) of the hyperdiffusion applied to chemical composition.

• difeta (default difeta=0.0) is a real. This is the amplitude \(D\) of the magnetic hyperdiffusion.

• ldif (default ldif=1) is an integer. This is the degree \(\ell_{hd}\) where hyperdiffusion starts to act.

• ldifexp (default ldifexp=-1) is an integer. This is the exponent \(\beta\) of hyperdiffusion.

Angular momentum correction

In case of the use of stress-free boundary conditions at both boundaries, it is safer to ensure that the angular momentum is correctly conserved. This can be enforced through the following input variables:

• l_correct_AMe (default l_correct_AMe=.false.) is a logical. This is used to correct the equatorial angular momentum.

• l_correct_AMz (default l_correct_AMz=.false.) is a logical. This is used to correct the axial angular momentum.

Radial scheme and mapping of the Gauss-Lobatto grid

In MagIC, one can either use finite differences or Chebyshev polynomials for the radial integration scheme. This choice is controlled by the following input parameter:

• radial_scheme (default radial_scheme='CHEB') is a character string.

  radial_scheme=’CHEB’	Use Chebyshev polynomials
  radial_scheme=’FD’	Use finite differences

When Chebyshev polynomials are used, it is also possible to use a non-linear mapping function to concentrate/diperse grid points around a point inside the domain.

• l_newmap (default l_newmap=.false.) is a logical. A radial mapping can be applied to the Chebyshev grid when l_newmap is set to .true.. The radial profile of the mapping function is then stored during the initialisation of the code in the file rNM.TAG.

• map_function (default map_function='arcsin') is a character string. This allows to select which mapping function is used:

  map_function=’TAN’	Use a tangent mapping (see Bayliss and Turkel 1992)
  map_function=’ARCSIN’	Use an arcsin mapping (see Kosloff and Tal-Ezer 1993)
  map_function=’TT’	Use the mapping by Tee and Trefethen 2006
  map_function=’JAFARI’	Use the mapping by Jafari-Varzaneh and Hosseini 2014

If the tangent mapping is used, the function that re-distributes the collocation points is expressed by

\[r=\frac{1}{2}\left(\alpha_2+\frac{\textrm{tan}\left[\lambda(x_{cheb}-x_0)\right]}{\alpha_1}\right) + \frac{r_i+r_o}{2} \textrm{ ,}\]

where the Gauss-Lobatto collocation points are

\[x_{cheb}=\textrm{cos}\left( \frac{\pi(k-1)}{N_r} \right) \textrm{ , }\;\; k=1,2,...,n_r \textrm{ , }\; n_r=n\_r\_max\]

and \(r\!\in\![r_i,r_o]\), \(x_{cheb}\!\in\![-1.0,1.0]\). The parameters to calculate \(r\) are

\[\begin{split}\lambda&=\frac{\textrm{tan}^{-1}\left(\alpha_1(1-\alpha_2)\right)}{1-x_0} \\ x_0&=\frac{K-1}{K+1} \\ K&=\frac{\textrm{tan}^{-1}\left(\alpha_1(1+\alpha_2)\right)}{\textrm{tan}^{-1}\left(\alpha_1(1-\alpha_2)\right)} \textrm{ .}\end{split}\]

The coefficient \(\alpha_1\) determines the degree of concentration/dispersion of the grid points around \(x_{cheb}\!=\!\alpha_2\). If \(\alpha_1\) is too high, the \(r\) function becomes nearly discontinuous. To avoid numerical problems, \(\alpha_1\) should remain close to unity.

If the arcsin mapping is used, the function that re-distributes the collocation points is given by

\[r=\frac{1}{2}\left[ \frac{\textrm{arcin}\left(\alpha_1 x_{cheb}\right)}{\textrm{arcsin} \alpha_1} \right]+\frac{r_i+r_o}{2} \textrm{ ,}\]

In the Kosloff and Tal-Ezer mapping, \(\alpha_1\) transforms the Gauss-Lobatto grid into a more regularly-spaced grid. When \(\alpha_1 \rightarrow 0\) one recovers the Gauss-Lobatto grid, while \(\alpha_1 \rightarrow 1\) yields a regular grid.

Warning

The Kosloff-Tal-Ezer mapping becomes singular when \(\alpha_1=1\). Acceptable values are \(0<\alpha_1<1\). Note that the error increases as \(\epsilon=\left(\frac{1-\sqrt{1-\alpha_1^2}}{\alpha_1}\right)^{N_r}\).

If the Tee and Trefethen sinh mapping is employed, the grid points are redistributed in the following manner

\[r=\frac{1}{2}\left(\alpha_2+\frac{\textrm{sinh}\left[A(x_{cheb}-1)+B\right]}{\alpha_1}\right) + \frac{r_i+r_o}{2} \textrm{ ,}\]

where

\[A=\frac{1}{2}\left[\textrm{sinh}(\alpha_1(1-\alpha_2))+\textrm{sinh}(\alpha_1(1+\alpha_2)) \right], \quad B = \textrm{sinh}(\alpha_1(1-\alpha_2))\]

With this mapping, \(\alpha_1\) is directly related to the stiffness of the transition.

If the Jafari-Varzaneh and Hosseini mapping is employed, similarly to the tangent mapping, \(\alpha_1\) determines the degree of concentration of the grid points around \(x_{cheb}\!=\!\alpha_2\). This is expected to do a better job than the tangent mapping, both in terms of matrix conditioning and in terms of reducing the Gibbs phenomenon around a steep change (Allen-Cahn type of equations involved in the phase field model comes to mind).

• alph1 (default alph1=0.8) is a real. This is a control parameter of the mapping function.

• alph2 (default alph2=0.0) is a real. This is a control parameter of the mapping function. The default value of \(0\) corresponds to the center of the grid.

Miscellaneous

• l_non_rot (default l_non_rot=.false.) is a logical. Use it when you want to do non-rotating numerical simulations.

• anelastic_flavour (default anelastic_flavour="None") is a character string. This allows to change the thermal diffusion operator used within the anelastic approximation. Possible values are:

  anelastic_flavour=’LBR’	Entropy diffusion
  anelastic_flavour=’ENT’	Entropy diffusion
  anelastic_flavour=’ALA’	Anelastic liquid approximation
  anelastic_flavour=’TDIFF’	Temperature diffusion
  anelastic_flavour=’TEMP’	Temperature diffusion

• polo_flow_eq (default polo_flow_eq="WP") is a character string. This allows to change how the equation for the poloidal flow potential is constructed. One can either use the radial component of the Navier-Stokes equation and hence keep a coupled system that involve the poloidal potential \(W\) and the pressure \(p\), or take the radial component of the double-curl of the Navier-Stokes equation to suppress pressure.

  polo_flow_eq=’WP’	Use the pressure formulation
  polo_flow_eq=’DC’	Use the double-curl formulation

• mpi_transp (default mpi_transp="auto") is a character string. It allows to change the way the global MPI transposes are handled by the code. By default, the code tries to determine by itself the fastest method. One can nevertheless force the code to use local communicators (such as Isend/Irecv/waitall), make use of the native alltoallv MPI variant or choose the alltoallw variant instead.

  mpi_transp=’auto’	Automatic determination of the fastest transpose
  mpi_transp=’p2p’	Use Isend/Irecv/Waitall communicators
  mpi_transp=’a2av’	Use alltoallv communicators
  mpi_transp=’a2aw’	Use alltoallw communicators

• mpi_packing (default mpi_packing="packed") is a character string. It allows to change the size of the global MPI transposes. One can choose between some packing of the fields into buffers (default) or a sequence of single field transposes. There is a possible automatic detection but testing unfortunately reveals frequent false detection.

  mpi_packing=’auto’	Automatic determination of the fastest transpose
  mpi_packing=’packed’	Pack some fields into buffers
  mpi_packing=’single’	Transpose each field individually

• l_adv_curl (default l_adv_curl=.true.) is a logical. When set to True, the advection term is treated as \(\vec{u}\times\vec{\omega}\) instead of \(\vec{u}\vec{\nabla}\vec{u}\). The practical consequence of that is to reduce the number of spectral/spatial Spherical Harmonic Transforms and hence to speed-up the code. Because of the treatment of the viscous heating term in the anelastic approximation, this is only an option when considering Boussinesq models.