Example configuration files can be found in the problems directory.
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.
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:
is the velocity field.
and are respectively the temperature and concentration fields, that contribute to the buoyancy in the Boussinesq formulation.
is the magnetic field ( and 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.
is the magnetic diffusivity of the fluid ( is its conductivity) and is set by the variable eta in xshells.par. Diffusivity depending on shell radius are also supported (see sec. 2.2.3).
and 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 plane ().
is the gravity potential (independent of time), controlled by field phi0 in xshells.par.
and are the imposed base temperature and concentration profiles, controlled by fields tp0 and c0 in xshells.par.
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.
If the following lines are found in xshells.par:
nu = 1.0 eta = sqrt(10) Omega0 = 2*pi*1e3 u = 0 b = 0 #tp = 0then the viscosity is set to , the magnetic diffusivity is set to , and the rotation rate is set to . 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.
The Lorentz force can be turned off. Just add no_jxb = 1 in the xshells.par file.
Vector fields are represented internally using a poloidal/toroidal decomposition:
where is the radial position vector, and and are the toroidal and poloidal scalars respectively. This decomposition ensures that the vector field is divergence-free.
The scalar fields and for each radial shell are then decomposed on the basis of spherical harmonics.
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).
boundary conditions can be fixed value (defined by the or profile), fixed flux (defined by or ), or Robin (such as ). 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).
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.
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
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).
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)
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).
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.
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 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.
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.
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.
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
XSHELLS uses spherical harmonics to represent fields on the sphere:
where is the spherical harmonic of degree and order . The expansion uses a -fold symmetry in longitude () and is truncated at maximum degree and order . If and , it is the standard triangular truncation.
, and are set in xshells.par using the Lmax, Mmax and Mres variables respectively. You must ensure that .
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 and Nphi ).
These lines limit the spherical harmonic degree to 170. A 3-fold symmetry is used, and the maximum harmonic order is .
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.
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..
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 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 sub_iter dt.
Finally, iter_save controls the number of iterations before a (partial) snapshot is saved to disk.
The following lines in xshells.par will use a time-step of 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 writesOutput will occur every time units (an iteration). The program will stop after iter_max=300 outputs (or iterations), spanning a total physical time of . 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.
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.
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.
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.
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
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).
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 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.
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.
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.
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.
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.
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 , , , , for the perturbations ,, and in the equations (see 2.1.1).
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)
All the following settings can be found in xshells.hpp. You have to recompile the program if you change this file.
Enable by uncommenting:
In addition to the total energy of the three fields , and , 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:
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 . To define a conductivity profile , uncomment:
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 properly and automatically.
In order to compute in a full sphere and avoid problems near , the spherical harmonic expansion must be truncated at low degree near . 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 ( in the line above, which is a good choice for full-sphere computations) is used as in the formula to determine the truncation degree :
where is defined by Lmax in xshells.par and is the radius of the last shell. Note that cannot exceed .
Note also that 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:
The formulae used to compute the viscosity is the one proposed by Nataf & Schaeffer (2015):
with . The cutoff and the ratio 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 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 is the same for all diffusivities.
Amplitude and frequency are set at runtime by a_forcing and w_forcing in the xshells.par file.
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).
You can impose arbitrary stationary flows at the solid boundaries. Uncomment:
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.
To replace the equations with their linearized version (no , no , no ), uncomment:
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:
And provide the alpha field using the init_alpha function. See the example located in problems/kinematic_dynamos.