XSHELLS 2.5 : User Manual

Chapter 2 Setting up the simulation

Example configuration files can be found in the problems directory.

2.1 Run-time options: xshells.par

The file xshells.par is a simple text file. Each line may contain a single statement like var = expression, or a comment starting with #. A simple math parser allows to use convenient expressions like sqrt(4*pi/3).

All the following features can be set in xshells.par. There is no need to recompile if this file is changed, as it is read and interpreted at program startup.

2.1.1 Equations and controlling parameters

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:

t𝐮+(2𝛀𝟎+×𝐮)×𝐮 =-p*+ν2𝐮+(×𝐛)×𝐛+(c+T)Φ0 (2.1)
t𝐛 =×(𝐮×𝐛-η×𝐛) (2.2)
tT+𝐮.(T+T0) =κ2T (2.3)
tc+𝐮.(c+C0) =κc2c (2.4)
𝐮 =0 (2.5)
𝐛 =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.

  • μ0ρ𝐛 is the magnetic field (ρ and μ0 are the fluid density and magnetic permeability, but are not relevant for this equation set).

  • ν is the kinematic viscosity of the fluid and is set by the variable nu in xshells.par.

  • η=(μ0σ)-1 is the magnetic diffusivity of the fluid (σ is its conductivity) and is set by the variable eta in xshells.par. Diffusivity η(r) depending on shell radius are also supported (see sec. 2.2.3).

  • κ and κc are respectively the thermal and chemical diffusivities of the fluid and are set by the variable kappa and kappa_c in xshells.par.

  • 𝛀𝟎 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 𝛀𝟎 is set by Omega0_angle (0 by default). Note that 𝛀𝟎 is always in the x-z plane (ϕ=0).

  • Φ0 is the gravity potential (independent of time), controlled by field phi0 in xshells.par.

  • T0 and C0 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 ν=1.0, the magnetic diffusivity is set to η=10, and the rotation rate is set to Ω0=2π×103. 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.

MHD without Lorentz force (e.g. kinematic dynamos)

The Lorentz force can be turned off. Just add no_jxb = 1 in the xshells.par file.

Internal representation of vector fields

Vector fields are represented internally using a poloidal/toroidal decomposition:

𝐮=×(T𝐫)+××(P𝐫) (2.7)

where 𝐫 is the radial position vector, and T and P are the toroidal and poloidal scalars respectively. This decomposition ensures that the vector field 𝐮 is divergence-free.

The scalar fields T and P for each radial shell are then decomposed on the basis of spherical harmonics.

2.1.2 Boundary conditions

Magnetic field

boundary conditions are that of an electrical insulator outside the computation domain, with or without external sources of magnetic field (see section 2.1.3 for externally imposed magnetic fields).

Temperature or concentration

boundary conditions can be fixed value (defined by the T0 or C0 profile), fixed flux (defined by rT0 or rC0), or Robin (such as rT+βT=0). For Robin, the constant β (related to the Biot number) is set using Biot_Ti and Biot_To for the inner and outer boundaries respectively (or Biot_Ci and Biot_Co for the composition field).

Velocity

boundary conditions can be zero, no-slip (with arbitrary prescribed velocity at the boundary) or stress-free.

The inner and outer boundary conditions can be chosen independently. The BC_U (for velocity), BC_T (for temperature) and BC_C (for concentration) entries in xshells.par allow to select the appropriate boundary conditions.

Example The following lines in xshells.par define zero velocity and fixed temperature boundary condition at the inner boundary, and no-slip and fixed flux boundary condition at the outer boudnary.
BC_U = 0,1    # inner,outer boundary conditions
              # (0=zero velocity, 1=no-slip, 2=free-slip)
BC_T = 1,2    # 1=fixed temperature, 2=fixed flux, 6=Robin
#Biot_To = 0.5   # Set the constant for a Robin boundary condition at the outer boundary

2.1.3 Initial conditions and imposed fields

Predefined fields

