Getting started

Here we provide general information on setting up and running PyLag. For specific applications - e.g., running PyLag with data generated by a given hydrodynamic model, such as FVCOM - please see PyLag’s tutorial examples.

After PyLag has been installed, a typical workflow involves three steps: 1) Obtaining and/or generating the required input files; 2) Running the simulation; and 3) Analysing the results. Further details regarding these three steps are provided below.

Step 1: Obtaining and/or generating the set of input files required by PyLag

For typical use cases, PyLag requires four types of input:

  1. NetCDF data files containing input data (e.g. velocity components) required by PyLag.

  2. A grid metrics file, which describes the underlying grid on which input data are defined.

  3. An initial positions file, which gives the starting coordinates of particles in Cartesian or Geographic coordinates.

  4. A run configuration file in which all run parameters and paths are set.

1) NetCDF data files

NetCDF data files containing fluid velocity components, eddy diffusivities etc are required inputs. All ocean specific variables for a given time point or set of time points should be provided in the same file. The dataset may be split into multiple files through time, with all files following a common naming convention (see below). When support for including the impact of waves and direct atmospheric forcing is included, the set of input files will be extended to include files specific to wave and astmopheric products respectively.

If the set of ocean specific variables have been generated by one of the named models supported by PyLag (e.g. ROMS, FVCOM) then the data files can be fed to PyLag as is. This will typically be the case for data obtained from other supported sources too (e.g. the CMEMS catalogue). However, some manipulation may be needed if, say, the set of ocean variables have been saved in separate files.

All input files must contain a time variable with a units attribute, allowing times to be converted into datetime objects. The files should also contain grid information (e.g. depth, latitude and longitude arrays). Any number of input data files can be provided. With respect to file name conventions, the primary requirement is that they are named using a pattern that allows them to be sorted sequentially through time. For example, two data files for January and February 2002 could be named file_001.nc and file_002.nc. Using appropriately constructed datestrings in the file name is usually sufficient to guarentee sequential ordering. PyLag will automatically check that the simulation start and end times lie within the time range covered by the input files, and that at least one of the particles in the initial particle seed lies within the spatial grid.

As an example, the header for a NetCDF file downloaded from the CMEMS catalogue is given below. The file contains information from the GLOBAL_REANALYSIS_PHY_001_031 product. The file contains grid information and u and v velocity components for a global reanalysis run at 1/4 degree resolution. From the global metadata, it can be seed data from multiple months has been concatenated into a single file for convenience.

NetCDF header for CMEMS data product GLOBAL_REANALYSIS_PHY_001_031

