4. Running Aronnax

To run a simulation with Aronnax one needs to have a Python session active in the folder for the simulation. This folder should contain a file, aronnax.conf that contains the configuration choices for the simulation. Some, or all, of these choices may be overridden by arguments passed to the simulate function, but it is likely simpler to specify many of the choices in a configuration file.

As described above, it is possible to define functions that can be passed to aronnax.driver.simulate and used to create input or forcing fields. The test suite, found in the ‘test’ folder uses this functionality to create the zonal wind stress for the \(\beta\)-plane gyre tests. The relevant code is shown below:

def wind(_, Y):
    return 0.05 * (1 - np.cos(2*np.pi * Y/np.max(grid.y)))

with working_directory(p.join(self_path, "beta_plane_gyre_red_grav")):
    drv.simulate(zonal_wind_file=[wind],
                 nx=10, ny=10, exe="aronnax_test", dx=xlen/10, dy=ylen/10)

Warning

Parameters cannot be set to 0 or 1 in the call to drv.simulate because the Python wrapper gets confused between numerical and logical variables, as described in https://github.com/edoddridge/aronnax/issues/132

4.1. Executables

Aronnax includes multiple targets for the Makefile. These produce executables intended either for testing or production runs. The different options for exe are discussed below. For a comparison on the execution speed of these executables see Benchmarking.

For production runs, and simulations that span multiple processors, use either aronnax_core (for \(n+1/2\)-layer mode) or aronnax_external_solver (for \(n\)-layer mode).

  • aronnax_core

    • Uses the internal Fortran code and aggressive compiler optimisation. This is the best executable to use for reduced gravity (\(n+1/2\) layer) simulations. It may also be used for \(n\) layer simulations but is considerably slower than aronnax_external_solver since it does not use the external Hypre library for the pressure solve. aronnax_core cannot run on multiple processors in \(n\)-layer mode.

  • aronnax_test

    • Uses the internal Fortran code and no compiler optimisations. This is the best executable to use for assessing code coverage when running tests. The Fortran error messages might be helpful.

  • aronnax_prof

    • Executable intended solely for profiling the Fortran core. Unlikely to be of general use.

  • aronnax_external_solver_test

    • Uses the external Hypre library to solve the matrix inversion in the \(n\) layer mode. Uses no optimisations for the internal Fortran code and is intended for assessing code coverage when running tests.

  • aronnax_external_solver

    • Uses the external Hypre library and aggressive compiler optimisations. This is the fastest executable to use in \(n\) layer mode. It can also be used in \(n+1/2\) layer mode, but there is no advantage over aronnax_core.

4.2. Parameters

Parameters can be passed to the model in two ways. Either they can be included in a file called aronnax.conf in the working directory, or they can be passed as keyword arguments to aronnax.driver.simulate(). The main directory of the repository includes an example aronnax.conf file.

Note

Aronnax takes a deliberately strict approach towards specifying parameter choices. The model contains very few default values, except in situations where a default of zero can sensibly be applied. As such, you will need to specify your parameter choices, either through the configuration file, or the function call. However, parameters that are not required, for example, bottom drag in n+1/2 layer mode, need not be set.

The example file is reproduced here with comments describing the parameters and their units. All possible parameters are included, but they are not all assigned values. After modifying this file for your simulation, any unassigned parameters should be deleted.

# Aronnax configuration file. Change the values, but not the names.

