Ways to reduce computing time for a long, multi-batch 3D run with OFES (Parcels v3.1.4)?

[imariaguiar]

Question

Hi,

I'm a PhD student studying the AABW downstream of the Vema Channel, using Parcels v3.1.4 with OFES (0.1°, daily U, V, W, T, S). Before I start the production run, I'd like to ask if there is anything I'm not doing that could reduce the computing time.

I want to follow the AABW pathways across the South Atlantic and quantify its transport and transformation. For that, particles have to be integrated for a very long time: 40 years per particle, while the OFES data cover 28 years (1996-2023), so I use time_periodic to loop the fields.

Current setup (without diffusion)

• Particles are released at every OFES grid point with valid v along 31°S (40°W to 38.6°W, below 3250 m), about 50 release points per release date. Each particle carries a fixed Lagrangian transport weight (v x cell area).
• 100 release dates ("batches"), each one a separate ParticleSet, spread over 10 years.
• Kernels: AdvectionRK4_3D, a custom kernel that pushes particles away from the seafloor (it samples a bathymetry field and the UV at the 4 cardinal and 4 diagonal neighbours when the next step would hit the bottom), then SampleUV, SampleTS, CheckDomainBounds, CheckSurfaceExit.
• dt = 3 h, outputdt = 1 day, float32 positions.
• Batches run in parallel with ProcessPoolExecutor (4 workers). Each worker builds its own FieldSet.from_netcdf once (U, V, W, T, S, daily files, default chunking) and reuses it for its batches.
• Machine: Mac with 10 cores and 32 GB RAM, data on an external drive.

kernels_list = [parcels.AdvectionRK4_3D, CardinalDiagonalGridStep,
                SampleUV, SampleTS, CheckDomainBounds, CheckSurfaceExit]

pset = parcels.ParticleSet(fieldset=fieldset, pclass=SampleParticle,
                           lon=lon_points, lat=np.full(n_points, lat_used),
                           depth=lev_points, time=release_time, transport=transport_points)
pset.execute(kernels_list, runtime=timedelta(days=365 * 40), dt=timedelta(hours=3),
             output_file=pset.ParticleFile(name=..., outputdt=timedelta(days=1)))

Timing
From my tests, building the FieldSet takes about 20-25 min per worker, and a batch takes roughly 1h40 for 2 years of tracking. Extrapolating to 100 batches x 40 years gives about a month of computing.

My suspicion (not tested yet) is that the time is dominated by reading the daily fields, not by the number of particles.

What comes next
I'm preparing a version with diffusion: constant Kh plus a vertical diffusivity Kz taken from the OFES output, read every day instead of a time mean. I know this will make it slower. I'm also testing whether I can use the model fields every 3 days instead of every day, by checking how much variance is lost with a linear interpolation.

Questions

1. Is there something I'm not doing that would reduce the time? For example, pre-concatenating the daily files into bigger files or zarr, or setting field_chunksize / deferred_load.
2. Does time_periodic with 28 years of daily fields add overhead, or is there a better way to run 40 years on 28 years of data?
3. Would it be more efficient to release all dates in a single ParticleSet (with an array of release times) so the fields are read once, instead of one ParticleSet per batch? Or is it better to keep separate processes?
4. Each worker builds its own FieldSet. Is there a better way to parallelize on a single machine (MPI, sharing the fields)?
5. Is the cost of my custom seafloor kernel (several Bathymetry and UV samples per step when blocked) likely to matter, and can dt be larger than 3 h with daily fields?
6. Would moving to Parcels v4 help for this kind of run?

I can share the full script if useful. Thank you very much for any suggestion!

[erikvansebille]

Thanks for reaching out @imariaguiar, and for the very extensive and detailed write-up.

My short answer would be: I don't think there are any easy ways to significantly speed up your simulation, except to use indices (in v3) to subset a part of the region. You may then get OutofBounds errors when particles hit the boundaries of your subdomain, though...

If the number of particles per node is less than a few thousands, it could be much fast in Parcels v4, because chunking works better there. Use the ChunkCachedArrays backend then, and hope that the chunking in your OFES files aligns well with the particles spreading you expect. However, I can't guarantee that (we're still actively benchmarking the performance of v4) and I do realise it might be a bit of effort to change your code. I leave it to you to decide if it's worthwhile