Chapter 2 Setting up the simulation

XSHELLS can time-step the Navier-Stokes equation in a rotating reference frame. Optionally it can include (i) a buoyancy force in the Boussinesq approximation, where the buoyancy has one or two sources obeying distinct advection-diffusion equations; and (ii) a Lorentz (or Laplace) force for conducting fluids where the magnetic field obeys the induction equation. Precisely, the following equations can be time-stepped by XSHELLS:

	${\partial }_t𝐮+(2{𝛀}_{\mathrm{𝟎}}+\nabla \times 𝐮)\times 𝐮$	$=-\nabla p^{*}+\nu {\nabla }^2𝐮+(\nabla \times 𝐛)\times 𝐛+(c+T)\nabla {\mathrm{\Phi }}_0$		(2.1)
	${\partial }_t𝐛$	$=\nabla \times (𝐮\times 𝐛-\eta \nabla \times 𝐛)$		(2.2)
	${\partial }_{t}T+𝐮.\nabla (T+T_0)$	$=\kappa {\nabla }^{2}T$		(2.3)
	${\partial }_{t}c+𝐮.\nabla (c+C_0)$	$={\kappa }_c{\nabla }^{2}c$		(2.4)
	$\nabla 𝐮$	$=0$		(2.5)
	$\nabla 𝐛$	$=0$		(2.6)

where

• •

  $𝐮$ is the velocity field.

• •

  $T$ and $c$ are respectively the temperature and concentration fields, that contribute to the buoyancy in the Boussinesq formulation.

• •

  $\sqrt{{\mu }_0\rho }𝐛$ is the magnetic field ($\rho$ and ${\mu }_0$ are the fluid density and magnetic permeability, but are not relevant for this equation set).

• •

  $\nu$ is the kinematic viscosity of the fluid and is set by the variable nu in xshells.par.

• •

  $\eta ={({\mu }_0\sigma )}^{-1}$ is the magnetic diffusivity of the fluid ($\sigma$ is its conductivity) and is set by the variable eta in xshells.par. Diffusivity $\eta (r)$ depending on shell radius are also supported (see sec. 2.2.3).

• •

  $\kappa$ and ${\kappa }_c$ are respectively the thermal and chemical diffusivities of the fluid and are set by the variable kappa and kappa_c in xshells.par.

• •

  ${𝛀}_{\mathrm{𝟎}}$ is the rotation vector of the reference frame, which is usually along the vertical axis ${𝐞}_{𝐳}$. It is set by the variable Omega0 in xshells.par, while the angle (in radians) between ${𝐞}_{𝐳}$ and ${𝛀}_{\mathrm{𝟎}}$ is set by Omega0_angle (0 by default). Note that ${𝛀}_{\mathrm{𝟎}}$ is always in the $x-z$ plane ($\varphi =0$).

• •

  ${\mathrm{\Phi }}_0$ is the gravity potential (independent of time), controlled by field phi0 in xshells.par.

• •

  $T_0$ and $C_0$ are the imposed base temperature and concentration profiles, controlled by fields tp0 and c0 in xshells.par.

• •

  $p^{*}$ is the dynamic pressure deviation (from hydrostatic equilibrium), which is eliminated by taking the curl of equation 2.1.

Equation 2.1 (respectively 2.2, 2.3 and 2.4) is time-stepped when u (respectively b, tp and c) is set to an initial condition in xshells.par. Disabling an equation is as easy as removing or commenting out the corresponding initial condition in xshells.par. Note that currently, the concentration equation cannot be time-stepped without the temperature equation.

Example If the following lines are found in xshells.par: nu = 1.0 eta = sqrt(10) Omega0 = 2*pi*1e3 u = 0 b = 0 #tp = 0 then the viscosity is set to $\nu =1.0$, the magnetic diffusivity is set to $\eta =\sqrt{10}$, and the rotation rate is set to ${\mathrm{\Omega }}_0=2\pi \times {10}^3$. The Navier-Stokes (2.1) and the induction equation (2.2) will be time-stepped, but not the temperature and concentration equations (2.3), simulating an isothermal fluid.