#------------------------------------------------------------------------------
# au is the lateral friction coefficient in m^2 / s
# ar is linear drag between layers in 1/s
# kh is thickness diffusivity in m^2 / s
# kv is vertical thickness diffusivity in m^2/s
# dt is time step in seconds
# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
# niter0: the timestep at which the simulation begins. If not zero, then there 
#   must be checkpoint files in the 'checkpoints' directory.
# n_time_steps: number of timesteps before stopping
# dump_freq: time between snapshot outputs in seconds
# av_freq: time between averaged output in seconds
# checkpiontFreq: time between checkpoints in seconds 
#   (these are used for restarting simulations)
# diag_freq: time between dumping layerwise diagnostics of the simulation. These
#   are mean, min, max, and std of h, u, and v in each layer.
# hmin: minimum layer thickness allowed by model (for stability) in metres
# maxits: maximum iterations for the pressure solver algorithm. Should probably
#   be at least max(nx,ny), and possibly nx*ny
# eps: convergence tolerance for pressure solver. Algorithm stops when error is 
#   less than eps*initial_error
# freesurf_fac: 1. = linear implicit free surface, 0. = rigid lid.
# bot_drag is the linear bottom friction in 1/s
# thickness_error is the  discrepancy between the summed layer thicknesses and 
#   the depth above which the model emits a warning. 1e-2 is a 1% discrepancy.
# debug_level controls the level of output produced by the model. When set to
#   zero, or not specified (it defaults to zero), the model outputs h, u, and v
#   at the frequency controlled by dump_freq and av_freq. Specifying larger
#   integer values results in progressively more output more frequently. See
#   the documentation for details.
# h_advec_scheme selects which advection scheme to use when advecting the
#   thickness field. Current options are:
#   1 first-order centered differencing
#   2 first-order upwind differencing
# ts_algorithm selects the time stepping algorithm used by the model
#   3 (default) is third-order Adams-Bashfort

[numerics]
au = 500.
ar = 0.0
kh = 0.0
kv = 0.0
dt = 600.
slip = 0.0
niter0 = 0
n_time_steps = 502
dump_freq = 1.2e5
av_freq = 1.2e5
checkpoint_freq = 1.2e5
diag_freq = 6e3
hmin = 100
maxits = 1000
eps = 1e-5
freesurf_fac = 0.
bot_drag = 1e-6
thickness_error = 1e-2
debug_level = 0
h_advec_scheme = 2
ts_algorithm = 3
#------------------------------------------------------------------------------

# red_grav selects whether to use n+1/2 layer physics (red_grav=yes), or n-layer 
#   physics with an ocean floor (red_grav=no)
# active_lower_layer is used in reduced gravity mode to swap from a quiescent abyss
# (the default) to an active lower layer with a quiescent upper layer by setting
# active_lower_layer = yes. You will need to provide bathymetry via `depth_file`.
# depth_file defines the depth of the ocean bottom below the sea surface in metres.
# hmean is a list of initial thicknesses for the layers in metres. Each value is 
#   separated by a comma. This input was a useful short cut for specifying 
#   initial conditions with constant layer thicknesses, but has been superseded 
#   and may be removed in the future.
# h0 is the depth of the ocean basin and is only required in n-layer mode. This 
#   input was a useful short cut for specifying a flat bottomed ocean, but has 
#   been superseded and may be removed in the future.

[model]
red_grav = no
active_lower_layer = no
depth_file
hmean = 400.,1600.
h0 = 2000.
#------------------------------------------------------------------------------

# these variables set the number of processors to use in each direction. 
#   You should ensure that the grid divides evenly into the number of tiles.
#   nx/n_proc_x and ny/n_proc_y must be integers.

[pressure_solver]
n_proc_x = 1
n_proc_y = 1
#------------------------------------------------------------------------------

# g_vec is the reduced gravity at interfaces in m/s^2. g_vec must have as many 
#   entries as there are layers. The values are given by the delta_rho*g/rho_0. 
#   In n-layer mode the first entry applies to the surface, i.e. the top of the 
#   upper layer. In n+1/2 layer mode the first entry applies to the bottom of 
#   the upper layer.
# rho0 is the reference density in kg/m^3, as required by the Boussinesq assumption.

[physics]
g_vec = 9.8, 0.01
rho0 = 1035.
#------------------------------------------------------------------------------

