The flexibility of Parcels allows for a wide range of applications and to build complex simulations. In order to help structuring your code, this tutorial describes the structure that a Parcels script uses.
Code that uses Parcels is generally build up from four different components:
Variables
can be added to the particles (e.g. temperature, to keep track of the temperature that particles experience).We discuss each component in more detail below.
from IPython.display import SVG
# Illustration of Parcels code structure
SVG(filename='parcels_user_diagram.svg')
Parcels provides a framework to simulate the movement of particles within an existing flow field environment. To start a parcels simulation we must define this environment with the FieldSet
class. The minimal requirements for this Fieldset are that it must contain the 'U'
and 'V'
fields: the 2D hydrodynamic data that will move the particles in a horizontal direction. Additionally, it can contain e.g. a temperature or vertical flow field.
A fieldset can be loaded with FieldSet.from_netcdf
, if the model output with the fields is written in NetCDF files. This function requires filenames
, variables
and dimensions
.
In this example, we only load the 'U'
and 'V'
fields, which represent the zonal and meridional flow velocity. First, fname
points to the location of the model output.
fname = 'GlobCurrent_example_data/*.nc'
filenames = {'U': fname, 'V': fname}
Second, we have to specify the 'U'
and 'V'
variable names, as given by the NetCDF files.
variables = {'U': 'eastward_eulerian_current_velocity', 'V': 'northward_eulerian_current_velocity'}
Third, we specify the names of the variable dimensions, as given by the NetCDF files.
# In the GlobCurrent data the dimensions are also called 'lon', 'lat' and 'time
dimensions = {'U': {'lat': 'lat', 'lon': 'lon', 'time': 'time'},
'V': {'lat': 'lat', 'lon': 'lon', 'time': 'time'}}
Finally, we load the fieldset.
from parcels import FieldSet
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
FieldSets
:¶Once the environment has a FieldSet
object, you can start defining your particles in a ParticleSet
object. This object requires:
from parcels import Variable, JITParticle, ParticleSet
# Define a new particle class
class AgeParticle(JITParticle): # It is a JIT particle
age = Variable('age',
initial=0) # Variable 'age' is added with initial value 0.
pset = ParticleSet(fieldset=fieldset, # the fields that the particleset uses
pclass=AgeParticle, # define the type of particle
lon=29, # release longitude
lat=-33) # release latitude
ParticleSet
:¶For more information on how to implement Particle
types with specific behaviour, see the section on writing your own kernels
Kernels are little snippets of code, which are applied to every particle in the ParticleSet
, for every time step during a simulation.
Basic kernels are included in Parcels, among which several types of advection kernels. AdvectionRK4
is the main advection kernel for two-dimensional advection, which is also used in this example.
One can also write custom kernels, to add certain types of 'behaviour' to the particles.
# load basic advection kernels
from parcels import AdvectionRK4
# Create a custom kernel which displaces each particle southward
def NorthVel(particle, fieldset, time):
if time > 10*86400 and time < 10.2*86400:
vvel = -1e-4
particle.lat += vvel * particle.dt
# Create a custom kernel which keeps track of the particle age (minutes)
def Age(particle, fieldset, time):
particle.age += particle.dt / 3600
# define all kernels to be executed on particles
kernels = pset.Kernel(Age) + pset.Kernel(NorthVel) + AdvectionRK4
Some key limitations exist to the Kernels that everyone who wants to write their own should be aware of:
Every Kernel must be a function with the following (and only those) arguments: (particle, fieldset, time)
In order to run successfully in JIT mode, Kernel definitions can only contain the following types of commands:
+
, -
, *
, /
, **
) and assignments (=
).<
, ==
, !=
, >
, &
, |
). Note that you can use a statement like particle.lon != particle.lon
to check if particle.lon
is NaN (since math.nan != math.nan
).if
and while
loops, as well as break
statements. Note that for
-loops are not supported in JIT mode.Field
from the FieldSet
at a [time, depth, lat, lon]
point, using square brackets notation.For example, to interpolate the zonal velocity (U) field at the particle location, use the following statement:
value = fieldset.U[time, particle.depth, particle.lat, particle.lon]
Functions from the math
standard library and from the custom ParcelsRandom
library at parcels.rng
Simple print
statements, such as:
print("Some print")
print(particle.lon)
print("particle id: %d" % particle.id)
print("lon: %f, lat: %f" % (particle.lon, particle.lat))
Although note that these print
statements are not shown in Jupyter notebooks in JIT mode, see this long-standing Issue.
Local variables can be used in Kernels, and these variables will be accessible in all concatenated Kernels. Note that these local variables are not shared between particles, and also not between time steps.
Note that one has to be careful with writing kernels for vector fields on Curvilinear grids. While Parcels automatically rotates the U and V field when necessary, this is not the case for for example wind data. In that case, a custom rotation function will have to be written.
The final part executes the simulation, given the ParticleSet
, FieldSet
and Kernels
, that have been defined in the previous steps. If you like to store the particle data generated in the simulation, you define the ParticleFile
to which the output of the kernel execution will be written. Then, on the ParticleSet
you have defined, you can use the method ParticleSet.execute()
which requires the following arguments:
runtime
defining how long the execution loop runs. Alternatively, you may define the endtime
at which the execution loop stops.dt
at which to execute the kernels.ParticleFile
object to write the output to.output_file = pset.ParticleFile(name="GCParticles.nc", outputdt=3600) # the file name and the time step of the outputs
pset.execute(kernels, # the kernel (which defines how particles move)
runtime=86400*24, # the total length of the run in seconds
dt=300, # the timestep of the kernel in seconds
output_file=output_file)
INFO: Compiled ArrayAgeParticleAgeNorthVelAdvectionRK4 ==> /var/folders/s5/gxtkk3c12yqgd7hkt1b_s0vr0000gq/T/parcels-503/lib0e6e33cfac25b37419e8147f88e31705_0.so
While executing the ParticleSet
, parcels stores the data in npy files in an output folder. To take all the data and store them in a netcdf file, you can use ParticleFile.export()
if you want to keep the folder with npy files; or ParticleFile.close()
if you only want to keep the netcdf file:
output_file.export()
output_file.close()
After running your simulation, you probably want to analyse the output. Although there is some simple plotting functionality built into Parcels, we recommend you write your own code to analyse your specific output and you can probably separate the analysis from the simulation.