Several predefined fields are defined in xshells_init.cpp. The command ./list_fields prints a list of these predefined fields, with their name in the first column. You can simply use this name in the xshells.par file to define an initial condition. You can also add your own by editing xshells_init.cpp.

Imposed fields are only supported for the gravity potential phi0 and for the basic state of temperature tp0 and concentration c0. Imposed magnetic fields can be obtained through the appropriate boundary conditions (magnetic fields generated by currents outside the computation domain only). Some predefined magnetic field include these boundary conditions, making them actually imposed fields (and are labeled as such). Note also that a linear mode exist which support arbitrary base fields (see §2.2.7).

Example The following lines in xshells.par set up the geodyanmo benchmark initial conditions.
E = 1e-3
Pm = 5
Ra = 100

u = 0                 # initial velocity field
b = bench2001*5/sqrt(Pm*E)  # initial magnetic field (scaled by
                            #   sqrt(1/(Pm*E)))
tp = bench2001*0.1    # initial temperature field
tp0 = delta*-1        # imposed (base) temperature field
phi0 = radial*Ra/E    # radial gravity field (multiplied by
                      #  Ra/E to match geodynamo benchmark)

Field files as initial conditions

In addition, any field file can be given as initial condition. If the radial grid is not the same, the field must be interpolated on the new grid. To avoid mistakes, interpolation is disabled by default and must be enabled by interp = 1 (often found near the end of the xshells.par file).

Example The following lines start from the velocity field saved in file fieldU.previous_job, which was performed at different parameters and with a different number of radial grid points.
u = fieldU.previous_job   # initial velocity field
interp = 1     # allow interpolation, to be able to use fields
       # defined on a different radial grid as initial condition.

2.1.4 Forcing

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).

2.1.5 Spatial discretization

Radial grid

XSHELLS uses second order finite differences in radius. The total number of radial grid points is defined in xshells.par by the variable NR. The radial extent of each field is set using the corresponding R_X variable, which stores a pair of increasing positive real numbers defining the radial extent of the field. The NR grid points will be distributed between radii corresponding to the minimum and maximum of these values. Currently, only the magnetic field can extend beyond the velocity field, modeling conducting solid layers.

Example The following lines in xshells.par define the radial extent of the fields:
R_U = 7/13 : 20/13
R_B = 0.0 : 20/13
R_T = 7/13 : 20/13

The default grid refines the number of points in the boundary layers, and this refinement can be controlled by the variable N_BL that stores a pair of integers, the first and second being the number of points reserved for the inner and outer boundary layer respectively, reinforcing the normal refinement. The code generating the grid can be found in the grid.cpp file.

Example The following lines in xshells.par define a grid with a total of 100 radial grid points, with 10 and 5 points reserved to the refinement of the inner and outer boundary layer respectively.
NR = 100
N_BL = 10,5

Alternatively, a radial grid can be loaded:

  • from a text file containing the radius of each grid point (increasing) in a separate line.

  • from a previously saved field (see section 1.7).

In both cases, simply indicate the filename in the r variable. It will override the NR and N_BL variables.

Example The following line in xshells.par will use the same grid as the field stored in file fieldU_0001.previous
r = fieldU_0001.previous      # load grid from file

Angular grid and spherical harmonic truncation

XSHELLS uses spherical harmonics to represent fields on the sphere:

f(θ,ϕ)=m=0M=mKLfmKYmK(θ,ϕ) (2.8)

where Ym is the spherical harmonic of degree and order m. The expansion uses a K-fold symmetry in longitude (ϕ) and is truncated at maximum degree L and order MK. If K=1 and M=L, it is the standard triangular truncation.

L, M and K are set in xshells.par using the Lmax, Mmax and Mres variables respectively. You must ensure that LMK.

The angular grid (spanning the co-latitude θ and longitude ϕ) consists of Nphi regularly spaced points in longitude, and Nlat gauss nodes in latitude. If these are not specified, XSHELLS will choose the values for Nlat and Nphi in order to ensure best performance and no aliasing of modes (Nlat >3L/2 and Nphi >3M).