netcdf grepv2_monthly_surface_2016_2017 {
dimensions:
    depth = 1 ;
    latitude = 681 ;
    longitude = 1440 ;
    time = UNLIMITED ; // (24 currently)
variables:
    float depth(depth) ;
        depth:valid_min = 0.50576f ;
        depth:valid_max = 5902.058f ;
        depth:units = "m" ;
        depth:positive = "down" ;
        depth:unit_long = "Meters" ;
        depth:long_name = "Depth" ;
        depth:standard_name = "depth" ;
        depth:axis = "Z" ;
    float latitude(latitude) ;
        latitude:valid_min = -80.f ;
        latitude:valid_max = 90.f ;
        latitude:step = 0.25f ;
        latitude:units = "degrees_north" ;
        latitude:unit_long = "Degrees North" ;
        latitude:long_name = "Latitude" ;
        latitude:standard_name = "latitude" ;
        latitude:axis = "Y" ;
    float longitude(longitude) ;
        longitude:valid_min = -180.f ;
        longitude:valid_max = 179.75f ;
        longitude:step = 0.25f ;
        longitude:units = "degrees_east" ;
        longitude:unit_long = "Degrees East" ;
        longitude:long_name = "Longitude" ;
        longitude:standard_name = "longitude" ;
        longitude:axis = "X" ;
    float time(time) ;
        time:axis = "T" ;
        time:long_name = "Time" ;
        time:standard_name = "time" ;
        time:calendar = "gregorian" ;
        time:units = "days since 1950-01-01 00:00:00" ;
        time:cell_methods = "time: mean" ;
    float uo_foam(time, depth, latitude, longitude) ;
        uo_foam:long_name = "Eastward velocity" ;
        uo_foam:standard_name = "eastward_sea_water_velocity" ;
        uo_foam:units = "m s-1" ;
        uo_foam:unit_long = "Meters per second" ;
        uo_foam:_FillValue = 9.96921e+36f ;
        uo_foam:cell_methods = "area: mean time: mean" ;
        uo_foam:valid_min = -3.f ;
        uo_foam:valid_max = 3.f ;
    float vo_foam(time, depth, latitude, longitude) ;
        vo_foam:long_name = "Northward velocity" ;
        vo_foam:standard_name = "northward_sea_water_velocity" ;
        vo_foam:units = "m s-1" ;
        vo_foam:unit_long = "Meters per second" ;
        vo_foam:_FillValue = 9.96921e+36f ;
        vo_foam:cell_methods = "area: mean time: mean" ;
        vo_foam:valid_min = -3.f ;
        vo_foam:valid_max = 3.f ;
// global attributes:
        :product = "GLOBAL_REANALYSIS_PHY_001_031" ;
        :producer = "CMEMS - Global Monitoring and Forecasting Centre" ;
        :area = "Global" ;
        :quality_information_document = "http://marine.copernicus.eu/documents/QUID/CMEMS-GLO-QUID-001_031.pdf" ;
        :Conventions = "CF-1.6" ;
        :credit = "E.U. Copernicus Marine Service Information (CMEMS)" ;
        :contact = "servicedesk.cmems@mercator-ocean.eu" ;
        :references = "http://marine.copernicus.eu" ;
        :licence = "http://marine.copernicus.eu/services-portfolio/service-commitments-and-licence/" ;
        :institution = "Mercator Ocean" ;
        :product_user_manual = "http://marine.copernicus.eu/documents/PUM/CMEMS-GLO-PUM-001-031.pdf" ;
        :title = "Monthly mean fields for product GLOBAL_REANALYSIS_PHY_001_031" ;
        :history = "Tue Mar  9 08:11:13 2021: ncrcat grepv2_monthly_surface_201601.nc grepv2_monthly_surface_201602.nc grepv2_monthly_surface_201603.nc grepv2_monthly_surface_201604.nc grepv2_monthly_surface_201605.nc grepv2_monthly_surface_201606.nc grepv2_monthly_surface_201607.nc grepv2_monthly_surface_201608.nc grepv2_monthly_surface_201609.nc grepv2_monthly_surface_201610.nc grepv2_monthly_surface_201611.nc grepv2_monthly_surface_201612.nc grepv2_monthly_surface_201701.nc grepv2_monthly_surface_201702.nc grepv2_monthly_surface_201703.nc grepv2_monthly_surface_201704.nc grepv2_monthly_surface_201705.nc grepv2_monthly_surface_201706.nc grepv2_monthly_surface_201707.nc grepv2_monthly_surface_201708.nc grepv2_monthly_surface_201709.nc grepv2_monthly_surface_201710.nc grepv2_monthly_surface_201711.nc grepv2_monthly_surface_201712.nc grepv2_monthly_surface_2016_2017.nc\n",
            "Tue Mar  9 08:07:09 2021: ncks -v vo_foam,uo_foam -d depth,0,0 /data/thaumus2/scratch/jcl/CMEMS_global_phy_reanalysis_1_4_deg/ftp/monthly/all/grepv2_monthly_201601.nc ./grepv2_monthly_surface_201601.nc\n",
            "Creation 2019-Jan-16 08:00:00 GMT+0200" ;
        :dataset = "global-reanalysis-phy-001-031-grepv2-monthly" ;
        :source = "Copernicus Marine Service" ;
        :NCO = "4.6.8" ;
        :nco_openmp_thread_number = 1 ;
}

Example files for ROMS and FVCOM, which are used in PyLag’s tutorial examples, can be downloaded here. You will see the header information for the different types of input data varies considerably. For the most part, PyLag deals with this variation internally, or through run configuration settings.

The grid metrics file

