GridSearchingError for interior-ocean particles on MOM6 curvilinear grid (ParcelsV4-alpha)

[leilanepassos]

Question

Question

Hi everyone,

I'm trying to run parcels in a MOM6 curvilinear grid using ParcelsV4 alpha (Parcels: 3.1.3.dev2070, Python 3.14, Linux x86_64 (conda)).

I used the script that converts NEMO outputs to SGRID format to create a similar script for MOM6. I tested with North Eastern Pacific domain, 816×342, built with mom6_to_sgrid + FieldSet.from_sgrid_conventions, mesh="spherical"), pset.execute raises GridSearchingError for particles at interior, open-ocean locations, far from land and far from the grid's degenerate cells. One failing particle aborts the whole run.

The static file is land-masked with NaN and I filled by row/column linear interpolation, so all coordinate arrays are NaN-free, and lon is continuous (157–255°E, no dateline jump inside the domain). Even with this smooth fill, Parcels still detects 67 degenerate faces. The simulation was run with 100 particles released on a 10×10 grid over 210–220°E, 40–56°N, AdvectionRK4, dt = 2 min.

Is this the recommended approach for reading a MOM6 curvilinear grid in Parcels? Based on this description, does anything in my grid conversion or SGRID setup look incorrect or potentially problematic?
If not, could you point me toward the most useful things to investigate to understand why GridSearchingError is occurring for particles that appear to be well inside the valid ocean domain?

I can provide the mom6_to_sgrid script, the resulting SGRID files if that would be helpful.

Thanks!
Leilane

Supporting code/error messages

import argparse
import os
import warnings
from datetime import timedelta

import numpy as np
import xarray as xr
import parcels
from mom6_tools import mom6_to_sgrid

from benchmark_utils import Benchmark

warnings.filterwarnings("ignore")

# ---- fixed inputs: identical to the V3 benchmarks -------------------------
COORDS_FILE = "cefi_ocean_static_v4.nc"
U_FILE = "ssu.nep.full.hcast.daily.raw.r20241015.199301-201912.nc"
V_FILE = "ssv.nep.full.hcast.daily.raw.r20241015.199301-201912.nc"

# Al's release region (210-220E / 40-56N) in -180..180 convention
# RELEASE_LON_MIN, RELEASE_LON_MAX = -150.0, -140.0
RELEASE_LON_MIN, RELEASE_LON_MAX = 210.0, 220.0
RELEASE_LAT_MIN, RELEASE_LAT_MAX = 40.0, 56.0


def main():
    p = argparse.ArgumentParser(description="Passive benchmark of Parcels V4 (advection-only)")
    p.add_argument("--nx", type=int, default=10, help="release grid points in longitude")
    p.add_argument("--ny", type=int, default=10, help="release grid points in latitude")
    p.add_argument("--sim-days", type=float, default=10.0, help="integration length in days")
    p.add_argument("--dt", type=float, default=2.0, help="advection timestep, MINUTES (default 2)")
    p.add_argument("--outputdt", type=float, default=24.0, help="output interval, HOURS (default 24)")
    p.add_argument("--name", default="v4_passive_metrics", help="benchmark name")
    p.add_argument("--out", default=None, help="output parquet path")
    p.add_argument("--metrics-out", default=None, help="metrics JSON path")
    args = p.parse_args()

    out_path = args.out or f"output/{args.name}.parquet"
    metrics_out = args.metrics_out or f"{args.name}.json"
    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)

    n_days_forcing = int(np.ceil(args.sim_days)) + 2   

    print("=== PASSIVE benchmark of Parcels V4 native grid, advection-only ===")
    print(f"Output:  {out_path}")
    print(f"Metrics: {metrics_out}")
    print(f"Particles: {args.nx}x{args.ny} = {args.nx * args.ny}, "
          f"runtime {args.sim_days} d, dt {args.dt} min, "
          f"forcing days loaded: {n_days_forcing}")

    bench = Benchmark(args.name)

    # ---- fieldset -------------------------------------------------------
    with bench.phase("fieldset_build"):
        u_ds = xr.open_dataset(U_FILE, decode_timedelta=False)
        v_ds = xr.open_dataset(V_FILE, decode_timedelta=False)
        dg = xr.open_dataset(COORDS_FILE)

        u = u_ds["ssu"].sortby("time").isel(time=slice(0, n_days_forcing))
        v = v_ds["ssv"].sortby("time").isel(time=slice(0, n_days_forcing))

        sgrid_ds = mom6_to_sgrid(fields={"U": u, "V": v}, coords=dg)
        fieldset = parcels.FieldSet.from_sgrid_conventions(sgrid_ds, mesh="spherical")
    print("FieldSet built.")

    # ---- seed particles -------------------------------------------------
    with bench.phase("seeding"):
        lons = np.linspace(RELEASE_LON_MIN, RELEASE_LON_MAX, args.nx)
        lats = np.linspace(RELEASE_LAT_MIN, RELEASE_LAT_MAX, args.ny)
        lon_grid, lat_grid = np.meshgrid(lons, lats)
        pset = parcels.ParticleSet(
            fieldset=fieldset,
            pclass=parcels.Particle,
            lon=lon_grid.ravel(),
            lat=lat_grid.ravel(),
        )
    n0 = len(pset)
    print(f"Released {n0} particles ({args.nx}x{args.ny} grid) in "
          f"lon [{RELEASE_LON_MIN}, {RELEASE_LON_MAX}], "
          f"lat [{RELEASE_LAT_MIN}, {RELEASE_LAT_MAX}]")

    # ---- run ------------------------------------------------------------
    with bench.phase("advection"):
        output_file = parcels.ParticleFile(
            path=out_path,
            outputdt=timedelta(hours=args.outputdt),
        )
        pset.execute(
            parcels.kernels.AdvectionRK4,
            runtime=timedelta(days=args.sim_days),
            dt=timedelta(minutes=args.dt),
            output_file=output_file,
        )
    print(f"Simulation complete. Output written to: {out_path}")

    n_steps = args.sim_days * 24 * 60 / args.dt
    bench.finalize(n_particles_start=n0, n_particles_end=len(pset), n_steps=n_steps)
    bench.report()
    bench.save(metrics_out)
    print(f"\nMetrics written to {metrics_out}")


if __name__ == "__main__":
    main()

Integration time: 1993-01-01T12:00:00   0%|                                                                                                         | [00:00<?, ?it/s]/home/miniconda3/envs/parcels-v4-alpha/lib/python3.14/site-packages/parcels/_core/spatialhash.py:592: RuntimeWarning: invalid value encountered in cast
  xq = np.clip((xn * bitwidth).astype(np.uint32), 0, bitwidth)
/home/miniconda3/envs/parcels-v4-alpha/lib/python3.14/site-packages/parcels/_core/spatialhash.py:593: RuntimeWarning: invalid value encountered in cast
  yq = np.clip((yn * bitwidth).astype(np.uint32), 0, bitwidth)
/home/miniconda3/envs/parcels-v4-alpha/lib/python3.14/site-packages/parcels/_core/spatialhash.py:594: RuntimeWarning: invalid value encountered in cast
  zq = np.clip((zn * bitwidth).astype(np.uint32), 0, bitwidth)
FAIL  west of dateline       (170.0, 50.0) -> GridSearchingError