Example These lines limit the spherical harmonic degree to 170. A 3-fold symmetry is used, and the maximum harmonic order is 56×3=168.
Lmax = 170   # max degree of spherical harmonics
Mmax = 56    # max fourier mode (phi)
Mres = 3     # phi-periodicity.
Most likely, 180 regularly spaced points in longitude and 256 gauss nodes in latitude will be used here.

2.1.6 Time-stepping

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-solution11 1 before version 2.0, a classical CNAB scheme was used.

  • CNAB2 is the classical second order Crank-Nicolson-Adams-Bashforth (CNAB) scheme22 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 scheme33 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 order44 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 accuracy55 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 × 2× 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 ΔT=2× sub_iter × 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 ΔT=0.01×25×2=0.5 time units (an iteration). The program will stop after iter_max=300 outputs (or iterations), spanning a total physical time of tend-tstart=150.0. Partial fields are saved every 10 iterations, or every 5.0 physical time units, if movie = 1 is set (see below).

2.1.7 Real time plotting

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.

2.1.8 Time lapse field snapshots

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.

2.1.9 Checkpointing and restarting capabilities

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.

2.1.10 Advanced options

Fine tuning of the automatic time-step selection

is possible through some variables.

C_u is a safety factor for the standard CFL (based on the velocity and the grid size), which depends on the time-scheme. For the PC2 time-scheme used by xshells, the default value is C_u = 0.65. You may need a more stringent value for a given problem, but if your simulation still does not work for a value lower than 0.1, the problem is probably elsewhere. C_alfv control the time-step adjustment (active if dt_adjust=1), regarding Alfvén criteria (due to magnetic field). The lower the values of C_u and C_alfv, the smaller the adjusted time step will be.

In addition, to prevent too many time step adjustments, if dt_tol_lo < dt/dt_target < dt_tol_hi, no time-step adjustment is done.

Example
C_u = 0.65         # default: 0.65
C_cori = 1         # default: 1
C_alfv = 2.8       # default: 2.8
dt_tol_lo = 0.9    # default: 0.9
dt_tol_hi = 1.05   # default: 1.05
Variable spherical harmonic degree truncation

is controlled by setting the variable rsat_ltr in xshells.par. This is possible only if XSHELLS was compiled with variable truncation enabled (see section 2.2.4).

Spectral convergence

of fields is checked by comparing the maximum of the end of the spectrum (four last modes) with the maximum of the spectrum. The first few modes are ignored (because their amplitude can be related directly to forcing or imposed fields). The variables sconv_lmin and sconv_mmin that can be set in xshells.par define the first spherical harmonic degree and order respectively that are considered when checking spectral convergence in the or m spectra. The default value is 3 for both variables. Setting a value larger than Lmax disables spectral convergence checks. Note also that the convergence checks are skipped when only a few modes are computed.

In addition, the maximum allowed value for Sconv can be set using variable sconv_stop (default value is 0.3). When the computed Sconv is larger than this limit, the program stops. Setting sconv_stop = 0 disables this safety check.

Example
sconv_lmin = 3     # ignore l=0-2 for spectral convergence.
sconv_mmin = 1     # ignore m=0 for spectral convergence.
sconv_stop = 0.3   # set the maximum value of Sconv before the program stops.
The SHTns library can be controlled

in terms of algorithm used for transforms, and in terms of polar optimization threshold.

The sht_type variable allows to constrain the transform method used:

  • sht_type = 0 : select fastest method using a classic Gauss-Legendre grid (default setting).

  • sht_type = 1 : select fastest method, allowing also regular grids (with DCT) which may be faster for small Mmax.

  • sht_type = 2 : impose a regularly spaced grid (not recommended as it is often slower).

  • sht_type = 3 : force a regularly spaced grid using DCT (not recommended as it is often slower).

  • sht_type = 4 : debug mode; initialization time is reduced, but a default method is used (no selection of fastest method).

  • sht_type = 6 : use a Gauss-Legendre grid with on-the-fly computation (preferred when parallel execution or big resolutions).