The second input requirement is a grid metrics file in NetCDF format. It is required for all runs involving 2D or 3D spatial data. Within it are a set of variables describing an unstructured triangular mesh which has been constructed from the underlying grid. In the case of FVCOM, which already employs an unstructured grid, the grid metrics file saves variables describing the same grid, but with node and neighbour elements reordered to match PyLag’s requirements. For regular grids, a new unstructured grid is created from scratch.

The grid metrics is a requirement as it saves on having to compute the same information at the start of each run. It provides all the information PyLag needs to correctly read in and interpret the gridded input data. Special functions are provided in the module `pylag.grid_metrics <../api/pylag.grid_metrics.rst>`__ which can be used to create grid metrics files that are specific to the type of input data being used. As an example, a grid metrics file for the above CMEMS data can be created so:

>>> import pylag
>>> pylag.grid_metrics.create_arakawa_a_grid_metrics_file('./grepv2_monthly_surface_2016_2017.nc', reference_var_name='uo_foam', is_global=True, surface_only=True)

A full explanation of each argument can be found in the module documentation. The above function will generate a file called grid_metrics.nc. The corresponding NetCDF header is:

NetCDF header for a PyLag grid metrics file constructed from CMEMS input data

netcdf grepv2_grid_metrics_surface_only {
dimensions:
    three = 3 ;
    longitude = 1440 ;
    latitude = 680 ;
    node = 979200 ;
    element = 1958396 ;
variables:
    double longitude(node) ;
        longitude:units = "degrees_east" ;
        longitude:standard_name = "longitude" ;
        longitude:long_name = "Longitude" ;
    double longitude_c(element) ;
        longitude_c:units = "degrees_east" ;
        longitude_c:standard_name = "longitude" ;
        longitude_c:long_name = "Longitude" ;
    double latitude(node) ;
        latitude:units = "degrees_north" ;
        latitude:standard_name = "latitude" ;
        latitude:long_name = "Latitude" ;
    double latitude_c(element) ;
        latitude_c:units = "degrees_north" ;
        latitude_c:standard_name = "latitude" ;
        latitude_c:long_name = "Latitude" ;
    int64 trim_first_latitude ;
        trim_first_latitude:long_name = "0 - no, 1 - yes" ;
    int64 trim_last_latitude ;
        trim_last_latitude:long_name = "0 - no, 1 - yes" ;
    int64 permutation(node) ;
        permutation:long_name = "node permutation" ;
    int64 nv(three, element) ;
        nv:long_name = "nodes surrounding each element" ;
    int64 nbe(three, element) ;
        nbe:long_name = "elements surrounding each element" ;
    int64 mask_c(element) ;
        mask_c:standard_name = "sea_binary_mask" ;
        mask_c:units = "1" ;
        mask_c:long_name = "Land-sea mask: sea = 0, land = 1, boundary element = 2" ;
    int64 mask(node) ;
        mask:standard_name = "sea_binary_mask" ;
        mask:units = "1" ;
        mask:long_name = "Land-sea mask: sea = 0, land = 1" ;
    double area(element) ;
        area:standard_name = "areas" ;
        area:units = "m^2" ;
        area:long_name = "Element areas" ;

// global attributes:
        :Conventions = "CF-1.7" ;
        :title = "PyLag grid metrics file" ;
        :institution = "Plymouth Marine Laboratory (PML)" ;
        :contact = "James R. Clark (jcl@pml.ac.uk)" ;
        :netcdf-version-id = "netCDF-4" ;
        :pylag-version-id = "6a56a34d0649e2c029f5d38584548c2ea01eeb1f" ;
        :is_global = "True" ;
}

Most of the variables are self explanatory, given the variable long names. Possible exceptions are the variables trim_first_latitude and trim_last_latitude. These are binary flags, which indicate whether the first and/or last latitude points should be trimmed when PyLag reads input data from file. They are included as global data products are sometimes interpolated onto a regular grid which includes a row of points lying directly over one or both poles. On a spherical Earth, all these points refer to the same location, and the row is frequently masked. To prevent problems associated with handling such points, PyLag will trim them if they exist. These flags indicate whether this has happened or not. Other variables of note are perturbation, gives the order of nodes in the grid. It is used to sort input arrays so that data points are correctly matched to nodes that make up the unstructured grid. Lastly, the _c name qualifier is used to flag variables specific to element centres (as opposed to element nodes).