# nx is the number of grid points in the x direction
# ny is the number of grid points in the y direction
# layers is the number of active layers
# OL is the width of the halo region around tiles
#   must be set to >= 3 if using multiple tiles
# dx is the x grid spacing in metres
# dy is the y grid spacing in metres
# f_u_file defines the Coriolis parameter on the u grid points in 1/s
# f_v_file defines the Coriolis parameter on the v grid points in 1/s
# wet_mask_file defines the computational domain - which grid points are ocean, 
#   with wetmask=1, and which are land, with wetmask=0. The wetmask is defined 
#   at the centre of each grid cell, the same location as thickness.

[grid]
nx = 10
ny = 10
layers = 2
OL = 3
dx = 2e4
dy = 2e4
f_u_file = :beta_plane_f_u:1e-5,2e-11
f_v_file = :beta_plane_f_v:1e-5,2e-11
wet_mask_file = :rectangular_pool:
#------------------------------------------------------------------------------

# These files define the values towards which the model variables are relaxed 
#   (in metres or m/s), and the timescale for the relaxation, in 1/s.
[sponge]
sponge_h_time_scale_file
sponge_u_time_scale_file
sponge_v_time_scale_file
sponge_h_file
sponge_u_file
sponge_v_file

#------------------------------------------------------------------------------

# These files define the initial values used in the simulation. If no values are
#   defined for the velocities (in m/s) or the free surface elevation (in m), 
#   they will be initialised with zeros. Layer thickness (in m) must be initialised, 
#   either by passing a file, or using the generator functions.

[initial_conditions]
init_u_file
init_v_file
init_h_file
init_eta_file

#------------------------------------------------------------------------------

# The wind files define the momentum forcing in N/m^2 or m/s
# wind_mag_time_series_file defines the constant factor by which the wind is 
#   multiplied by at each timestep.
# wind_depth specifies the depth over which the wind forcing is spread. This
#   only impacts the simulation if the surface layer is thinner than wind_depth.
#   If wind_depth is set to 0 (default) then the wind forcing acts only on the 
#   surface layer, no matter how thin it gets.
# dump_wind defines whether the model outputs the wind field when it outputs other 
#   variables (at the rate controlled by dump_freq).
# relative_wind selects whether the wind forcing is given in 
#   N/m^2 (relative_wind = no), or
#   m/s   (relative_wind = yes)
# Cd is the quadratic drag coefficient used if relative_wind = yes
# wind_n_records is the number of wind snapshots provided
# wind_period is the time (in seconds) between subsequent snapshots
# wind_loop_fields sets whether to loop through the provided snapshots, or
#   to stay with the final one for the rest of the simulation
# wind_interpolate sets whether to interpolate between snapshots, or swap 
#   from one to the next in a single timestep

[external_forcing]
zonal_wind_file = 'wind_x.bin'
meridional_wind_file = 'wind_y.bin'
wind_mag_time_series_file
wind_depth = 30
dump_wind = no
relative_wind = no
Cd = 0.
wind_n_records = 1
wind_period = 1000
wind_loop_fields = yes
wind_interpolate = yes

#------------------------------------------------------------------------------

Warning

The configuration file shown above includes all of the possible input parameters and fields since it forms part of the documentation. IT WILL NOT WORK AS IS. To use it in an actual simulation the file will need to be modified either by giving values to the parameters that are currently unspecified, or deleting them from the file. If you wish to see a configuration file that corresponds to a successful simulation, look in any of the test, benchmark, or reproduction directories.

4.2.1. debug_level

This parameter determines whether the model produces additional outputs. It should be set to an integer value greater than or equal to zero. The different values have the following effects:

  • 0: no additional outputs. Output frequency controlled by dump_freq and av_freq

  • 1: output tendencies at frequency given by dump_freq

  • 2: output tendencies and convergence diagnostics from the linear solve at frequency given by dump_freq (not implemented)

  • 3: output convergence diagnostics and tendencies before and after applying some or all of sponges, barotropic correction, winds, and boundary conditions at frequency controlled by dump_freq (not implemented)

  • 4: dump all of the above fields every time step (mostly implemented)

  • 5: dump everything every time step including the two initial RK4 steps (not implemented)

