There are many different ways in which to use Parcels for research. The flexibility of the parcels code enables this wide range of applicability and allows you to build complex simulations. But this also means that you have to know what you're doing. SO please take some time to learn how to use Parcels, starting with this tutorial.
For a smooth programming experience with Parcels, it is recommended you make a general structure in setting up your simulations. Here we give you an overview of the main components of a Parcels simulation. These components provide you with the basic requirements for your first simulation, while the structure allows you to keep track of more complex simulations later on. A good practice is to separate the different parts of your code into sections:
Variables
they have and what their initial conditions are.These four components are used in the python cell below. Further on in this notebook, we'll focus on each component separately.
from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4
# 1. Setting up the velocity fields in a FieldSet object
fname = 'GlobCurrent_example_data/*.nc'
filenames = {'U': fname, 'V': fname}
variables = {'U': 'eastward_eulerian_current_velocity', 'V': 'northward_eulerian_current_velocity'}
dimensions = {'U': {'lat': 'lat', 'lon': 'lon', 'time': 'time'},
'V': {'lat': 'lat', 'lon': 'lon', 'time': 'time'}}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
# 2. Defining the particles type and initial conditions in a ParticleSet object
pset = ParticleSet(fieldset=fieldset, # the fields on which the particles are advected
pclass=JITParticle, # the type of particles (JITParticle or ScipyParticle)
lon=28, # release longitudes
lat=-33) # release latitudes
# 3. Executing an advection kernel on the given fieldset
output_file = pset.ParticleFile(name="GCParticles.nc", outputdt=3600) # the file name and the time step of the outputs
pset.execute(AdvectionRK4, # the kernel (which defines how particles move)
runtime=86400*6, # the total length of the run
dt=300, # the timestep of the kernel
output_file=output_file)
# 4. Exporting the simulation output to a netcdf file
output_file.export()
WARNING: Casting lon data to np.float32 WARNING: Casting lat data to np.float32 WARNING: Casting depth data to np.float32 INFO: Compiled JITParticleAdvectionRK4 ==> C:\Users\GEBRUI~1\AppData\Local\Temp\parcels-tmp\35821cfe0ca062441b33915339ad4dac_0.dll WARNING: Casting field data to np.float32
When you start making the parcels simulation more complex, it is a good idea to keep these different steps separate to keep a clear overview and find bugs more easily
Parcels provides a framework to simulate the movement of particles within an existing flowfield 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.
The general method to use is FieldSet.from_netcdf
, which requires filenames
, variables
and dimensions
. Each of these is a dictionary, and variables
requires at least a U
and V
, but any other variable
can be added too (e.g. vertical velocity, temperature, mixedlayerdepth, etc). Note also that filenames
can contain wildcards. For example, the GlobCurrent data that can be downloaded using the parcels_get_examples
script (see step 4 of the installation guide) can be read with:
fname = 'GlobCurrent_example_data/*.nc'
filenames = {'U': fname, 'V': fname}
variables = {'U': 'eastward_eulerian_current_velocity', 'V': 'northward_eulerian_current_velocity'}
dimensions = {'U': {'lat': 'lat', 'lon': 'lon', 'time': 'time'}, # In the GlobCurrent data the dimensions are also called 'lon', 'lat' and 'time'
'V': {'lat': 'lat', 'lon': 'lon', 'time': 'time'}}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
FieldSets
:¶Once you have set up the environment with the FieldSet object, you can start defining your particles in a ParticleSet
object. This object requires:
pset = ParticleSet(fieldset=fieldset, # the fields on which the particles are advected
pclass=JITParticle, # the type of particles (JITParticle or ScipyParticle)
lon=28, # release longitude
lat=-33) # release latitude
The different Particle
types available are the JITParticle
and the ScipyParticle
, but it is very easy to create your own particle class which includes other Variables
:
from parcels import Variable
class PressureParticle(JITParticle): # Define a new particle class
p = Variable('p', initial=0) # Variable 'p' with initial value 0.
ParticleSet
:¶For more information on how to implement Particle
types with specific behavior, see the section on writing your own kernels
After defining the flowfield environment with FieldSet
and the particle information with ParticleSet
, we can move on to actually running the parcels simulation by using ParticleSet.execute()
. Running a simulation in parcels actually means executing kernels, little snippets of code that are run for each particle at each timestep. The most basic kernels are advection kernels which calculate the movement of each particle based on the FieldSet
in which the ParticleSet
lives. A few different advection kernels are included in Parcels.
If you want to store the particle data generated in the simulation, you usually first want to 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(AdvectionRK4, # the kernel (which defines how particles move)
runtime=86400*6, # the total length of the run in seconds
dt=300, # the timestep of the kernel
output_file=output_file)
INFO: Compiled JITParticleAdvectionRK4 ==> C:\Users\GEBRUI~1\AppData\Local\Temp\parcels-tmp\a23c56b16e735df652d3b3c4b045eabb_0.dll
One of the most powerful features of Parcels is the ability to write custom Kernels (see e.g. this example to add the vertical movements of an Argo float). You probably want to define these kernels here; after defining your ParticleSet
and before executing them. If your kernels become very large or complex, you might want to store them in another python file and import them into your simulation script.
def WestVel(particle, fieldset, time):
if time > 86400 and time < 2*86400:
uvel = -0.02
particle.lon += uvel * particle.dt
WV_kernel = pset.Kernel(WestVel)
output_file = pset.ParticleFile(name="GC_WestVel.nc", outputdt=3600)
pset.execute(AdvectionRK4 + WV_kernel, # simply add kernels using the + operator
runtime=86400*6, # the total length of the run in seconds
dt=300, # the timestep of the kernel
output_file=output_file)
INFO: Compiled JITParticleAdvectionRK4WestVel ==> C:\Users\GEBRUI~1\AppData\Local\Temp\parcels-tmp\2d3461f1f21e8976d23b053f3bf517e8_0.dll
However, there are some key limitations 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.
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.