The initial positions file

PyLag must be given the initial positions of particles at the start of each new simulation. When running PyLag from the command line, particle initial positions can be specified in one of two ways. In the first approach, particle initial positions are read in from a text file with a specified format, as outlined below. In the second approach, particle initial positions are read from a restart file which has been created at the end of or during a previous simulation. The choice is specified in the run configuration file.

Initial position text files must have the following format:

n
group_id x1 x2 x3
group_id x1 x2 x3
...

The entries have the following meaning:

  • n is the total number of particles. It must always be the first entry in the file.

  • group_id is an integer number that can be used to group particles. This can be useful if, say, you want to look at a group of particles that started off from a similar location. It does not affect particle dynamics. It is recorded in the simulation output files.

  • x1 is the particle’s first position coordinate. In a cartesian coordinate system, this should be x in meters. In a geographic coordinates, it should be the particle’s longitude in degrees east.

  • x2 is the particle’s second position coordinate. In a cartesian coordinate system, this should be y in meters. In a geographic coordinates, it should be the particle’s latitude in degrees north.

  • x3 is the particle’s third position coordinate. This is always the particle’s depth, z, in meters. z can be given relative to the moving free surface or the sea floor. The choice is set in the run configuration file using the option depth_coordinates. Depth is positive up. Thus, if depth_coordinates is set to depth_below_surface, a depth, \(z\), of -0.1 would specify a starting z position 0.1 m below the free surface. The model will interpolate the free surface height to the particle’s position in time and space, so you don’t need to worry about what the sea surface height is at a particlular point in time (unless you set a starting depth that lies below the sea floor). There is also an option to start particles off from a specified distance above the sea floor.

For example, an initial positions file, initial_positions.dat, for a simulation involving two particles with the same group ID, whose lateral coordinates are given in degrees and initial vertical coordinates are coincident with the free surface, might look like:

2
1 -4.17 50.0 0.0
1 -4.17 50.1 0.0

For simulations involving tens of particles or more, the initial positions file is most easily created programatically.

The run configuration file

The final input required by PyLag is the run configuration file, where all settings for a given run are specified. Example run configuration files for each of the different types of input data that PyLag works with have been provided in the resources directory of the code repository. They can simply be downloaded and adapted to specific use cases. The template file for running with CMEMS like data on an Arakawa A-grid is given below. The file contains a detailed set of comments for each setting.

Template run configuration file for CMEMS like data

[GENERAL]
# Logging
log_level = INFO

# Directory containing input files
in_dir = ./input

# Directory in which simulation output will be saved
out_dir = ./output

# Name of the netCDF data file to be created (without the nc extension)
output_file = %(out_dir)s/pylag

[SIMULATION]
# Simulation type
simulation_type = trace

# Initialisation method
#   : init_file: Initialise particle state data from an ascii file. The
#     name of the file is given by the parameter `initial_positions_file`.
#   : restart_file: Initialise particle state data from a restart file. The
#     name of the restart file is given by the parameter `restart_file_name`
#     in section `RESTART`.
#   : rectangular_grid: Initialise particle state data by creating a set
#     of particles on a regular rectangular grid. (NOT IMPLEMENTED)
initialisation_method = init_file

# Initial positions data file
initial_positions_file = ./input/initial_positions.dat

# Flag for surface only transport. If set to True, PyLag will only read
# surface fields. In this case, initial vertical positions should be
# set to 0.0 m below the free surface. Depth restoring should be set to
# False for efficiency reasons.
surface_only = False

# The coordiate system. Options include cartesian (x, y) and geographic (lat, lon)
# coordinate systems. The choice should reflect the grid on which the underlying
# input data (e.g. u, v velocity components) are defined. PyLag will adopt the
# same coordinate system for particle positions. This choice does not impact
# the vertical coordinate system, where at the current time, regular (z) depth
# levels are assumed. Note that some models such as FVCOM can output particle
# coordinates in both cartesian (x, y) and geographic (lat, lon) coordinates.
# In these cases, and in paricular for small domains where grid distortions
# due to the projection are small, it is generally faster and better to use
# cartesian coordinates.
#     : cartesian - Cartesian (x, y) coordinates.
#     : geographic - Geographic (lat, lon) coordinates.
coordinate_system = geographic