4.2.2. niter0

This parameter allows a simulation to be restarted from the given timestep. It requires that the appropriate files are in the ‘checkpoints’ directory. All parameters, except for the number of grid points in the domain, may be altered when restarting a simulation. This is intended for breaking long simulations into shorter, more manageable chunks, and for running perturbation experiments.

4.2.3. wet_mask_file

The wetmask defines which grid points within the computational domain contain fluid. The wetmask is defined on the tracer points, and a value of 1 defines fluid, while a value of 0 defines land. The domain is doubly periodic in x and y by default. To produce a closed domain the wetmask should be set to 0 along the edges of the domain. Placing a barrier at either the southern or the northern boundary will remove periodicity in the meridional direction. Similarly, a barrier along either the eastern or western boundary will prevent the model from being periodic in the zonal direction.

4.2.4. relative_wind

If this is false, then the wind input file is given in N m–2. If true, then the wind input file is in m s–1 and a quadratic drag law is used to accelerate the fluid with Cd as the quadratic drag coefficient.

4.2.5. h_advec_scheme

h_advec_scheme is an integer that selects the advection scheme used for thickness in the continuity equation. Currently two options are implemented:

  • h_advec_scheme = 1 (default) uses a first-order centered stencil

  • h_advec_scheme = 2 uses a first-order upwind stencil

4.2.6. ts_algorithm

ts_algorithm is an integer that selects the timestepping algorithm used by the model. The default behaviour is to use a third-order Adams-Bashforth scheme (ts_algorithm = 3), with the initialisation performed by a second-order Runge-Kutta method (ts_algorithm = 12).

  • ts_algorithm = 1: Forward Euler

  • ts_algorithm = 2: Second-order Adams-Bashfort

  • ts_algorithm = 3: Third-order Adams-Bashfort (default)

  • ts_algorithm = 4: Fourth-order Adams-Bashfort

  • ts_algorithm = 5: Fifth-order Adams-Bashfort

  • ts_algorithm = 12: Second-order Runge-Kutta

  • ts_algorithm = 13: Third-order Runge-Kutta (not implemented)

  • ts_algorithm = 14: Fourth-order Runge-Kutta (not implemented)

4.3. More complicated forcing fields

It is possible to use surface forcing fields that vary both in space and time.

If the pattern of the surface forcing is always the same, the wind_mag_time_series_file parameter can be used to make the magnitude of the forcing vary. The provided wind_mag_time_series_file needs an entry for each timestep of the simulation. Like all inputs, this can be specified as a Python function that will be evaluated when the model initialises and the output saved to disk so that the Fortran core can obtain the input.

Spatially varying forcing is slightly more complex. In this situation zonal_wind_file and meridional_wind_file will now point to three-dimensional arrays. These arrays can be generated by the Python wrapper if a function is passed in the call to aronnax.driver.simulate(). This function must depend only on X, Y, wind_n_records. The parameters that must be specified are:

  • wind_n_records: the number of time slices in the provided wind field. An integer.

  • wind_period: the time (in seconds) between subsequent wind forcing patterns. A float.

  • wind_loop_fields: whether the wind forcing should loop back through the provided fields. A Boolean.

  • wind_interpolate: whether the wind forcing should be interpolated between slices (provides a continuously varying field), or jump discontinuously between the slices. A Boolean.

4.4. Discussion

The following points may be helpful when using the model but don’t fit above.

4.4.1. Outcropping

Aronnax allows for layers to become very thin, which simulates outcropping of isopycnals at the surface, and grounding of isopycnals at the sea floor. Unfortunately outcropping adversely impacts layerwise volume conservation. Over a long simulation the change in volume of a layer may be substantial.

To allow for outcropping the following parameters should be set:

  • hmin very small

  • wind_depth set to >10 m