Note that it is up to the user to choose dimensional or non-dimensional control parameters. A notable exception is the magnetic field, which is always scaled to the same unit as the velocity field.

Besides thermal convection, mechanical forcing can be imposed at the boundaries.

Predefined variable a_forcing and w_forcing define the amplitude and frequency of a forcing. The precise nature of the forcing (e.g. differential rotation) must be defined in the xshells.hpp file before compilation (see section 2.2.6).

XSHELLS uses semi-implicit (or IMEX) time steppers, where the diffusive terms are treated implicitly, while all other terms (including non-linear terms) are treated explicitly. Since v2.4, XSHELLS offers a selection of time steppers:

• •

  PC2 (the default) is a second order predictor-corrector scheme. Although this scheme requires three times more work per time-step than the classical Crank-Nicolson-Adams-Bashforth (CNAB) scheme, it allows time-steps from three to four times that of CNAB, leading to shorter time-to-solution¹¹ 1 before version 2.0, a classical CNAB scheme was used.

• •

  CNAB2 is the classical second order Crank-Nicolson-Adams-Bashforth (CNAB) scheme²² 2 CNAB2: see eq. 2.9 of Wang & Ruuth (2008) https://www.jstor.org/stable/43693484..

• •

  SBDF2 is the semi-implicit Backward Difference Formula of order 2. This scheme³³ 3 SBDF2: see eq. 2.8 of Wang & Ruuth (2008) https://www.jstor.org/stable/43693484. has a stringent time-step restriction with the Coriolis force.

• •

  SBDF3 is the semi-implicit Backward Difference Formula of third order⁴⁴ 4 SBDF3: see eq. 2.14 of Wang & Ruuth (2008) https://www.jstor.org/stable/43693484..

• •

  BPR353 is a third order Runge-Kutta scheme featuring good stability and high accuracy⁵⁵ 5 BPR353: see page A49 of Boscarino, Pareschi, Russo (2013) https://doi.org/10.1137/110842855 https://arxiv.org/abs/1110.4375..

Simply add stepper = SBDF3 in the xshells.par file to use the SBDF3 scheme. Note that the allowable time steps can vary largely between these schemes. The automatic time-step selection (see below) should handle most problems well.

Both PC2 and CNAB2 use the Crank-Nicolson scheme for the implicit part, and as a consequence are not very accurate when computing decay rates (e.g. magnetic dipole decay rate). On the other hand SBDF2 and SBDF3 are accurate for decay rates, but non-linear terms like advection impose small time-steps for the scheme to remain stable.

The time-step of the numerical integration is set by dt in xshells.par.

The sub_iter variable is half the number of time-steps taken before any diagnostic is computed and written to file energy.job or displayed on screen. For example, if sub_iter = 50, then 100 time-steps will be performed before computing and printing some diagnostics. This is then called an iteration.

The iter_max variable is the total number of iterations, so that the total number of time-steps before the code will stop is iter_max $\times \mathrm{\hspace{0.17em}2}\times$ sub_iter.

By setting dt_adjust = 1, an automatic time-step adjustment can be turned on. The time-step adjustment can be controlled through advanced options (see §2.1.10). In that case, the number of sub-iterations sub_iter is also adjusted so that an iteration is a constant time span, and thus the outputs happen at fixed time intervals $\mathrm{\Delta }T=2\times$ sub_iter $\times$ dt.

Finally, iter_save controls the number of iterations before a (partial) snapshot is saved to disk.

Example The following lines in xshells.par will use a time-step of $0.01$ for the numerical integration: dt_adjust = 0 # 0: dt fixed (default), 1: variable time-step dt = 0.01 # time step iter_max = 300 # iteration number (total number of text and # energy file ouputs) sub_iter = 25 # sub-iterations (the time between outputs # = 2*dt*sub_iter is fixed even with variable dt) iter_save = 10 # number of iterations between field writes Output will occur every $\mathrm{\Delta }T=0.01\times 25\times 2=0.5$ time units (an iteration). The program will stop after iter_max=300 outputs (or iterations), spanning a total physical time of $t_{end}-t_{start}=150.0$. Partial fields are saved every 10 iterations, or every 5.0 physical time units, if movie = 1 is set (see below).