# Depth coordinates
#   : depth_below_surface: Starting depth is given as the depth below the
#     (moving) free surface in meters. Positive up, meaning a starting depth
#     of -0.1 is 0.1 m below the free surface. NB A starting depth of 0.1 would
#     be above the free surface and flagged as an error.
#   : hieght_above_bottom: Depth is given as the height above the sea floor
#     in meters. Positive up, meaning a starting depth of 0.1 is 0.1 m above
#     the sea floor.
depth_coordinates = depth_below_surface

# Keep particles at a fixed depth below the surface by restoring to a fixed depth each time step
depth_restoring = False

# The fixed depth below the surface that particles are held at if `depth_restoring' is True. A value
# 0f 0.0 corresponds to the sea surface; a value of -1.0 to 1m below the free surface.
fixed_depth = 0.0

# The time and date at which the first set of particles is released. The format is: yyyy-mm-dd hh:mm:ss
start_datetime = 2016-01-12 00:00:00

# Simulation end time - only used when running a single particle release
end_datetime = 2016-03-12 00:00:00

# Time direction
#  : forward: Forward in time integration
#  : reverse: Backward in time integration
time_direction = forward

# Number of particle releases
number_of_particle_releases = 1

# Number of hours between particle releases
particle_release_interval_in_hours = 1.0

# Simulation duration - used to compute simulation end times with multiple particle releases
duration_in_days = 0.5

# Period at which output data should be saved (s)
output_frequency = 900.0

# Period at which output data is synced to disk (s)
sync_frequency = 900.0

[RESTART]
# The name of the restart file that will be used to initialise the model if
# initialisation_method == restart_file
restart_file_name = ./restart.nc

# Create restart files?
create_restarts = True

# The directory in which restart files will be created. Full or relative path.
restart_dir = ./restart

# Period at which restart files will be created (s)
restart_frequency = 3600.0

[OCEAN_DATA]
# Name of the ocean circulation model or data source from which velocity data
# have been generated. Supported options include:
#   : ArakawaA: Model outputs on an Arakawa-a grid with CF compliant variable names
#   : FVCOM: The Finite Volume Community Ocean Model
#   : GOTM: The General Ocean Turbulence Model
name = ArakawaA

# Directory containing model output files
data_dir = /data/arakawa_a

# File containing grid info, in particular the adjacency. If not given, this
# data is extracted from the first data file, which will increase run times.
# The full path is required.
grid_metrics_file = %(data_dir)s/grid_metrics.nc

# Data file name stem (e.g. "data_"). Files should be stored at the top level of
# `data_dir'. The model will automatically look for and concatenate across files
# whose names begin with this string (e.g. data_0001.nc, data_0002.nc etc etc).
# If a grid metrics file is given make sure its file name does not match
# data_file_stem, otherwise the model will attempt to read this too.
data_file_stem = data_

# Rounding interval (s) used for rounding datetime objects constructed from the
# model time variable, which may have been written to file with only limited
# numerical precision. e.g. a rounding interval of 3600s will round
# times to the nearest hour.
rounding_interval = 3600

# Constant value for the horizontal eddy diffusivity (units: m^2/s). This value is
# only used with the iterative method `Diff_Const_2D'. See section `NUMERICS'.
horizontal_eddy_diffusivity_constant = 10.0

# Dimension names
# ---------------

# Time dimension
time_dim_name = time

# Depth dimension name
depth_dim_name = depth

# Latitude dimension name
latitude_dim_name = latitude

# Longitude dimension name
longitude_dim_name = longitude

# Variable names
# --------------

# Time (mandatory)
time_var_name = time

# u-velocity component (mandatory)
uo_var_name = uo_foam

# v-velocity component (mandatory)
vo_var_name = vo_foam

# w-velocity component (leave blank or comment if not included)
wo_var_name =

# Sea surface elevation (leave blank or comment if not included)
zos_var_name =

# Temperature (leave blank or comment if not included)
thetao_var_name =

# Salinity (leave blank or comment if not included)
so_var_name =