Finally, the polar optimization threshold can be adjusted with sht_polar_opt_max, the value below which coefficients near the pole are neglected. To give the reader some more insight, here are a some possible values and their impact:

  • sht_polar_opt_max = 0 : no polar optimization.

  • sht_polar_opt_max = 1e-14 : very safe optimization (default).

  • sht_polar_opt_max = 1e-10 : safe optimization.

  • sht_polar_opt_max = 1e-6 : aggressive optimization.

OpenMP tasks with dependencies

are used by xsbig_hyb whenever the number of shells cannot be distributed evenly enough between the available threads. However, OpenMP tasks with dependencies can cause problems with some compilers (e.g intel compilers prior to version 18), and can be turned off by setting omp_tasks = -1, or forced with omp_tasks = 1 in the xshells.par file.

2.1.11 Beta features

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.

Non-linear terms in linear mode

can be included to compute the saturation of an instability growing on a base flow. The program must be compiled in linear mode (see 2.2.7), and each non-linear term can be activated separately using a comma separated list in xshells.par:

nonlin = ugu,uxb,jxb,ugt,ugc

where ugu, uxb, jxb, ugt and ugc activate respectively 𝐮.𝐮, ×(𝐮×𝐛), (×𝐛)×𝐛, 𝐮T, 𝐮C for the perturbations 𝐮,𝐛,T and C in the equations (see 2.1.1).

Example To compute a kinematic dynamo growing on a saturated hydrodynamic instability, use:
nonlin = ugu,uxb     # include only u.grad(u) and curl(u x b)

2.2 Compile-time settings: xshells.hpp

All the following settings can be found in xshells.hpp. You have to recompile the program if you change this file.

2.2.1 Custom diagnostics

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.

2.2.2 Variable time-step

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).

2.2.3 Variable conductivity

In equation 2.2, conductivity can depend on radius r. To define a conductivity profile η(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 η(r) properly and automatically.

2.2.4 Variable spherical harmonic degree truncation

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 α in the formula to determine the truncation degree tr:

tr(r)=Lmaxrαrmax+1

where Lmax is defined by Lmax in xshells.par and rmax is the radius of the last shell. Note that tr cannot exceed Lmax.

Note also that α can be overridden by the rsat_ltr variable in xshells.par.

2.2.5 Hyper-diffusivities

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):

ν()={ν0for lcν0q-cfor >c (2.9)

with q=(νmax/ν0)1/(Lmax-c). The cutoff c and the ratio νmax/ν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 =Lmax 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 c is the same for all diffusivities.

2.2.6 Boundary forcing

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

Time dependent boundary forcing

are defined in the function calc_Uforcing(), found in the xshells.hpp file. In this function you must define a name for your forcing through the macro U_FORCING. The angular velocity of the solid bodies (defining the boundary of the fluid shell) can be set in this function. It will be used as a boundary condition for the flow if no-slip boundaries are used (see section 2.1.2).

See the example found in the problems/couette/ folder for more details, and uncomment the part of the function corresponding to your problem. In particular axial differential rotation of the inner or outer boundary can be used to simulate a spherical Couette flow; equatorial differential rotation (or spin-over) can be used to simulate precession or nutation.

Note that the rotation rate of the solid bodies is also used in the induction equation if the magnetic field extends into the solids (conducting solid shells).

Arbitrary stationary boundary flows

You can impose arbitrary stationary flows at the solid boundaries. Uncomment:

#define XS_SET_BC

and change the set_U_bc() function found in xshells.hpp according to your needs. Note that the boundary conditions for the poloidal velocity field is stored in the shell beyond the first or last fluid shells (respectively NG-1 and NM+1). See for example the xshells.hpp file in the problems/full_sphere/ folder. Note that the solid should not be conducting if this feature is used, as no corresponding solid flow will be generated.

2.2.7 Linearized equations and base fields

To replace the equations with their linearized version (no (×𝐮)×𝐮, no (×𝐛)×𝐛, no 𝐮.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.

2.2.8 Mean-field dynamo

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.