At each iteration, XSHELLS can plot the kinetic and magnetic energies as a function of time, using gnuplot. Note that the plots are refreshed every iteration, but no more than once every two seconds. This allows to follow program execution in real-time, but might not be useful for high performance distributed jobs. The interaction with gnuplot must be turned on by passing the --enable-gnuplot option to ./configure.

The variable plot in the file xshells.par then allows some control:

• •

  plot = 0: disables plotting.

• •

  plot = 1: shows plot on display; if no display found, write to png file instead. This is the default.

• •

  plot = 2: saves plot to png file only.

• •

  plot = 3: shows plot on display (if available) and also saves plot to png file.

The parameter movie controls the field snapshots, saved every iter_save iterations (see above).

• •

  movie = 1 : the initial field is saved to fieldX_0000.job, after iter_save iterations the fields are saved to fieldX_0001.job, then fieldX_0002.job, and so on.

• •

  movie = 0 : no such fields are saved.

• •

  movie = 2 : in addition to the snapshots of the fields, the time-average since the last snapshot is also computed and saved.

The parameter prec_out controls the precision (single or double precision) of the snapshot files. The snapshots are saved in double precision by default (prec_out = 2), which allows restarting or computing gradients. If you need to save disk space, you can use single precision instead by setting set prec_out = 1, thus writing field files of half size. To save further disk space, snapshots can be truncated at lower spherical harmonic degree and order, using the parameters lmax_out and mmax_out, respectively.

Sometimes, single precision is not enough but double precision is not needed. Setting set prec_out = 3 uses the custom fp48 format which takes 25% less space than double precision, while keeping high accuracy.

The snapshots can then be post-processed with xspp to produce plots or movies (see 3). Note that field in fp48 format are read seamlessly by xspp which can also convert between formats. However, pyxshells cannot read the fp48 format.

Example The following lines in xshells.par instruct the program to output snapshots and time-averages of the axisymmetric component of the fields, every iter_save iterations: movie = 2 # 0=field output at the end only (default), # 1=output every iter_save, 2=also writes # time-averaged fields lmax_out = -1 # lmax for movie output (-1 = same as Lmax, which # is also the default) mmax_out = 0 # mmax for movie output (-1 = same as Mmax, which # is also the default) prec_out = 2 # write double precision snapshots.

By default after initialization the job starts at the beginning (iteration 0). It is easy to start a new job by using as input fields the field files written by a previous job, effectively continuing that job.

Sometimes, it is useful to automatically continue a stopped or killed job (e.g. in batch execution environments found in high-performance computing machines). By default, a full resolution snapshot is written to disk every four hours. Parameters found in xshells.par allow to tune that interval, and enable restart from these checkpoint automatically when the program is run again.

For increased safety, when writing a new checkpoint (or backup) to file fieldX.jobname, the previous one is first renamed to fieldX_back.jobname. This file may allow to continue a simulation in case of an unexpected termination of the program while writing the new checkpoint.

If restart = 1, the program will start by looking in the current directory for checkpoint files that have been saved by a previous run with the same job name, and use these to resume that job.

Example Suppose that on a supercomputer, the wall time of the jobs is limited to 24 hours. In order to run a job that spans several days, the following lines in xshells.par allow a job to be resumed by simply resubmitting it: restart = 1 # 1: try to restart from a previous run with # same name. 0: no auto-restart (default). backup_time = 470 # ensures that full fields are saved to disk at # least every backup_time minutes, for restart. nbackup = 3 # number of backups before terminating program # (useful for time-limited jobs). # 0 = no limit (default) In addition, a checkpoint (or backup) is written to disk every 470 minutes, and the program will stop after writing the third backup, thus leaving 30 minutes of safety time for program initialization and file writing time. In case of a system failure, no more than 470 minutes of computing time will be lost.

The following features have not been thoroughly tested and may not work flawlessly in all situations. If you want to use them, please test and report bugs.