# Does the mask change with wetting and drying?
has_is_wet = False

# Options for vertical eddy diffusivity
# -------------------------------------

# Method used to calculate the vertical eddy diffusivity (Kz). NB values are case insensitive.
#     : none - Kz cannot be retrieved or calculated from the data set.
#     : file - Kz is read in from file. Users must provide a value for `Kz_var_name`.
Kz_method = none

# Vertical eddy diffusivity variable name (leave blank or comment if `Kz_method` is `none`)
Kz_var_name =

# Options for horizontal eddy diffusivity
# ---------------------------------------

# Method used to calculate the horizontal eddy diffusivity (Ah). NB values are case insensitive.
#     : none        - Ah cannot be retrieved or calculated from the data set. Setting this value
#                     to none does not preclude the use of a constant value for Ah. This is viewed
#                     as being a property of the numerical scheme, and is set in the `NUMERICS`
#                     section of the config.
#     : smagorinsky - Compute Ah using Smagorinsky formulation. Parameters for the Smagorinksy
#                     relationship should be set in the section [SMAGORINSKY].
#     : file        - Ah is read in from file. With this option users must provide a value for `Ah_var_name`.
Ah_method = none

# Horizontal eddy diffusivity variable name (leave blank or comment if `Ah_method` is `none`)
Ah_var_name =

# Interpolation scheme used when computing the vertical eddy diffusivity. At
# present, this option is only supported by GOTM. If set when working with
# FVCOM model outputs, it will have absolutely no effect on the interpolation
# scheme used. See PyLag's documentation for more details.
# linear - Linear interpolation
# cubic_spline = Cubic spline interpolation
# vertical_interpolation_scheme = cubic_spline

[SMAGORINSKY]
# Constant to use with the Smagoriksy estimate of horizontal eddy diffusivities.
constant = 0.2

[NUMERICS]
# The numerical method used to compute changes in particle positions. This option
# controls whether operator splitting is used to combine the contributions of
# advection and diffusion.
# standard          - The iterative method is set through the option `iterative_method'.
#                     The specified iterative method may correspond to cases of pure
#                     advection, pure diffusion or both. In the case of both advection
#                     and diffusion, operator splitting is not performed; rather, the effects
#                     of both processes are computed at the same time using the same time step
#                     and using the specified iterative method. StdNumMethod objects require no
#                     knowledge or the type of iterative method being used; they simply
#                     use it compute the change in a particle's position.
# operator_split_0  - A form of operator splitting in which the advection step is
#                     performed first followed by the diffusion step. The iterative
#                     methods used for the advection and diffusion steps are set through
#                     the options `adv_iterative_method' and `diff_iterative_method'
#                     respectively. The two processes can use different time steps; these
#                     are set through the parameters `time_step_adv' and `time_step_diff'
#                     respectively. `time_step_adv' must be greater than or equal to
#                     `time_step_diff' and `time_step_adv'%`time_step_diff' must equal 0.
#                     An exception will be raised if either of these conditions is not
#                     met.
# operator_split_1  - A form of operator splitting in which a half step for diffusion
#                     is performed first followed by a full advection step followed
#                     by a half diffusion step. The iterative methods used for the advection
#                     and diffusion steps are set through the options `adv_iterative_method'
#                     and `det_iterative_method' respectively. The two processes use different
#                     time steps; these are set through the parameters `time_step_adv' and
#                     `time_step_diff' respectively. `time_step_adv' must be exactly equal to
#                     2 * `time_step_diff'. An exception will be raised if this condition is not
#                     met.
num_method = standard