nonlin = ugu,uxb     # include only u.grad(u) and curl(u x b)

Enable by uncommenting:

#define XS_CUSTOM_DIAGNOSTICS

In addition to the total energy of the three fields $U$, $B$ and $T$, which are saved every 2 sub_iter time steps (see section 2.1.6), custom diagnostics can be defined in the custom_diagnostic() function, found in the xshells.hpp file. They are computed every iteration and saved in energy.job after the energies. The best is to look at the existing diagnostics defined in the custom_diagnostic() function to add your own.

Enable compilation of variable time-step code by uncommenting:

#define XS_ADJUST_DT

In addition,variable time-step must also be allowed by setting dt_adjust = 1 in file xshells.par (see also section 2.1.6).

In equation 2.2, conductivity can depend on radius $r$. To define a conductivity profile $\eta (r)$, uncomment:

#define XS_ETA_PROFILE

and then define your profile in the calc_eta() function, found in the xshells.hpp file. The purpose of calc_eta() is simply to fill the array etar with values of the magnetic diffusivity for every radial shell. The program handles continuous profiles as well as discontinuities in $\eta (r)$ properly and automatically.

In order to compute in a full sphere and avoid problems near $r=0$, the spherical harmonic expansion must be truncated at low degree near $r=0$. XSHELLS can truncate the spherical harmonic expansions at a different degree for each shell, when the following line is uncommented in xshells.hpp:

#define VAR_LTR 0.5

The value of VAR_LTR ($0.5$ in the line above, which is a good choice for full-sphere computations) is used as $\alpha$ in the formula to determine the truncation degree ${\mathrm{\ell }}_{tr}$:

$${\mathrm{\ell }}_{tr}(r)=L_{max}\sqrt{\frac{r}{\alpha r_{max}}}+1$$

where $L_{max}$ is defined by Lmax in xshells.par and $r_{max}$ is the radius of the last shell. Note that ${\mathrm{\ell }}_{tr}$ cannot exceed $L_{max}$.

Note also that $\alpha$ can be overridden by the rsat_ltr variable in xshells.par.

To avoid computing small-scales that may not be dynamically relevant, the easy way is to employ hyper-diffusivity, which increase the diffusivity of small-scales. XSHELLS can be compiled with hyper-diffusivities when the following line is uncommented in xshells.hpp:

#define XS_HYPER_DIFF

The formulae used to compute the viscosity is the one proposed by Nataf & Schaeffer (2015):

(2.9)

with $q={({\nu }_{max}/{\nu }_0)}^{1/(L_{max}-{\mathrm{\ell }}_c)}$. The cutoff ${\mathrm{\ell }}_c$ and the ratio ${\nu }_{max}/{\nu }_0$ are set respectively by the hyper_diff_l0 and hyper_nu parameters in xshells.par. If hyper_diff_l0 = 0, hyper-diffusivity is disabled. Similarly, hyper_nu = 1 gives uniform viscosity, while hyper_nu = 100 leads to a viscosity at $\mathrm{\ell }=L_{max}$ that is 100 times larger than at large scales.

Note that thermal and magnetic diffusivities can also be treated similarly, using respectively the hyper_kappa and hyper_eta parameters. The cutoff parameter ${\mathrm{\ell }}_c$ is the same for all diffusivities.

Amplitude and frequency are set at runtime by a_forcing and w_forcing in the xshells.par file.

#define XS_SET_BC

To replace the equations with their linearized version (no $(\nabla \times 𝐮)\times 𝐮$, no $(\nabla \times 𝐛)\times 𝐛$, no $𝐮.\nabla c$), uncomment:

#define XS_LINEAR

Note that spherical harmonic transforms are not needed anymore if there are no base fields, leading to a much faster program. Base fields are also supported, and can be set using u0, b0 and tp0 in the xshells.par file. See also the example located in problems/taw.

To include an axisymmetric alpha effect in the induction equation, add the following in the xshells.hpp file:

#define XS_MEAN_FIELD

And provide the alpha field using the init_alpha function. See the example located in problems/kinematic_dynamos.