# The iterative method used to simulate pure advection, pure diffusion or both combined
# (in the absence of operator splitting). For deterministic methods that solve for advection
# only, "Adv" is appended to the beginning of the name. For stochastic methods that solve for
# diffusion only, "Diff" is appended to the beginning of the name. Lastly, for methods that do
# both advection and diffusion, "AdvDiff" is appended to the beginning of the name. The
# dimensionality of each scheme is implicit within the name. 1D schemes compute changes
# in position in e_k only; 2D schemes in e_i and e_j only; and 3D scehemes in e_i, e_j and e_k.
# Adv_RK4_2D - 2D fourth order Runge-Kutta (uses horizontal velocities only)
# Adv_RK4_3D - 3D fourth order Runge-Kutta (uses horizontal and vertical velocities)
# Diff_Const_2D          - 2D method that uses a constant horizontal eddy diffusivity. The constant
#                          is set using the parameter `horizontal_eddy_diffusivity_constant'
#                          in the section OCEAN_DATA.
# Diff_Naive_[1,2]D      - 1- or 2-D Naive scheme; only suitable for homogeneous diffusivity/
#                          diffusivity fields.
# Diff_Euler_1D          - 1D Euler sceheme which includes a pseudo velocity correction term for
#                          inhomogeneous diffusivity fields. The sceheme converges with
#                          delta t and sqrt(delta t) in the weak and strong sense respectively.
# Diff_Visser_1D         - 1D Visser scheme which includes a pseudo velocity correction for
#                          inhomgeneous diffusivity fields. Diffusivities are evaluated at a z
#                          position offset from Z_n. The scheme converges with delta t and
#                          sqrt(delta t) in the weak and strong sense respectively.
# Diff_Milstein_[1,2,3]D - 1-, 2- or 3-D Milstein schemes that converge with delta t in both the
#                          strong and the weak sense.
# AdvDiff_Milstein_3D    - 3D Milstein scheme that converges with delta t in both the strong and the
#                          weak sense.
iterative_method = Adv_RK4_3D

# The iterative method used for advection when operator splitting is used:
adv_iterative_method = Adv_RK4_3D

# The iterative method used for diffusion when operator splitting is used:
diff_iterative_method = Diff_Milstein_3D

# The time step used by iterative methods that deal with advection only (s)
time_step_adv = 100

# The time step used by iterative methods that deal with either diffusion only
# or advection and diffusion combined (s).
time_step_diff = 1

[BOUNDARY_CONDITIONS]

# Horizontal boundary condition calculator
# NB These are specific to the coordinate system (cartesian or geographic) in which
# the model is run.
# None - Run without a horizontal boundary condition calculator. This will only
# work successfully if there are no land boundaries, or if the setup is such that
# horizontal boundary crossings are impossible. If neither of these are true,
# and this option is set to `None', an exception will be raised.
# restoring - Restore the particle to its last known good position.
# reflecting - Apply reflecting horizontal boundary conditions
horiz_bound_cond = restoring

# Vertical boundary condition calculator
# None - Run without a vertical boundary condition calculator. It is okay to set
# this to `None' if it is possible to guarantee a vertical boundadry will not be
# crossed (e.g. if the input data is 2D). An exception will be raised if a vertical
# boundary is crossed and this option is set to `None'.
# reflecting - Apply reflecting vertical boundary conditions
vert_bound_cond = reflecting

[OUTPUT]

# List of environmental variables to be saved to file along with particle data.
# This facility allows one to analyse changes in a given environmental variable
# (e.g. temperature) along a particle's pathline. Values are output along with
# particle position data. For now, support for the following environmental
# variables exists:
#
# ArakawaA
# --------
# thetao - Sea water potential temperature
# so - Salinity
#
# FVCOM
# -----
# thetao - Sea water potential temperature
# so - Salinity
#
# GOTM
# ----
# N/A
#
# Variables should be given as a comma separated list.
#environmental_variables = thetao, so

Step 2: Running PyLag

With PyLag installed and all required inputs generated, you can start a new simulation. PyLag can be run in serial or parallel modes. A simulation can be run in serial using the following commands:

$ python -m pylag.main -c pylag.cfg

where pylag.cfg is the name of the run configuration file.

To launch a parallel run, simply type:

$ mpiexec -np n_proc python -m pylag.parallel.main -c pylag.cfg

where n_proc is the number of processors the run will be distributed over.

Note:

For the time being, the number of processors used must be an exact divisor of the total number of particles. If this is not the case, an error will be raised and the simulation halted.

If the model run completed successfully, a set of outputs will have been saved in the specified output directory.

Step3: Analysing simulation outputs

The pylag.processing sub-package provides a number of tools to assist with analysing simulation outputs. Examples of their usage can be found in the tutorial examples.