More xgrid generated datasets by erikvansebille · Pull Request #2135 · Parcels-code/Parcels

parcels/_datasets/structured/generated.py

-Original file line number
+Diff line change
@@ -1,3 +1,5 @@
+import math
 import numpy as np
 import xarray as xr
@@ Expand All @@
  "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}),
  },
  )
+def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes.
+ lon = np.linspace(0, 60, xdim, dtype=np.float32)
+ lat = np.linspace(0, 60, ydim, dtype=np.float32)
+ x0 = 30.0 # Define the origin to be the centre of the Field.
+ y0 = 30.0
+ U = np.zeros((2, 1, ydim, xdim), dtype=np.float32)
+ V = np.zeros((2, 1, ydim, xdim), dtype=np.float32)
+ omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day.
+ for i in range(lon.size):
+ for j in range(lat.size):
+ r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2)
+ assert r >= 0.0
+ assert r <= np.sqrt(x0**2 + y0**2)
+ theta = np.arctan2((lat[j] - y0), (lon[i] - x0))
+ assert abs(theta) <= np.pi
+ U[:, :, j, i] = r * np.sin(theta) * omega
+ V[:, :, j, i] = -r * np.cos(theta) * omega
+ return xr.Dataset(
+ {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)},
+ coords={
+ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}),
+ "depth": (["depth"], np.array([0.0]), {"axis": "Z"}),
+ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}),
+ "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}),
+ "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}),
+ "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}),
+ "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}),
+ "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}),
+ },
+ )
+def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1
+ """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow."""
+ f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants
+ lon = np.linspace(0, 25000, xdim, dtype=np.float32)
+ lat = np.linspace(0, 25000, ydim, dtype=np.float32)
+ time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m"))
+ U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32)
+ V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32)
+ for t in range(len(time)):
+ U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s")))
+ V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s")))
+ return xr.Dataset(
+ {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)},
+ coords={
+ "time": (["time"], time, {"axis": "T"}),
+ "depth": (["depth"], np.array([0.0]), {"axis": "Z"}),
+ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}),
+ "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}),
+ "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}),
+ "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}),
+ "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}),
+ "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}),
+ },
+ attrs={
+ "u_0": u_0,
+ "u_g": u_g,
+ "f": f,
+ },
+ )
+def decaying_moving_eddy_dataset(xdim=2, ydim=2):
+ """Simulate an ocean that accelerates subject to Coriolis force
+ and dissipative effects, upon which a geostrophic current is
+ superimposed.
+ The original test description can be found in: N. Fabbroni, 2009,
+ Numerical Simulation of Passive tracers dispersion in the sea,
+ Ph.D. dissertation, University of Bologna
+ http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf
+ """
+ u_g = 0.04 # Geostrophic current
+ u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0
+ gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity.
+ gamma_g = 1.0 / (28.9 * 86400)
+ f = 1.0e-4 # Coriolis parameter.
+ time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m"))
+ lon = np.linspace(0, 20000, xdim, dtype=np.float32)
+ lat = np.linspace(5000, 12000, ydim, dtype=np.float32)
+ U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32)
+ V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32)
+ for t in range(time.size):
+ t_float = time[t] / np.timedelta64(1, "s")
+ U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float)
+ V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float)
+ return xr.Dataset(
+ {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)},
+ coords={
+ "time": (["time"], time, {"axis": "T"}),
+ "depth": (["depth"], np.array([0.0]), {"axis": "Z"}),
+ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}),
+ "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}),
+ "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}),
+ "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}),
+ "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}),
+ "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}),
+ },
+ attrs={
+ "u_0": u_0,
+ "u_g": u_g,
+ "f": f,
+ "gamma": gamma,
+ "gamma_g": gamma_g,
+ },
+ )
+def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"):
+ """Construct a fieldset encapsulating the flow field around an idealised peninsula.
+ Parameters
+ ----------
+ xdim :
+ Horizontal dimension of the generated fieldset
+ ydim :
+ Vertical dimension of the generated fieldset
+ mesh : str
+ String indicating the type of mesh coordinates and
+ units used during velocity interpolation:
+. spherical: Lat and lon in degree, with a
+ correction for zonal velocity U near the poles.
+. flat (default): No conversion, lat/lon are assumed to be in m.
+ grid_type :
+ Option whether grid is either Arakawa A (default) or C
+ The original test description can be found in Fig. 2.2.3 in:
+ North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of
+ recommended practices for modelling physical - biological
+ interactions during fish early life.
+ ICES Cooperative Research Report No. 295. 111 pp.
+ http://archimer.ifremer.fr/doc/00157/26792/24888.pdf
+ """
+ domainsizeX, domainsizeY = (1.0e5, 5.0e4)
+ La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32)
+ Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32)
+ u0 = 1
+ x0 = domainsizeX / 2
+ R = 0.32 * domainsizeX / 2
+ # Create the fields
+ P = np.zeros((ydim, xdim), dtype=np.float32)
+ U = np.zeros_like(P)
+ V = np.zeros_like(P)
+ x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy")
+ P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y
+ # Set land points to zero
+ landpoints = P >= 0.0
+ P[landpoints] = 0.0
+ if grid_type == "A":
+ U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2)
+ V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2)
+ U[landpoints] = 0.0
+ V[landpoints] = 0.0
+ Udims = ["YC", "XG"]
+ Vdims = ["YG", "XC"]
+ elif grid_type == "C":
+ U = np.zeros(P.shape)
+ V = np.zeros(P.shape)
+ V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0])
+ U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0])
+ Udims = ["YG", "XG"]
+ Vdims = ["YG", "XG"]
+ else:
+ raise RuntimeError(f"Grid_type {grid_type} is not a valid option")
+ # Convert from m to lat/lon for spherical meshes
+ lon = La / 1852.0 / 60.0 if mesh == "spherical" else La
+ lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa
+ return xr.Dataset(
+ {
+ "U": (Udims, U),
+ "V": (Vdims, V),
+ "P": (["YG", "XG"], P),
+ },
+ coords={
+ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}),
+ "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}),
+ "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}),
+ "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}),
+ "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}),
+ "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}),
+ },
+ )
+def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"):
+ """Simulate a periodic current along a western boundary, with significantly
+ larger velocities along the western edge than the rest of the region
+ The original test description can be found in: N. Fabbroni, 2009,
+ Numerical Simulation of Passive tracers dispersion in the sea,
+ Ph.D. dissertation, University of Bologna
+ http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf
+ """
+ a = b = 10000 * 1e3
+ scalefac = 0.05 # to scale for physically meaningful velocities
+ dx, dy = a / xdim, b / ydim
+ # Coordinates of the test fieldset (on A-grid in deg)
+ lon = np.linspace(0, a, xdim, dtype=np.float32)
+ lat = np.linspace(0, b, ydim, dtype=np.float32)
+ # Define arrays U (zonal), V (meridional) and P (sea surface height)
+ U = np.zeros((lat.size, lon.size), dtype=np.float32)
+ V = np.zeros((lat.size, lon.size), dtype=np.float32)
+ P = np.zeros((lat.size, lon.size), dtype=np.float32)
+ beta = 2e-11
+ r = 1 / (11.6 * 86400)
+ es = r / (beta * a)
+ for j in range(lat.size):
+ for i in range(lon.size):
+ xi = lon[i] / a
+ yi = lat[j] / b
+ P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac
+ if grid_type == "A":
+ U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac
+ V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac
+ if grid_type == "C":
+ V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a
+ U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b
+ Udims = ["YC", "XG"]
+ Vdims = ["YG", "XC"]
+ else:
+ Udims = ["YG", "XG"]
+ Vdims = ["YG", "XG"]
+ return xr.Dataset(
+ {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)},
+ coords={
+ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}),
+ "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}),
+ "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}),
+ "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}),
+ "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}),
+ "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}),
+ },
+ )

parcels/application_kernels/advection.py

-Original file line number
+Diff line change
@@ Expand Up @@
  Time-step dt is halved if error is larger than fieldset.RK45_tol,
  and doubled if error is smaller than 1/10th of tolerance.
  """
- dt = min(
- particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt
- ) # TODO: improve API for converting dt to seconds
+ dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds
+ if dt > fieldset.RK45_max_dt:
+ dt = fieldset.RK45_max_dt
+ particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s")
+ if dt < fieldset.RK45_min_dt:
+ particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s")
+ return StatusCode.Repeat
+ particle.dt = particle.next_dt
  c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]
  A = [
  [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0],
@@ Expand DownExpand Up @@
  if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)):
  particle.dlon += lon_4th
  particle.dlat += lat_4th
- if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)):
+ if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)):
  particle.next_dt *= 2
  else:
  particle.next_dt /= 2
@@ Expand Down @@

parcels/application_kernels/interpolation.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -40,7 +40,10 @@ def XBiLinear( @@
  zi, _ = position["Z"]
  data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2]
- data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :]
+ if tau > 0:
+ data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :]
+ else:
+ data = data[ti, :, :]
  return (
  (1 - xsi) * (1 - eta) * data[0, 0]
@@ Expand Down @@

parcels/fieldset.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -164,7 +164,8 @@ def add_constant(self, name, value): @@
  """
  if name in self.constants:
  raise ValueError(f"FieldSet already has a constant with name '{name}'")
+ if not isinstance(value, (float, np.floating, int, np.integer)):
+ raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}")
  self.constants[name] = np.float32(value)
  @property
@@ Expand Down @@

parcels/kernel.py

Original file line numberDiff line numberDiff line change
@@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime):
 res_tmp = f(p, self._fieldset, p.time_nextloop)
 if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode
 res = res_tmp
 if res == StatusCode.StopExecution:
 if res in [StatusCode.StopExecution, StatusCode.Repeat]:
 break

 if res is None:
Expand Down

parcels/xgrid.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -271,7 +271,10 @@ def _gtype(self): @@
  def search(self, z, y, x, ei=None):
  ds = self.xgcm_grid._ds
- zi, zeta = _search_1d_array(ds.depth.values, z)
+ if "Z" in self.axes:
+ zi, zeta = _search_1d_array(ds.depth.values, z)
+ else:
+ zi, zeta = 0, 0.0
  if zi == -1:
  if zeta < 0:
  raise FieldOutOfBoundError(
@@ Expand Down @@

blocks\n(function() {\n function addCopyButtons() {\n document.querySelectorAll('pre code').forEach(function(codeBlock) {\n if (codeBlock.parentElement.hasAttribute('data-copy-added')) return;\n codeBlock.parentElement.setAttribute('data-copy-added', 'true');\n \n var btn = document.createElement('button');\n btn.textContent = 'Copy';\n btn.style.cssText = 'position:absolute;top:4px;right:4px;padding:2px 8px;font-size:11px;background:#4ecdc4;border:none;border-radius:4px;color:#1a1a2e;cursor:pointer;opacity:0.7;transition:opacity 0.2s;';\n btn.onmouseover = function() { this.style.opacity = '1'; };\n btn.onmouseout = function() { this.style.opacity = '0.7'; };\n btn.onclick = function() {\n navigator.clipboard.writeText(codeBlock.textContent).then(function() {\n btn.textContent = 'Copied!';\n setTimeout(function() { btn.textContent = 'Copy'; }, 1500);\n });\n };\n codeBlock.parentElement.style.position = 'relative';\n codeBlock.parentElement.appendChild(btn);\n });\n }\n \n addCopyButtons();\n \n // Re-run on dynamic content\n var observer = new MutationObserver(addCopyButtons);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Add Copy Buttons to Code Blocks"); } } catch(__e) { console.warn('[Userscript:Add Copy Buttons to Code Blocks]', __e); } })(); (function(){ try { var __m = "github.com"; var __re = new RegExp('^' + "github\\.com" + ' Skip to contentNavigation Menu Sign inAppearance settings AI CODE CREATIONGitHub CopilotWrite better code with AI GitHub Copilot appDirect agents from issue to merge MCP RegistryIntegrate external tools DEVELOPER WORKFLOWSActionsAutomate any workflow CodespacesInstant dev environments IssuesPlan and track work Code ReviewManage code changes Code QualityEnforce quality at merge APPLICATION SECURITYGitHub Advanced SecurityFind and fix vulnerabilities Code securitySecure your code as you build Secret protectionStop leaks before they start EXPLOREWhy GitHub Documentation Blog Changelog Marketplace View all features BY COMPANY SIZEEnterprises Small and medium teams Startups Nonprofits BY USE CASEApp Modernization DevSecOps DevOps CI/CD View all use cases BY INDUSTRYHealthcare Financial services Manufacturing Government View all industries View all solutions EXPLORE BY TOPICAI Software Development DevOps Security View all topics EXPLORE BY TYPECustomer stories Events & webinars Ebooks & reports Business insights GitHub Skills SUPPORT & SERVICESDocumentation Customer support Community forum Trust center Partners View all resources COMMUNITYGitHub SponsorsFund open source developers PROGRAMSSecurity Lab Maintainer Community GitHub Stars Archive Program REPOSITORIESTopics Trending Collections ENTERPRISE SOLUTIONSEnterprise platformAI-powered developer platform AVAILABLE ADD-ONSGitHub Advanced SecurityEnterprise-grade security features Copilot for BusinessEnterprise-grade AI features Premium SupportEnterprise-grade 24/7 support Pricing Sign in Sign upAppearance settings You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert Uh oh! There was an error while loading. Please reload this page. Parcels-code /ParcelsPublic Notifications You must be signed in to change notification settings Fork 188 Star 360 Code Issues109 Pull requests13 Discussions Actions Projects Wiki Security and quality0 Insights Additional navigation options Code Issues Pull requests Discussions Actions Projects Wiki Security and quality Insights More xgrid generated datasets#2135 New issue Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community. Sign up for GitHub By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails. Already on GitHub? Sign in to your account Merged erikvansebille merged 16 commits into v4-dev from more_xgrid_generated_datasetsAug 6, 2025 Conversation 13 Commits 16 Checks 7 Files changed Uh oh! There was an error while loading. Please reload this page. Merged More xgrid generated datasets #2135 Show file tree Hide file tree Changes from all commits Commits Show all changes 16 commits Select commit Hold shift + click to select a range b8177fd Adding radial rotation flow erikvansebille Jul 29, 20253dfd8ac Adding radial rotation unit test erikvansebille Jul 31, 20256e192f8 Adding extra challenge of later particle starttime to test_radialrota… erikvansebille Aug 1, 2025c8b8bfd Refactoring moving_eddy test to generated dataset erikvansebille Aug 1, 202592d2a45 Making kernel loop work again for RK45 erikvansebille Aug 1, 2025c3bb197 Adding Stommel gyre example and test erikvansebille Aug 4, 20254dd3340 Adding test for multple variables setting on Particle erikvansebille Aug 4, 2025a0c02b4 Adding unit test for time-varying (decaying) moving eddy erikvansebille Aug 6, 2025f734195 Adding peninsula dataset to generated and advection tests erikvansebille Aug 6, 202539d1d63 Using np.testing.assert_allclose throughout advection tests erikvansebille Aug 6, 2025cc720ce Update parcels/xgrid.py erikvansebille Aug 6, 2025c58d7f3 Fixing review suggestion erikvansebille Aug 6, 20250287d7c Throwing error when add_constant is not called with float or int erikvansebille Aug 6, 2025568f133 Removing unit test that is now part of #2134 erikvansebille Aug 6, 20252218e8f Update parcels/fieldset.py erikvansebille Aug 6, 20254810e29 Merge branch 'v4-dev' into more_xgrid_generated_datasets erikvansebille Aug 6, 2025 File filter Filter by extension Filter by extension .py (8) All 1 file type selected Viewed files Clear filters Conversations Failed to load comments. Loading Uh oh! There was an error while loading. Please reload this page. Jump to Jump to file Failed to load files. Loading Uh oh! There was an error while loading. Please reload this page. Diff view Unified Split Hide whitespace Diff view Unified Split Hide whitespace parcels/_datasets/structured/generated.py generated.py parcels/application_kernels/advection.py advection.py parcels/application_kernels/interpolation.py interpolation.py parcels/fieldset.py fieldset.py parcels/kernel.py kernel.py parcels/xgrid.py xgrid.py tests/v4/test_advection.py test_advection.py tests/v4/test_diffusion.py test_diffusion.py There are no files selected for viewing 266 changes: 266 additions & 0 deletions 266 parcels/_datasets/structured/generated.py Show comments View file Open in desktop Original file line number Diff line number Diff line change @@ -1,3 +1,5 @@ import math import numpy as np import xarray as xr Expand All @@ -18,3 +20,267 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh_type="spherical"): "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes. lon = np.linspace(0, 60, xdim, dtype=np.float32) lat = np.linspace(0, 60, ydim, dtype=np.float32) x0 = 30.0 # Define the origin to be the centre of the Field. y0 = 30.0 U = np.zeros((2, 1, ydim, xdim), dtype=np.float32) V = np.zeros((2, 1, ydim, xdim), dtype=np.float32) omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day. for i in range(lon.size): for j in range(lat.size): r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2) assert r >= 0.0 assert r <= np.sqrt(x0**2 + y0**2) theta = np.arctan2((lat[j] - y0), (lon[i] - x0)) assert abs(theta) <= np.pi U[:, :, j, i] = r * np.sin(theta) * omega V[:, :, j, i] = -r * np.cos(theta) * omega return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1 Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow.""" f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants lon = np.linspace(0, 25000, xdim, dtype=np.float32) lat = np.linspace(0, 25000, ydim, dtype=np.float32) time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m")) U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) for t in range(len(time)): U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s"))) V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s"))) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, }, ) def decaying_moving_eddy_dataset(xdim=2, ydim=2): """Simulate an ocean that accelerates subject to Coriolis force and dissipative effects, upon which a geostrophic current is superimposed. The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ u_g = 0.04 # Geostrophic current u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0 gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity. gamma_g = 1.0 / (28.9 * 86400) f = 1.0e-4 # Coriolis parameter. time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m")) lon = np.linspace(0, 20000, xdim, dtype=np.float32) lat = np.linspace(5000, 12000, ydim, dtype=np.float32) U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) for t in range(time.size): t_float = time[t] / np.timedelta64(1, "s") U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float) V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, "gamma": gamma, "gamma_g": gamma_g, }, ) def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): """Construct a fieldset encapsulating the flow field around an idealised peninsula. Parameters ---------- xdim : Horizontal dimension of the generated fieldset ydim : Vertical dimension of the generated fieldset mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. 2. flat (default): No conversion, lat/lon are assumed to be in m. grid_type : Option whether grid is either Arakawa A (default) or C The original test description can be found in Fig. 2.2.3 in: North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of recommended practices for modelling physical - biological interactions during fish early life. ICES Cooperative Research Report No. 295. 111 pp. http://archimer.ifremer.fr/doc/00157/26792/24888.pdf """ domainsizeX, domainsizeY = (1.0e5, 5.0e4) La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32) Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32) u0 = 1 x0 = domainsizeX / 2 R = 0.32 * domainsizeX / 2 # Create the fields P = np.zeros((ydim, xdim), dtype=np.float32) U = np.zeros_like(P) V = np.zeros_like(P) x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy") P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y # Set land points to zero landpoints = P >= 0.0 P[landpoints] = 0.0 if grid_type == "A": U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2) V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2) U[landpoints] = 0.0 V[landpoints] = 0.0 Udims = ["YC", "XG"] Vdims = ["YG", "XC"] elif grid_type == "C": U = np.zeros(P.shape) V = np.zeros(P.shape) V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0]) U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0]) Udims = ["YG", "XG"] Vdims = ["YG", "XG"] else: raise RuntimeError(f"Grid_type {grid_type} is not a valid option") # Convert from m to lat/lon for spherical meshes lon = La / 1852.0 / 60.0 if mesh == "spherical" else La lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa return xr.Dataset( { "U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): """Simulate a periodic current along a western boundary, with significantly larger velocities along the western edge than the rest of the region The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ a = b = 10000 * 1e3 scalefac = 0.05 # to scale for physically meaningful velocities dx, dy = a / xdim, b / ydim # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(0, a, xdim, dtype=np.float32) lat = np.linspace(0, b, ydim, dtype=np.float32) # Define arrays U (zonal), V (meridional) and P (sea surface height) U = np.zeros((lat.size, lon.size), dtype=np.float32) V = np.zeros((lat.size, lon.size), dtype=np.float32) P = np.zeros((lat.size, lon.size), dtype=np.float32) beta = 2e-11 r = 1 / (11.6 * 86400) es = r / (beta * a) for j in range(lat.size): for i in range(lon.size): xi = lon[i] / a yi = lat[j] / b P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "A": U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "C": V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b Udims = ["YC", "XG"] Vdims = ["YG", "XC"] else: Udims = ["YG", "XG"] Vdims = ["YG", "XG"] return xr.Dataset( {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)}, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) 14 changes: 10 additions & 4 deletions 14 parcels/application_kernels/advection.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -115,9 +115,15 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ dt = min( particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt ) # TODO: improve API for converting dt to seconds dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds if dt > fieldset.RK45_max_dt: dt = fieldset.RK45_max_dt particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s") if dt < fieldset.RK45_min_dt: particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s") return StatusCode.Repeat particle.dt = particle.next_dt Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0], Expand DownExpand Up @@ -162,7 +168,7 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)): particle.dlon += lon_4th particle.dlat += lat_4th if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): particle.next_dt *= 2 else: particle.next_dt /= 2 Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/application_kernels/interpolation.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -40,7 +40,10 @@ def XBiLinear( zi, _ = position["Z"] data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2] data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] if tau > 0: data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] else: data = data[ti, :, :] return ( (1 - xsi) * (1 - eta) * data[0, 0] Expand Down 3 changes: 2 additions & 1 deletion 3 parcels/fieldset.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -164,7 +164,8 @@ def add_constant(self, name, value): """ if name in self.constants: raise ValueError(f"FieldSet already has a constant with name '{name}'") if not isinstance(value, (float, np.floating, int, np.integer)): raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}") self.constants[name] = np.float32(value) @property Expand Down 2 changes: 1 addition & 1 deletion 2 parcels/kernel.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime): res_tmp = f(p, self._fieldset, p.time_nextloop) if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode res = res_tmp if res == StatusCode.StopExecution: if res in [StatusCode.StopExecution, StatusCode.Repeat]: break if res is None: Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/xgrid.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -271,7 +271,10 @@ def _gtype(self): def search(self, z, y, x, ei=None): ds = self.xgcm_grid._ds zi, zeta = _search_1d_array(ds.depth.values, z) if "Z" in self.axes: zi, zeta = _search_1d_array(ds.depth.values, z) else: zi, zeta = 0, 0.0 if zi == -1: if zeta < 0: raise FieldOutOfBoundError( Expand Down Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later. Footer © 2026 GitHub, Inc. Footer navigation Terms Privacy Security Status Community Docs Contact You can’t perform that action at this time. , 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Force GitHub README to respect dark mode\n(function() {\n var style = document.createElement('style');\n style.textContent = '\n .markdown-body {\n color-scheme: dark light;\n }\n .markdown-body pre { background: #161b22 !important; }\n .markdown-body code { background: rgba(110, 118, 129, 0.4) !important; }\n .markdown-body table th, .markdown-body table td { border-color: #30363d !important; }\n .markdown-body img { background: #0d1117; }\n .markdown-body blockquote { border-left-color: #8b949e; }\n .markdown-body hr { border-color: #30363d; }\n ';\n document.head.appendChild(style);\n})();", "GitHub Dark Mode README Fix"); } } catch(__e) { console.warn('[Userscript:GitHub Dark Mode README Fix]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + ' Skip to contentNavigation Menu Sign inAppearance settings AI CODE CREATIONGitHub CopilotWrite better code with AI GitHub Copilot appDirect agents from issue to merge MCP RegistryIntegrate external tools DEVELOPER WORKFLOWSActionsAutomate any workflow CodespacesInstant dev environments IssuesPlan and track work Code ReviewManage code changes Code QualityEnforce quality at merge APPLICATION SECURITYGitHub Advanced SecurityFind and fix vulnerabilities Code securitySecure your code as you build Secret protectionStop leaks before they start EXPLOREWhy GitHub Documentation Blog Changelog Marketplace View all features BY COMPANY SIZEEnterprises Small and medium teams Startups Nonprofits BY USE CASEApp Modernization DevSecOps DevOps CI/CD View all use cases BY INDUSTRYHealthcare Financial services Manufacturing Government View all industries View all solutions EXPLORE BY TOPICAI Software Development DevOps Security View all topics EXPLORE BY TYPECustomer stories Events & webinars Ebooks & reports Business insights GitHub Skills SUPPORT & SERVICESDocumentation Customer support Community forum Trust center Partners View all resources COMMUNITYGitHub SponsorsFund open source developers PROGRAMSSecurity Lab Maintainer Community GitHub Stars Archive Program REPOSITORIESTopics Trending Collections ENTERPRISE SOLUTIONSEnterprise platformAI-powered developer platform AVAILABLE ADD-ONSGitHub Advanced SecurityEnterprise-grade security features Copilot for BusinessEnterprise-grade AI features Premium SupportEnterprise-grade 24/7 support Pricing Sign in Sign upAppearance settings You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert Uh oh! There was an error while loading. Please reload this page. Parcels-code /ParcelsPublic Notifications You must be signed in to change notification settings Fork 188 Star 360 Code Issues109 Pull requests13 Discussions Actions Projects Wiki Security and quality0 Insights Additional navigation options Code Issues Pull requests Discussions Actions Projects Wiki Security and quality Insights More xgrid generated datasets#2135 New issue Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community. Sign up for GitHub By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails. Already on GitHub? Sign in to your account Merged erikvansebille merged 16 commits into v4-dev from more_xgrid_generated_datasetsAug 6, 2025 Conversation 13 Commits 16 Checks 7 Files changed Uh oh! There was an error while loading. Please reload this page. Merged More xgrid generated datasets #2135 Show file tree Hide file tree Changes from all commits Commits Show all changes 16 commits Select commit Hold shift + click to select a range b8177fd Adding radial rotation flow erikvansebille Jul 29, 20253dfd8ac Adding radial rotation unit test erikvansebille Jul 31, 20256e192f8 Adding extra challenge of later particle starttime to test_radialrota… erikvansebille Aug 1, 2025c8b8bfd Refactoring moving_eddy test to generated dataset erikvansebille Aug 1, 202592d2a45 Making kernel loop work again for RK45 erikvansebille Aug 1, 2025c3bb197 Adding Stommel gyre example and test erikvansebille Aug 4, 20254dd3340 Adding test for multple variables setting on Particle erikvansebille Aug 4, 2025a0c02b4 Adding unit test for time-varying (decaying) moving eddy erikvansebille Aug 6, 2025f734195 Adding peninsula dataset to generated and advection tests erikvansebille Aug 6, 202539d1d63 Using np.testing.assert_allclose throughout advection tests erikvansebille Aug 6, 2025cc720ce Update parcels/xgrid.py erikvansebille Aug 6, 2025c58d7f3 Fixing review suggestion erikvansebille Aug 6, 20250287d7c Throwing error when add_constant is not called with float or int erikvansebille Aug 6, 2025568f133 Removing unit test that is now part of #2134 erikvansebille Aug 6, 20252218e8f Update parcels/fieldset.py erikvansebille Aug 6, 20254810e29 Merge branch 'v4-dev' into more_xgrid_generated_datasets erikvansebille Aug 6, 2025 File filter Filter by extension Filter by extension .py (8) All 1 file type selected Viewed files Clear filters Conversations Failed to load comments. Loading Uh oh! There was an error while loading. Please reload this page. Jump to Jump to file Failed to load files. Loading Uh oh! There was an error while loading. Please reload this page. Diff view Unified Split Hide whitespace Diff view Unified Split Hide whitespace parcels/_datasets/structured/generated.py generated.py parcels/application_kernels/advection.py advection.py parcels/application_kernels/interpolation.py interpolation.py parcels/fieldset.py fieldset.py parcels/kernel.py kernel.py parcels/xgrid.py xgrid.py tests/v4/test_advection.py test_advection.py tests/v4/test_diffusion.py test_diffusion.py There are no files selected for viewing 266 changes: 266 additions & 0 deletions 266 parcels/_datasets/structured/generated.py Show comments View file Open in desktop Original file line number Diff line number Diff line change @@ -1,3 +1,5 @@ import math import numpy as np import xarray as xr Expand All @@ -18,3 +20,267 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh_type="spherical"): "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes. lon = np.linspace(0, 60, xdim, dtype=np.float32) lat = np.linspace(0, 60, ydim, dtype=np.float32) x0 = 30.0 # Define the origin to be the centre of the Field. y0 = 30.0 U = np.zeros((2, 1, ydim, xdim), dtype=np.float32) V = np.zeros((2, 1, ydim, xdim), dtype=np.float32) omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day. for i in range(lon.size): for j in range(lat.size): r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2) assert r >= 0.0 assert r <= np.sqrt(x0**2 + y0**2) theta = np.arctan2((lat[j] - y0), (lon[i] - x0)) assert abs(theta) <= np.pi U[:, :, j, i] = r * np.sin(theta) * omega V[:, :, j, i] = -r * np.cos(theta) * omega return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1 Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow.""" f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants lon = np.linspace(0, 25000, xdim, dtype=np.float32) lat = np.linspace(0, 25000, ydim, dtype=np.float32) time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m")) U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) for t in range(len(time)): U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s"))) V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s"))) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, }, ) def decaying_moving_eddy_dataset(xdim=2, ydim=2): """Simulate an ocean that accelerates subject to Coriolis force and dissipative effects, upon which a geostrophic current is superimposed. The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ u_g = 0.04 # Geostrophic current u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0 gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity. gamma_g = 1.0 / (28.9 * 86400) f = 1.0e-4 # Coriolis parameter. time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m")) lon = np.linspace(0, 20000, xdim, dtype=np.float32) lat = np.linspace(5000, 12000, ydim, dtype=np.float32) U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) for t in range(time.size): t_float = time[t] / np.timedelta64(1, "s") U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float) V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, "gamma": gamma, "gamma_g": gamma_g, }, ) def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): """Construct a fieldset encapsulating the flow field around an idealised peninsula. Parameters ---------- xdim : Horizontal dimension of the generated fieldset ydim : Vertical dimension of the generated fieldset mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. 2. flat (default): No conversion, lat/lon are assumed to be in m. grid_type : Option whether grid is either Arakawa A (default) or C The original test description can be found in Fig. 2.2.3 in: North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of recommended practices for modelling physical - biological interactions during fish early life. ICES Cooperative Research Report No. 295. 111 pp. http://archimer.ifremer.fr/doc/00157/26792/24888.pdf """ domainsizeX, domainsizeY = (1.0e5, 5.0e4) La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32) Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32) u0 = 1 x0 = domainsizeX / 2 R = 0.32 * domainsizeX / 2 # Create the fields P = np.zeros((ydim, xdim), dtype=np.float32) U = np.zeros_like(P) V = np.zeros_like(P) x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy") P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y # Set land points to zero landpoints = P >= 0.0 P[landpoints] = 0.0 if grid_type == "A": U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2) V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2) U[landpoints] = 0.0 V[landpoints] = 0.0 Udims = ["YC", "XG"] Vdims = ["YG", "XC"] elif grid_type == "C": U = np.zeros(P.shape) V = np.zeros(P.shape) V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0]) U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0]) Udims = ["YG", "XG"] Vdims = ["YG", "XG"] else: raise RuntimeError(f"Grid_type {grid_type} is not a valid option") # Convert from m to lat/lon for spherical meshes lon = La / 1852.0 / 60.0 if mesh == "spherical" else La lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa return xr.Dataset( { "U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): """Simulate a periodic current along a western boundary, with significantly larger velocities along the western edge than the rest of the region The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ a = b = 10000 * 1e3 scalefac = 0.05 # to scale for physically meaningful velocities dx, dy = a / xdim, b / ydim # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(0, a, xdim, dtype=np.float32) lat = np.linspace(0, b, ydim, dtype=np.float32) # Define arrays U (zonal), V (meridional) and P (sea surface height) U = np.zeros((lat.size, lon.size), dtype=np.float32) V = np.zeros((lat.size, lon.size), dtype=np.float32) P = np.zeros((lat.size, lon.size), dtype=np.float32) beta = 2e-11 r = 1 / (11.6 * 86400) es = r / (beta * a) for j in range(lat.size): for i in range(lon.size): xi = lon[i] / a yi = lat[j] / b P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "A": U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "C": V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b Udims = ["YC", "XG"] Vdims = ["YG", "XC"] else: Udims = ["YG", "XG"] Vdims = ["YG", "XG"] return xr.Dataset( {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)}, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) 14 changes: 10 additions & 4 deletions 14 parcels/application_kernels/advection.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -115,9 +115,15 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ dt = min( particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt ) # TODO: improve API for converting dt to seconds dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds if dt > fieldset.RK45_max_dt: dt = fieldset.RK45_max_dt particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s") if dt < fieldset.RK45_min_dt: particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s") return StatusCode.Repeat particle.dt = particle.next_dt Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0], Expand DownExpand Up @@ -162,7 +168,7 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)): particle.dlon += lon_4th particle.dlat += lat_4th if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): particle.next_dt *= 2 else: particle.next_dt /= 2 Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/application_kernels/interpolation.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -40,7 +40,10 @@ def XBiLinear( zi, _ = position["Z"] data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2] data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] if tau > 0: data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] else: data = data[ti, :, :] return ( (1 - xsi) * (1 - eta) * data[0, 0] Expand Down 3 changes: 2 additions & 1 deletion 3 parcels/fieldset.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -164,7 +164,8 @@ def add_constant(self, name, value): """ if name in self.constants: raise ValueError(f"FieldSet already has a constant with name '{name}'") if not isinstance(value, (float, np.floating, int, np.integer)): raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}") self.constants[name] = np.float32(value) @property Expand Down 2 changes: 1 addition & 1 deletion 2 parcels/kernel.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime): res_tmp = f(p, self._fieldset, p.time_nextloop) if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode res = res_tmp if res == StatusCode.StopExecution: if res in [StatusCode.StopExecution, StatusCode.Repeat]: break if res is None: Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/xgrid.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -271,7 +271,10 @@ def _gtype(self): def search(self, z, y, x, ei=None): ds = self.xgcm_grid._ds zi, zeta = _search_1d_array(ds.depth.values, z) if "Z" in self.axes: zi, zeta = _search_1d_array(ds.depth.values, z) else: zi, zeta = 0, 0.0 if zi == -1: if zeta < 0: raise FieldOutOfBoundError( Expand Down Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later. Footer © 2026 GitHub, Inc. Footer navigation Terms Privacy Security Status Community Docs Contact You can’t perform that action at this time. , 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Highlight search terms from Google/DuckDuckGo/Bing referrer\n(function() {\n var ref = document.referrer;\n var terms = [];\n \n if (ref.includes('google.com') || ref.includes('duckduckgo.com') || ref.includes('bing.com')) {\n var url = new URL(ref);\n var q = url.searchParams.get('q') || url.searchParams.get('p');\n if (q) {\n terms = q.split(/\\s+/).filter(function(t) { return t.length > 2; });\n }\n }\n \n if (terms.length === 0) return;\n \n var style = document.createElement('style');\n style.textContent = '.userscript-highlight { background: #fbbf24; color: #1a1a2e; padding: 1px 3px; border-radius: 2px; }';\n document.head.appendChild(style);\n \n function highlight(node) {\n if (node.nodeType === 3) { // text node\n var text = node.textContent;\n var found = false;\n terms.forEach(function(term) {\n var regex = new RegExp('(' + term.replace(/[.*+?^${}()|[\\]\\\\]/g, '\\\\') + ')', 'gi');\n if (regex.test(text)) {\n found = true;\n var frag = document.createDocumentFragment();\n var parts = text.split(regex);\n parts.forEach(function(part, i) {\n if (i % 2 === 0) {\n frag.appendChild(document.createTextNode(part));\n } else {\n var span = document.createElement('span');\n span.className = 'userscript-highlight';\n span.textContent = part;\n frag.appendChild(span);\n }\n });\n node.parentNode.replaceChild(frag, node);\n }\n });\n } else if (node.nodeType === 1 && node.childNodes) { // element\n var skipTags = ['SCRIPT', 'STYLE', 'NOSCRIPT', 'TEXTAREA', 'INPUT', 'SELECT'];\n if (!skipTags.includes(node.tagName)) {\n Array.from(node.childNodes).forEach(highlight);\n }\n }\n }\n \n highlight(document.body);\n \n // Re-highlight on dynamic content\n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1 || node.nodeType === 3) highlight(node);\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Highlight Search Terms"); } } catch(__e) { console.warn('[Userscript:Highlight Search Terms]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + ' Skip to contentNavigation Menu Sign inAppearance settings AI CODE CREATIONGitHub CopilotWrite better code with AI GitHub Copilot appDirect agents from issue to merge MCP RegistryIntegrate external tools DEVELOPER WORKFLOWSActionsAutomate any workflow CodespacesInstant dev environments IssuesPlan and track work Code ReviewManage code changes Code QualityEnforce quality at merge APPLICATION SECURITYGitHub Advanced SecurityFind and fix vulnerabilities Code securitySecure your code as you build Secret protectionStop leaks before they start EXPLOREWhy GitHub Documentation Blog Changelog Marketplace View all features BY COMPANY SIZEEnterprises Small and medium teams Startups Nonprofits BY USE CASEApp Modernization DevSecOps DevOps CI/CD View all use cases BY INDUSTRYHealthcare Financial services Manufacturing Government View all industries View all solutions EXPLORE BY TOPICAI Software Development DevOps Security View all topics EXPLORE BY TYPECustomer stories Events & webinars Ebooks & reports Business insights GitHub Skills SUPPORT & SERVICESDocumentation Customer support Community forum Trust center Partners View all resources COMMUNITYGitHub SponsorsFund open source developers PROGRAMSSecurity Lab Maintainer Community GitHub Stars Archive Program REPOSITORIESTopics Trending Collections ENTERPRISE SOLUTIONSEnterprise platformAI-powered developer platform AVAILABLE ADD-ONSGitHub Advanced SecurityEnterprise-grade security features Copilot for BusinessEnterprise-grade AI features Premium SupportEnterprise-grade 24/7 support Pricing Sign in Sign upAppearance settings You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert Uh oh! There was an error while loading. Please reload this page. Parcels-code /ParcelsPublic Notifications You must be signed in to change notification settings Fork 188 Star 360 Code Issues109 Pull requests13 Discussions Actions Projects Wiki Security and quality0 Insights Additional navigation options Code Issues Pull requests Discussions Actions Projects Wiki Security and quality Insights More xgrid generated datasets#2135 New issue Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community. Sign up for GitHub By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails. Already on GitHub? Sign in to your account Merged erikvansebille merged 16 commits into v4-dev from more_xgrid_generated_datasetsAug 6, 2025 Conversation 13 Commits 16 Checks 7 Files changed Uh oh! There was an error while loading. Please reload this page. Merged More xgrid generated datasets #2135 Show file tree Hide file tree Changes from all commits Commits Show all changes 16 commits Select commit Hold shift + click to select a range b8177fd Adding radial rotation flow erikvansebille Jul 29, 20253dfd8ac Adding radial rotation unit test erikvansebille Jul 31, 20256e192f8 Adding extra challenge of later particle starttime to test_radialrota… erikvansebille Aug 1, 2025c8b8bfd Refactoring moving_eddy test to generated dataset erikvansebille Aug 1, 202592d2a45 Making kernel loop work again for RK45 erikvansebille Aug 1, 2025c3bb197 Adding Stommel gyre example and test erikvansebille Aug 4, 20254dd3340 Adding test for multple variables setting on Particle erikvansebille Aug 4, 2025a0c02b4 Adding unit test for time-varying (decaying) moving eddy erikvansebille Aug 6, 2025f734195 Adding peninsula dataset to generated and advection tests erikvansebille Aug 6, 202539d1d63 Using np.testing.assert_allclose throughout advection tests erikvansebille Aug 6, 2025cc720ce Update parcels/xgrid.py erikvansebille Aug 6, 2025c58d7f3 Fixing review suggestion erikvansebille Aug 6, 20250287d7c Throwing error when add_constant is not called with float or int erikvansebille Aug 6, 2025568f133 Removing unit test that is now part of #2134 erikvansebille Aug 6, 20252218e8f Update parcels/fieldset.py erikvansebille Aug 6, 20254810e29 Merge branch 'v4-dev' into more_xgrid_generated_datasets erikvansebille Aug 6, 2025 File filter Filter by extension Filter by extension .py (8) All 1 file type selected Viewed files Clear filters Conversations Failed to load comments. Loading Uh oh! There was an error while loading. Please reload this page. Jump to Jump to file Failed to load files. Loading Uh oh! There was an error while loading. Please reload this page. Diff view Unified Split Hide whitespace Diff view Unified Split Hide whitespace parcels/_datasets/structured/generated.py generated.py parcels/application_kernels/advection.py advection.py parcels/application_kernels/interpolation.py interpolation.py parcels/fieldset.py fieldset.py parcels/kernel.py kernel.py parcels/xgrid.py xgrid.py tests/v4/test_advection.py test_advection.py tests/v4/test_diffusion.py test_diffusion.py There are no files selected for viewing 266 changes: 266 additions & 0 deletions 266 parcels/_datasets/structured/generated.py Show comments View file Open in desktop Original file line number Diff line number Diff line change @@ -1,3 +1,5 @@ import math import numpy as np import xarray as xr Expand All @@ -18,3 +20,267 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh_type="spherical"): "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes. lon = np.linspace(0, 60, xdim, dtype=np.float32) lat = np.linspace(0, 60, ydim, dtype=np.float32) x0 = 30.0 # Define the origin to be the centre of the Field. y0 = 30.0 U = np.zeros((2, 1, ydim, xdim), dtype=np.float32) V = np.zeros((2, 1, ydim, xdim), dtype=np.float32) omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day. for i in range(lon.size): for j in range(lat.size): r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2) assert r >= 0.0 assert r <= np.sqrt(x0**2 + y0**2) theta = np.arctan2((lat[j] - y0), (lon[i] - x0)) assert abs(theta) <= np.pi U[:, :, j, i] = r * np.sin(theta) * omega V[:, :, j, i] = -r * np.cos(theta) * omega return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1 Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow.""" f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants lon = np.linspace(0, 25000, xdim, dtype=np.float32) lat = np.linspace(0, 25000, ydim, dtype=np.float32) time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m")) U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) for t in range(len(time)): U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s"))) V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s"))) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, }, ) def decaying_moving_eddy_dataset(xdim=2, ydim=2): """Simulate an ocean that accelerates subject to Coriolis force and dissipative effects, upon which a geostrophic current is superimposed. The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ u_g = 0.04 # Geostrophic current u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0 gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity. gamma_g = 1.0 / (28.9 * 86400) f = 1.0e-4 # Coriolis parameter. time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m")) lon = np.linspace(0, 20000, xdim, dtype=np.float32) lat = np.linspace(5000, 12000, ydim, dtype=np.float32) U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) for t in range(time.size): t_float = time[t] / np.timedelta64(1, "s") U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float) V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, "gamma": gamma, "gamma_g": gamma_g, }, ) def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): """Construct a fieldset encapsulating the flow field around an idealised peninsula. Parameters ---------- xdim : Horizontal dimension of the generated fieldset ydim : Vertical dimension of the generated fieldset mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. 2. flat (default): No conversion, lat/lon are assumed to be in m. grid_type : Option whether grid is either Arakawa A (default) or C The original test description can be found in Fig. 2.2.3 in: North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of recommended practices for modelling physical - biological interactions during fish early life. ICES Cooperative Research Report No. 295. 111 pp. http://archimer.ifremer.fr/doc/00157/26792/24888.pdf """ domainsizeX, domainsizeY = (1.0e5, 5.0e4) La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32) Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32) u0 = 1 x0 = domainsizeX / 2 R = 0.32 * domainsizeX / 2 # Create the fields P = np.zeros((ydim, xdim), dtype=np.float32) U = np.zeros_like(P) V = np.zeros_like(P) x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy") P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y # Set land points to zero landpoints = P >= 0.0 P[landpoints] = 0.0 if grid_type == "A": U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2) V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2) U[landpoints] = 0.0 V[landpoints] = 0.0 Udims = ["YC", "XG"] Vdims = ["YG", "XC"] elif grid_type == "C": U = np.zeros(P.shape) V = np.zeros(P.shape) V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0]) U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0]) Udims = ["YG", "XG"] Vdims = ["YG", "XG"] else: raise RuntimeError(f"Grid_type {grid_type} is not a valid option") # Convert from m to lat/lon for spherical meshes lon = La / 1852.0 / 60.0 if mesh == "spherical" else La lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa return xr.Dataset( { "U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): """Simulate a periodic current along a western boundary, with significantly larger velocities along the western edge than the rest of the region The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ a = b = 10000 * 1e3 scalefac = 0.05 # to scale for physically meaningful velocities dx, dy = a / xdim, b / ydim # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(0, a, xdim, dtype=np.float32) lat = np.linspace(0, b, ydim, dtype=np.float32) # Define arrays U (zonal), V (meridional) and P (sea surface height) U = np.zeros((lat.size, lon.size), dtype=np.float32) V = np.zeros((lat.size, lon.size), dtype=np.float32) P = np.zeros((lat.size, lon.size), dtype=np.float32) beta = 2e-11 r = 1 / (11.6 * 86400) es = r / (beta * a) for j in range(lat.size): for i in range(lon.size): xi = lon[i] / a yi = lat[j] / b P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "A": U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "C": V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b Udims = ["YC", "XG"] Vdims = ["YG", "XC"] else: Udims = ["YG", "XG"] Vdims = ["YG", "XG"] return xr.Dataset( {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)}, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) 14 changes: 10 additions & 4 deletions 14 parcels/application_kernels/advection.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -115,9 +115,15 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ dt = min( particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt ) # TODO: improve API for converting dt to seconds dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds if dt > fieldset.RK45_max_dt: dt = fieldset.RK45_max_dt particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s") if dt < fieldset.RK45_min_dt: particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s") return StatusCode.Repeat particle.dt = particle.next_dt Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0], Expand DownExpand Up @@ -162,7 +168,7 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)): particle.dlon += lon_4th particle.dlat += lat_4th if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): particle.next_dt *= 2 else: particle.next_dt /= 2 Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/application_kernels/interpolation.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -40,7 +40,10 @@ def XBiLinear( zi, _ = position["Z"] data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2] data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] if tau > 0: data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] else: data = data[ti, :, :] return ( (1 - xsi) * (1 - eta) * data[0, 0] Expand Down 3 changes: 2 additions & 1 deletion 3 parcels/fieldset.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -164,7 +164,8 @@ def add_constant(self, name, value): """ if name in self.constants: raise ValueError(f"FieldSet already has a constant with name '{name}'") if not isinstance(value, (float, np.floating, int, np.integer)): raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}") self.constants[name] = np.float32(value) @property Expand Down 2 changes: 1 addition & 1 deletion 2 parcels/kernel.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime): res_tmp = f(p, self._fieldset, p.time_nextloop) if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode res = res_tmp if res == StatusCode.StopExecution: if res in [StatusCode.StopExecution, StatusCode.Repeat]: break if res is None: Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/xgrid.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -271,7 +271,10 @@ def _gtype(self): def search(self, z, y, x, ei=None): ds = self.xgcm_grid._ds zi, zeta = _search_1d_array(ds.depth.values, z) if "Z" in self.axes: zi, zeta = _search_1d_array(ds.depth.values, z) else: zi, zeta = 0, 0.0 if zi == -1: if zeta < 0: raise FieldOutOfBoundError( Expand Down Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later. Footer © 2026 GitHub, Inc. Footer navigation Terms Privacy Security Status Community Docs Contact You can’t perform that action at this time. , 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Strip utm_, fbclid, gclid, etc. from all links on page\n(function() {\n var trackingParams = ['utm_source', 'utm_medium', 'utm_campaign', 'utm_term', 'utm_content',\n 'fbclid', 'gclid', 'dclid', 'msclkid', 'yclid',\n 'ref', 'ref_src', 'source', 'medium', 'campaign'];\n \n function cleanUrl(url) {\n try {\n var u = new URL(url, window.location.origin);\n var changed = false;\n trackingParams.forEach(function(p) {\n if (u.searchParams.has(p)) {\n u.searchParams.delete(p);\n changed = true;\n }\n });\n return changed ? u.toString() : url;\n } catch (e) {\n return url;\n }\n }\n \n function cleanLinks() {\n document.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n \n cleanLinks();\n \n var observer = new MutationObserver(function(mutations) {\n mutations.forEach(function(m) {\n m.addedNodes.forEach(function(node) {\n if (node.nodeType === 1) {\n if (node.tagName === 'A') cleanLinks();\n node.querySelectorAll('a[href]').forEach(function(a) {\n var clean = cleanUrl(a.href);\n if (clean !== a.href) a.href = clean;\n });\n }\n });\n });\n });\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "Remove Tracking Parameters from Links"); } } catch(__e) { console.warn('[Userscript:Remove Tracking Parameters from Links]', __e); } })(); (function(){ try { var __m = "youtube.com"; var __re = new RegExp('^' + "youtube\\.com" + ' Skip to contentNavigation Menu Sign inAppearance settings AI CODE CREATIONGitHub CopilotWrite better code with AI GitHub Copilot appDirect agents from issue to merge MCP RegistryIntegrate external tools DEVELOPER WORKFLOWSActionsAutomate any workflow CodespacesInstant dev environments IssuesPlan and track work Code ReviewManage code changes Code QualityEnforce quality at merge APPLICATION SECURITYGitHub Advanced SecurityFind and fix vulnerabilities Code securitySecure your code as you build Secret protectionStop leaks before they start EXPLOREWhy GitHub Documentation Blog Changelog Marketplace View all features BY COMPANY SIZEEnterprises Small and medium teams Startups Nonprofits BY USE CASEApp Modernization DevSecOps DevOps CI/CD View all use cases BY INDUSTRYHealthcare Financial services Manufacturing Government View all industries View all solutions EXPLORE BY TOPICAI Software Development DevOps Security View all topics EXPLORE BY TYPECustomer stories Events & webinars Ebooks & reports Business insights GitHub Skills SUPPORT & SERVICESDocumentation Customer support Community forum Trust center Partners View all resources COMMUNITYGitHub SponsorsFund open source developers PROGRAMSSecurity Lab Maintainer Community GitHub Stars Archive Program REPOSITORIESTopics Trending Collections ENTERPRISE SOLUTIONSEnterprise platformAI-powered developer platform AVAILABLE ADD-ONSGitHub Advanced SecurityEnterprise-grade security features Copilot for BusinessEnterprise-grade AI features Premium SupportEnterprise-grade 24/7 support Pricing Sign in Sign upAppearance settings You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert Uh oh! There was an error while loading. Please reload this page. Parcels-code /ParcelsPublic Notifications You must be signed in to change notification settings Fork 188 Star 360 Code Issues109 Pull requests13 Discussions Actions Projects Wiki Security and quality0 Insights Additional navigation options Code Issues Pull requests Discussions Actions Projects Wiki Security and quality Insights More xgrid generated datasets#2135 New issue Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community. Sign up for GitHub By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails. Already on GitHub? Sign in to your account Merged erikvansebille merged 16 commits into v4-dev from more_xgrid_generated_datasetsAug 6, 2025 Conversation 13 Commits 16 Checks 7 Files changed Uh oh! There was an error while loading. Please reload this page. Merged More xgrid generated datasets #2135 Show file tree Hide file tree Changes from all commits Commits Show all changes 16 commits Select commit Hold shift + click to select a range b8177fd Adding radial rotation flow erikvansebille Jul 29, 20253dfd8ac Adding radial rotation unit test erikvansebille Jul 31, 20256e192f8 Adding extra challenge of later particle starttime to test_radialrota… erikvansebille Aug 1, 2025c8b8bfd Refactoring moving_eddy test to generated dataset erikvansebille Aug 1, 202592d2a45 Making kernel loop work again for RK45 erikvansebille Aug 1, 2025c3bb197 Adding Stommel gyre example and test erikvansebille Aug 4, 20254dd3340 Adding test for multple variables setting on Particle erikvansebille Aug 4, 2025a0c02b4 Adding unit test for time-varying (decaying) moving eddy erikvansebille Aug 6, 2025f734195 Adding peninsula dataset to generated and advection tests erikvansebille Aug 6, 202539d1d63 Using np.testing.assert_allclose throughout advection tests erikvansebille Aug 6, 2025cc720ce Update parcels/xgrid.py erikvansebille Aug 6, 2025c58d7f3 Fixing review suggestion erikvansebille Aug 6, 20250287d7c Throwing error when add_constant is not called with float or int erikvansebille Aug 6, 2025568f133 Removing unit test that is now part of #2134 erikvansebille Aug 6, 20252218e8f Update parcels/fieldset.py erikvansebille Aug 6, 20254810e29 Merge branch 'v4-dev' into more_xgrid_generated_datasets erikvansebille Aug 6, 2025 File filter Filter by extension Filter by extension .py (8) All 1 file type selected Viewed files Clear filters Conversations Failed to load comments. Loading Uh oh! There was an error while loading. Please reload this page. Jump to Jump to file Failed to load files. Loading Uh oh! There was an error while loading. Please reload this page. Diff view Unified Split Hide whitespace Diff view Unified Split Hide whitespace parcels/_datasets/structured/generated.py generated.py parcels/application_kernels/advection.py advection.py parcels/application_kernels/interpolation.py interpolation.py parcels/fieldset.py fieldset.py parcels/kernel.py kernel.py parcels/xgrid.py xgrid.py tests/v4/test_advection.py test_advection.py tests/v4/test_diffusion.py test_diffusion.py There are no files selected for viewing 266 changes: 266 additions & 0 deletions 266 parcels/_datasets/structured/generated.py Show comments View file Open in desktop Original file line number Diff line number Diff line change @@ -1,3 +1,5 @@ import math import numpy as np import xarray as xr Expand All @@ -18,3 +20,267 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh_type="spherical"): "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes. lon = np.linspace(0, 60, xdim, dtype=np.float32) lat = np.linspace(0, 60, ydim, dtype=np.float32) x0 = 30.0 # Define the origin to be the centre of the Field. y0 = 30.0 U = np.zeros((2, 1, ydim, xdim), dtype=np.float32) V = np.zeros((2, 1, ydim, xdim), dtype=np.float32) omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day. for i in range(lon.size): for j in range(lat.size): r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2) assert r >= 0.0 assert r <= np.sqrt(x0**2 + y0**2) theta = np.arctan2((lat[j] - y0), (lon[i] - x0)) assert abs(theta) <= np.pi U[:, :, j, i] = r * np.sin(theta) * omega V[:, :, j, i] = -r * np.cos(theta) * omega return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1 Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow.""" f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants lon = np.linspace(0, 25000, xdim, dtype=np.float32) lat = np.linspace(0, 25000, ydim, dtype=np.float32) time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m")) U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) for t in range(len(time)): U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s"))) V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s"))) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, }, ) def decaying_moving_eddy_dataset(xdim=2, ydim=2): """Simulate an ocean that accelerates subject to Coriolis force and dissipative effects, upon which a geostrophic current is superimposed. The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ u_g = 0.04 # Geostrophic current u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0 gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity. gamma_g = 1.0 / (28.9 * 86400) f = 1.0e-4 # Coriolis parameter. time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m")) lon = np.linspace(0, 20000, xdim, dtype=np.float32) lat = np.linspace(5000, 12000, ydim, dtype=np.float32) U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) for t in range(time.size): t_float = time[t] / np.timedelta64(1, "s") U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float) V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, "gamma": gamma, "gamma_g": gamma_g, }, ) def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): """Construct a fieldset encapsulating the flow field around an idealised peninsula. Parameters ---------- xdim : Horizontal dimension of the generated fieldset ydim : Vertical dimension of the generated fieldset mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. 2. flat (default): No conversion, lat/lon are assumed to be in m. grid_type : Option whether grid is either Arakawa A (default) or C The original test description can be found in Fig. 2.2.3 in: North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of recommended practices for modelling physical - biological interactions during fish early life. ICES Cooperative Research Report No. 295. 111 pp. http://archimer.ifremer.fr/doc/00157/26792/24888.pdf """ domainsizeX, domainsizeY = (1.0e5, 5.0e4) La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32) Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32) u0 = 1 x0 = domainsizeX / 2 R = 0.32 * domainsizeX / 2 # Create the fields P = np.zeros((ydim, xdim), dtype=np.float32) U = np.zeros_like(P) V = np.zeros_like(P) x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy") P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y # Set land points to zero landpoints = P >= 0.0 P[landpoints] = 0.0 if grid_type == "A": U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2) V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2) U[landpoints] = 0.0 V[landpoints] = 0.0 Udims = ["YC", "XG"] Vdims = ["YG", "XC"] elif grid_type == "C": U = np.zeros(P.shape) V = np.zeros(P.shape) V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0]) U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0]) Udims = ["YG", "XG"] Vdims = ["YG", "XG"] else: raise RuntimeError(f"Grid_type {grid_type} is not a valid option") # Convert from m to lat/lon for spherical meshes lon = La / 1852.0 / 60.0 if mesh == "spherical" else La lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa return xr.Dataset( { "U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): """Simulate a periodic current along a western boundary, with significantly larger velocities along the western edge than the rest of the region The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ a = b = 10000 * 1e3 scalefac = 0.05 # to scale for physically meaningful velocities dx, dy = a / xdim, b / ydim # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(0, a, xdim, dtype=np.float32) lat = np.linspace(0, b, ydim, dtype=np.float32) # Define arrays U (zonal), V (meridional) and P (sea surface height) U = np.zeros((lat.size, lon.size), dtype=np.float32) V = np.zeros((lat.size, lon.size), dtype=np.float32) P = np.zeros((lat.size, lon.size), dtype=np.float32) beta = 2e-11 r = 1 / (11.6 * 86400) es = r / (beta * a) for j in range(lat.size): for i in range(lon.size): xi = lon[i] / a yi = lat[j] / b P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "A": U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "C": V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b Udims = ["YC", "XG"] Vdims = ["YG", "XC"] else: Udims = ["YG", "XG"] Vdims = ["YG", "XG"] return xr.Dataset( {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)}, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) 14 changes: 10 additions & 4 deletions 14 parcels/application_kernels/advection.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -115,9 +115,15 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ dt = min( particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt ) # TODO: improve API for converting dt to seconds dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds if dt > fieldset.RK45_max_dt: dt = fieldset.RK45_max_dt particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s") if dt < fieldset.RK45_min_dt: particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s") return StatusCode.Repeat particle.dt = particle.next_dt Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0], Expand DownExpand Up @@ -162,7 +168,7 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)): particle.dlon += lon_4th particle.dlat += lat_4th if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): particle.next_dt *= 2 else: particle.next_dt /= 2 Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/application_kernels/interpolation.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -40,7 +40,10 @@ def XBiLinear( zi, _ = position["Z"] data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2] data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] if tau > 0: data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] else: data = data[ti, :, :] return ( (1 - xsi) * (1 - eta) * data[0, 0] Expand Down 3 changes: 2 additions & 1 deletion 3 parcels/fieldset.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -164,7 +164,8 @@ def add_constant(self, name, value): """ if name in self.constants: raise ValueError(f"FieldSet already has a constant with name '{name}'") if not isinstance(value, (float, np.floating, int, np.integer)): raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}") self.constants[name] = np.float32(value) @property Expand Down 2 changes: 1 addition & 1 deletion 2 parcels/kernel.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime): res_tmp = f(p, self._fieldset, p.time_nextloop) if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode res = res_tmp if res == StatusCode.StopExecution: if res in [StatusCode.StopExecution, StatusCode.Repeat]: break if res is None: Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/xgrid.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -271,7 +271,10 @@ def _gtype(self): def search(self, z, y, x, ei=None): ds = self.xgcm_grid._ds zi, zeta = _search_1d_array(ds.depth.values, z) if "Z" in self.axes: zi, zeta = _search_1d_array(ds.depth.values, z) else: zi, zeta = 0, 0.0 if zi == -1: if zeta < 0: raise FieldOutOfBoundError( Expand Down Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later. Footer © 2026 GitHub, Inc. Footer navigation Terms Privacy Security Status Community Docs Contact You can’t perform that action at this time. , 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Auto-enable theater mode on YouTube\n(function() {\n function tryTheater() {\n var btn = document.querySelector('button[aria-label=\"Theater mode\"], ytd-player #player button[title=\"Theater mode\"]');\n if (btn && !btn.classList.contains('activated')) {\n btn.click();\n }\n }\n \n // Try immediately\n tryTheater();\n \n // Try after navigation (SPA)\n var lastUrl = location.href;\n setInterval(function() {\n if (location.href !== lastUrl) {\n lastUrl = location.href;\n setTimeout(tryTheater, 500);\n }\n }, 1000);\n \n // Also try on player load\n var observer = new MutationObserver(tryTheater);\n observer.observe(document.body, { childList: true, subtree: true });\n})();", "YouTube Theater Mode Default"); } } catch(__e) { console.warn('[Userscript:YouTube Theater Mode Default]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + ' Skip to contentNavigation Menu Sign inAppearance settings AI CODE CREATIONGitHub CopilotWrite better code with AI GitHub Copilot appDirect agents from issue to merge MCP RegistryIntegrate external tools DEVELOPER WORKFLOWSActionsAutomate any workflow CodespacesInstant dev environments IssuesPlan and track work Code ReviewManage code changes Code QualityEnforce quality at merge APPLICATION SECURITYGitHub Advanced SecurityFind and fix vulnerabilities Code securitySecure your code as you build Secret protectionStop leaks before they start EXPLOREWhy GitHub Documentation Blog Changelog Marketplace View all features BY COMPANY SIZEEnterprises Small and medium teams Startups Nonprofits BY USE CASEApp Modernization DevSecOps DevOps CI/CD View all use cases BY INDUSTRYHealthcare Financial services Manufacturing Government View all industries View all solutions EXPLORE BY TOPICAI Software Development DevOps Security View all topics EXPLORE BY TYPECustomer stories Events & webinars Ebooks & reports Business insights GitHub Skills SUPPORT & SERVICESDocumentation Customer support Community forum Trust center Partners View all resources COMMUNITYGitHub SponsorsFund open source developers PROGRAMSSecurity Lab Maintainer Community GitHub Stars Archive Program REPOSITORIESTopics Trending Collections ENTERPRISE SOLUTIONSEnterprise platformAI-powered developer platform AVAILABLE ADD-ONSGitHub Advanced SecurityEnterprise-grade security features Copilot for BusinessEnterprise-grade AI features Premium SupportEnterprise-grade 24/7 support Pricing Sign in Sign upAppearance settings You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert Uh oh! There was an error while loading. Please reload this page. Parcels-code /ParcelsPublic Notifications You must be signed in to change notification settings Fork 188 Star 360 Code Issues109 Pull requests13 Discussions Actions Projects Wiki Security and quality0 Insights Additional navigation options Code Issues Pull requests Discussions Actions Projects Wiki Security and quality Insights More xgrid generated datasets#2135 New issue Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community. Sign up for GitHub By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails. Already on GitHub? Sign in to your account Merged erikvansebille merged 16 commits into v4-dev from more_xgrid_generated_datasetsAug 6, 2025 Conversation 13 Commits 16 Checks 7 Files changed Uh oh! There was an error while loading. Please reload this page. Merged More xgrid generated datasets #2135 Show file tree Hide file tree Changes from all commits Commits Show all changes 16 commits Select commit Hold shift + click to select a range b8177fd Adding radial rotation flow erikvansebille Jul 29, 20253dfd8ac Adding radial rotation unit test erikvansebille Jul 31, 20256e192f8 Adding extra challenge of later particle starttime to test_radialrota… erikvansebille Aug 1, 2025c8b8bfd Refactoring moving_eddy test to generated dataset erikvansebille Aug 1, 202592d2a45 Making kernel loop work again for RK45 erikvansebille Aug 1, 2025c3bb197 Adding Stommel gyre example and test erikvansebille Aug 4, 20254dd3340 Adding test for multple variables setting on Particle erikvansebille Aug 4, 2025a0c02b4 Adding unit test for time-varying (decaying) moving eddy erikvansebille Aug 6, 2025f734195 Adding peninsula dataset to generated and advection tests erikvansebille Aug 6, 202539d1d63 Using np.testing.assert_allclose throughout advection tests erikvansebille Aug 6, 2025cc720ce Update parcels/xgrid.py erikvansebille Aug 6, 2025c58d7f3 Fixing review suggestion erikvansebille Aug 6, 20250287d7c Throwing error when add_constant is not called with float or int erikvansebille Aug 6, 2025568f133 Removing unit test that is now part of #2134 erikvansebille Aug 6, 20252218e8f Update parcels/fieldset.py erikvansebille Aug 6, 20254810e29 Merge branch 'v4-dev' into more_xgrid_generated_datasets erikvansebille Aug 6, 2025 File filter Filter by extension Filter by extension .py (8) All 1 file type selected Viewed files Clear filters Conversations Failed to load comments. Loading Uh oh! There was an error while loading. Please reload this page. Jump to Jump to file Failed to load files. Loading Uh oh! There was an error while loading. Please reload this page. Diff view Unified Split Hide whitespace Diff view Unified Split Hide whitespace parcels/_datasets/structured/generated.py generated.py parcels/application_kernels/advection.py advection.py parcels/application_kernels/interpolation.py interpolation.py parcels/fieldset.py fieldset.py parcels/kernel.py kernel.py parcels/xgrid.py xgrid.py tests/v4/test_advection.py test_advection.py tests/v4/test_diffusion.py test_diffusion.py There are no files selected for viewing 266 changes: 266 additions & 0 deletions 266 parcels/_datasets/structured/generated.py Show comments View file Open in desktop Original file line number Diff line number Diff line change @@ -1,3 +1,5 @@ import math import numpy as np import xarray as xr Expand All @@ -18,3 +20,267 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh_type="spherical"): "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes. lon = np.linspace(0, 60, xdim, dtype=np.float32) lat = np.linspace(0, 60, ydim, dtype=np.float32) x0 = 30.0 # Define the origin to be the centre of the Field. y0 = 30.0 U = np.zeros((2, 1, ydim, xdim), dtype=np.float32) V = np.zeros((2, 1, ydim, xdim), dtype=np.float32) omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day. for i in range(lon.size): for j in range(lat.size): r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2) assert r >= 0.0 assert r <= np.sqrt(x0**2 + y0**2) theta = np.arctan2((lat[j] - y0), (lon[i] - x0)) assert abs(theta) <= np.pi U[:, :, j, i] = r * np.sin(theta) * omega V[:, :, j, i] = -r * np.cos(theta) * omega return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1 Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow.""" f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants lon = np.linspace(0, 25000, xdim, dtype=np.float32) lat = np.linspace(0, 25000, ydim, dtype=np.float32) time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m")) U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) for t in range(len(time)): U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s"))) V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s"))) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, }, ) def decaying_moving_eddy_dataset(xdim=2, ydim=2): """Simulate an ocean that accelerates subject to Coriolis force and dissipative effects, upon which a geostrophic current is superimposed. The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ u_g = 0.04 # Geostrophic current u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0 gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity. gamma_g = 1.0 / (28.9 * 86400) f = 1.0e-4 # Coriolis parameter. time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m")) lon = np.linspace(0, 20000, xdim, dtype=np.float32) lat = np.linspace(5000, 12000, ydim, dtype=np.float32) U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) for t in range(time.size): t_float = time[t] / np.timedelta64(1, "s") U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float) V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, "gamma": gamma, "gamma_g": gamma_g, }, ) def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): """Construct a fieldset encapsulating the flow field around an idealised peninsula. Parameters ---------- xdim : Horizontal dimension of the generated fieldset ydim : Vertical dimension of the generated fieldset mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. 2. flat (default): No conversion, lat/lon are assumed to be in m. grid_type : Option whether grid is either Arakawa A (default) or C The original test description can be found in Fig. 2.2.3 in: North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of recommended practices for modelling physical - biological interactions during fish early life. ICES Cooperative Research Report No. 295. 111 pp. http://archimer.ifremer.fr/doc/00157/26792/24888.pdf """ domainsizeX, domainsizeY = (1.0e5, 5.0e4) La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32) Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32) u0 = 1 x0 = domainsizeX / 2 R = 0.32 * domainsizeX / 2 # Create the fields P = np.zeros((ydim, xdim), dtype=np.float32) U = np.zeros_like(P) V = np.zeros_like(P) x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy") P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y # Set land points to zero landpoints = P >= 0.0 P[landpoints] = 0.0 if grid_type == "A": U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2) V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2) U[landpoints] = 0.0 V[landpoints] = 0.0 Udims = ["YC", "XG"] Vdims = ["YG", "XC"] elif grid_type == "C": U = np.zeros(P.shape) V = np.zeros(P.shape) V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0]) U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0]) Udims = ["YG", "XG"] Vdims = ["YG", "XG"] else: raise RuntimeError(f"Grid_type {grid_type} is not a valid option") # Convert from m to lat/lon for spherical meshes lon = La / 1852.0 / 60.0 if mesh == "spherical" else La lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa return xr.Dataset( { "U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): """Simulate a periodic current along a western boundary, with significantly larger velocities along the western edge than the rest of the region The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ a = b = 10000 * 1e3 scalefac = 0.05 # to scale for physically meaningful velocities dx, dy = a / xdim, b / ydim # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(0, a, xdim, dtype=np.float32) lat = np.linspace(0, b, ydim, dtype=np.float32) # Define arrays U (zonal), V (meridional) and P (sea surface height) U = np.zeros((lat.size, lon.size), dtype=np.float32) V = np.zeros((lat.size, lon.size), dtype=np.float32) P = np.zeros((lat.size, lon.size), dtype=np.float32) beta = 2e-11 r = 1 / (11.6 * 86400) es = r / (beta * a) for j in range(lat.size): for i in range(lon.size): xi = lon[i] / a yi = lat[j] / b P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "A": U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "C": V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b Udims = ["YC", "XG"] Vdims = ["YG", "XC"] else: Udims = ["YG", "XG"] Vdims = ["YG", "XG"] return xr.Dataset( {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)}, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) 14 changes: 10 additions & 4 deletions 14 parcels/application_kernels/advection.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -115,9 +115,15 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ dt = min( particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt ) # TODO: improve API for converting dt to seconds dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds if dt > fieldset.RK45_max_dt: dt = fieldset.RK45_max_dt particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s") if dt < fieldset.RK45_min_dt: particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s") return StatusCode.Repeat particle.dt = particle.next_dt Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0], Expand DownExpand Up @@ -162,7 +168,7 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)): particle.dlon += lon_4th particle.dlat += lat_4th if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): particle.next_dt *= 2 else: particle.next_dt /= 2 Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/application_kernels/interpolation.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -40,7 +40,10 @@ def XBiLinear( zi, _ = position["Z"] data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2] data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] if tau > 0: data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] else: data = data[ti, :, :] return ( (1 - xsi) * (1 - eta) * data[0, 0] Expand Down 3 changes: 2 additions & 1 deletion 3 parcels/fieldset.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -164,7 +164,8 @@ def add_constant(self, name, value): """ if name in self.constants: raise ValueError(f"FieldSet already has a constant with name '{name}'") if not isinstance(value, (float, np.floating, int, np.integer)): raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}") self.constants[name] = np.float32(value) @property Expand Down 2 changes: 1 addition & 1 deletion 2 parcels/kernel.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime): res_tmp = f(p, self._fieldset, p.time_nextloop) if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode res = res_tmp if res == StatusCode.StopExecution: if res in [StatusCode.StopExecution, StatusCode.Repeat]: break if res is None: Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/xgrid.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -271,7 +271,10 @@ def _gtype(self): def search(self, z, y, x, ei=None): ds = self.xgcm_grid._ds zi, zeta = _search_1d_array(ds.depth.values, z) if "Z" in self.axes: zi, zeta = _search_1d_array(ds.depth.values, z) else: zi, zeta = 0, 0.0 if zi == -1: if zeta < 0: raise FieldOutOfBoundError( Expand Down Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later. Footer © 2026 GitHub, Inc. Footer navigation Terms Privacy Security Status Community Docs Contact You can’t perform that action at this time. , 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Remove or un-stick sticky/fixed headers that block content\n(function() {\n function unstick() {\n document.querySelectorAll('header, nav, [role=\"banner\"], .header, .navbar, .sticky, .fixed-top, [style*=\"position: fixed\"], [style*=\"position:sticky\"]').forEach(function(el) {\n if (el.style.position === 'fixed' || el.style.position === 'sticky' || \n getComputedStyle(el).position === 'fixed' || getComputedStyle(el).position === 'sticky') {\n el.style.position = 'static';\n el.style.top = 'auto';\n el.style.zIndex = 'auto';\n }\n });\n }\n \n unstick();\n \n var observer = new MutationObserver(unstick);\n observer.observe(document.body, { childList: true, subtree: true, attributes: true, attributeFilter: ['style', 'class'] });\n})();", "Kill Sticky Headers"); } } catch(__e) { console.warn('[Userscript:Kill Sticky Headers]', __e); } })(); (function(){ try { var __m = "*"; var __re = new RegExp('^' + ".*" + ' Skip to contentNavigation Menu Sign inAppearance settings AI CODE CREATIONGitHub CopilotWrite better code with AI GitHub Copilot appDirect agents from issue to merge MCP RegistryIntegrate external tools DEVELOPER WORKFLOWSActionsAutomate any workflow CodespacesInstant dev environments IssuesPlan and track work Code ReviewManage code changes Code QualityEnforce quality at merge APPLICATION SECURITYGitHub Advanced SecurityFind and fix vulnerabilities Code securitySecure your code as you build Secret protectionStop leaks before they start EXPLOREWhy GitHub Documentation Blog Changelog Marketplace View all features BY COMPANY SIZEEnterprises Small and medium teams Startups Nonprofits BY USE CASEApp Modernization DevSecOps DevOps CI/CD View all use cases BY INDUSTRYHealthcare Financial services Manufacturing Government View all industries View all solutions EXPLORE BY TOPICAI Software Development DevOps Security View all topics EXPLORE BY TYPECustomer stories Events & webinars Ebooks & reports Business insights GitHub Skills SUPPORT & SERVICESDocumentation Customer support Community forum Trust center Partners View all resources COMMUNITYGitHub SponsorsFund open source developers PROGRAMSSecurity Lab Maintainer Community GitHub Stars Archive Program REPOSITORIESTopics Trending Collections ENTERPRISE SOLUTIONSEnterprise platformAI-powered developer platform AVAILABLE ADD-ONSGitHub Advanced SecurityEnterprise-grade security features Copilot for BusinessEnterprise-grade AI features Premium SupportEnterprise-grade 24/7 support Pricing Sign in Sign upAppearance settings You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert Uh oh! There was an error while loading. Please reload this page. Parcels-code /ParcelsPublic Notifications You must be signed in to change notification settings Fork 188 Star 360 Code Issues109 Pull requests13 Discussions Actions Projects Wiki Security and quality0 Insights Additional navigation options Code Issues Pull requests Discussions Actions Projects Wiki Security and quality Insights More xgrid generated datasets#2135 New issue Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community. Sign up for GitHub By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails. Already on GitHub? Sign in to your account Merged erikvansebille merged 16 commits into v4-dev from more_xgrid_generated_datasetsAug 6, 2025 Conversation 13 Commits 16 Checks 7 Files changed Uh oh! There was an error while loading. Please reload this page. Merged More xgrid generated datasets #2135 Show file tree Hide file tree Changes from all commits Commits Show all changes 16 commits Select commit Hold shift + click to select a range b8177fd Adding radial rotation flow erikvansebille Jul 29, 20253dfd8ac Adding radial rotation unit test erikvansebille Jul 31, 20256e192f8 Adding extra challenge of later particle starttime to test_radialrota… erikvansebille Aug 1, 2025c8b8bfd Refactoring moving_eddy test to generated dataset erikvansebille Aug 1, 202592d2a45 Making kernel loop work again for RK45 erikvansebille Aug 1, 2025c3bb197 Adding Stommel gyre example and test erikvansebille Aug 4, 20254dd3340 Adding test for multple variables setting on Particle erikvansebille Aug 4, 2025a0c02b4 Adding unit test for time-varying (decaying) moving eddy erikvansebille Aug 6, 2025f734195 Adding peninsula dataset to generated and advection tests erikvansebille Aug 6, 202539d1d63 Using np.testing.assert_allclose throughout advection tests erikvansebille Aug 6, 2025cc720ce Update parcels/xgrid.py erikvansebille Aug 6, 2025c58d7f3 Fixing review suggestion erikvansebille Aug 6, 20250287d7c Throwing error when add_constant is not called with float or int erikvansebille Aug 6, 2025568f133 Removing unit test that is now part of #2134 erikvansebille Aug 6, 20252218e8f Update parcels/fieldset.py erikvansebille Aug 6, 20254810e29 Merge branch 'v4-dev' into more_xgrid_generated_datasets erikvansebille Aug 6, 2025 File filter Filter by extension Filter by extension .py (8) All 1 file type selected Viewed files Clear filters Conversations Failed to load comments. Loading Uh oh! There was an error while loading. Please reload this page. Jump to Jump to file Failed to load files. Loading Uh oh! There was an error while loading. Please reload this page. Diff view Unified Split Hide whitespace Diff view Unified Split Hide whitespace parcels/_datasets/structured/generated.py generated.py parcels/application_kernels/advection.py advection.py parcels/application_kernels/interpolation.py interpolation.py parcels/fieldset.py fieldset.py parcels/kernel.py kernel.py parcels/xgrid.py xgrid.py tests/v4/test_advection.py test_advection.py tests/v4/test_diffusion.py test_diffusion.py There are no files selected for viewing 266 changes: 266 additions & 0 deletions 266 parcels/_datasets/structured/generated.py Show comments View file Open in desktop Original file line number Diff line number Diff line change @@ -1,3 +1,5 @@ import math import numpy as np import xarray as xr Expand All @@ -18,3 +20,267 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh_type="spherical"): "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes. lon = np.linspace(0, 60, xdim, dtype=np.float32) lat = np.linspace(0, 60, ydim, dtype=np.float32) x0 = 30.0 # Define the origin to be the centre of the Field. y0 = 30.0 U = np.zeros((2, 1, ydim, xdim), dtype=np.float32) V = np.zeros((2, 1, ydim, xdim), dtype=np.float32) omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day. for i in range(lon.size): for j in range(lat.size): r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2) assert r >= 0.0 assert r <= np.sqrt(x0**2 + y0**2) theta = np.arctan2((lat[j] - y0), (lon[i] - x0)) assert abs(theta) <= np.pi U[:, :, j, i] = r * np.sin(theta) * omega V[:, :, j, i] = -r * np.cos(theta) * omega return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1 Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow.""" f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants lon = np.linspace(0, 25000, xdim, dtype=np.float32) lat = np.linspace(0, 25000, ydim, dtype=np.float32) time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m")) U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) for t in range(len(time)): U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s"))) V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s"))) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, }, ) def decaying_moving_eddy_dataset(xdim=2, ydim=2): """Simulate an ocean that accelerates subject to Coriolis force and dissipative effects, upon which a geostrophic current is superimposed. The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ u_g = 0.04 # Geostrophic current u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0 gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity. gamma_g = 1.0 / (28.9 * 86400) f = 1.0e-4 # Coriolis parameter. time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m")) lon = np.linspace(0, 20000, xdim, dtype=np.float32) lat = np.linspace(5000, 12000, ydim, dtype=np.float32) U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) for t in range(time.size): t_float = time[t] / np.timedelta64(1, "s") U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float) V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, "gamma": gamma, "gamma_g": gamma_g, }, ) def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): """Construct a fieldset encapsulating the flow field around an idealised peninsula. Parameters ---------- xdim : Horizontal dimension of the generated fieldset ydim : Vertical dimension of the generated fieldset mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. 2. flat (default): No conversion, lat/lon are assumed to be in m. grid_type : Option whether grid is either Arakawa A (default) or C The original test description can be found in Fig. 2.2.3 in: North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of recommended practices for modelling physical - biological interactions during fish early life. ICES Cooperative Research Report No. 295. 111 pp. http://archimer.ifremer.fr/doc/00157/26792/24888.pdf """ domainsizeX, domainsizeY = (1.0e5, 5.0e4) La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32) Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32) u0 = 1 x0 = domainsizeX / 2 R = 0.32 * domainsizeX / 2 # Create the fields P = np.zeros((ydim, xdim), dtype=np.float32) U = np.zeros_like(P) V = np.zeros_like(P) x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy") P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y # Set land points to zero landpoints = P >= 0.0 P[landpoints] = 0.0 if grid_type == "A": U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2) V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2) U[landpoints] = 0.0 V[landpoints] = 0.0 Udims = ["YC", "XG"] Vdims = ["YG", "XC"] elif grid_type == "C": U = np.zeros(P.shape) V = np.zeros(P.shape) V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0]) U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0]) Udims = ["YG", "XG"] Vdims = ["YG", "XG"] else: raise RuntimeError(f"Grid_type {grid_type} is not a valid option") # Convert from m to lat/lon for spherical meshes lon = La / 1852.0 / 60.0 if mesh == "spherical" else La lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa return xr.Dataset( { "U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): """Simulate a periodic current along a western boundary, with significantly larger velocities along the western edge than the rest of the region The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ a = b = 10000 * 1e3 scalefac = 0.05 # to scale for physically meaningful velocities dx, dy = a / xdim, b / ydim # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(0, a, xdim, dtype=np.float32) lat = np.linspace(0, b, ydim, dtype=np.float32) # Define arrays U (zonal), V (meridional) and P (sea surface height) U = np.zeros((lat.size, lon.size), dtype=np.float32) V = np.zeros((lat.size, lon.size), dtype=np.float32) P = np.zeros((lat.size, lon.size), dtype=np.float32) beta = 2e-11 r = 1 / (11.6 * 86400) es = r / (beta * a) for j in range(lat.size): for i in range(lon.size): xi = lon[i] / a yi = lat[j] / b P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "A": U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "C": V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b Udims = ["YC", "XG"] Vdims = ["YG", "XC"] else: Udims = ["YG", "XG"] Vdims = ["YG", "XG"] return xr.Dataset( {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)}, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) 14 changes: 10 additions & 4 deletions 14 parcels/application_kernels/advection.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -115,9 +115,15 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ dt = min( particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt ) # TODO: improve API for converting dt to seconds dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds if dt > fieldset.RK45_max_dt: dt = fieldset.RK45_max_dt particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s") if dt < fieldset.RK45_min_dt: particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s") return StatusCode.Repeat particle.dt = particle.next_dt Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0], Expand DownExpand Up @@ -162,7 +168,7 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)): particle.dlon += lon_4th particle.dlat += lat_4th if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): particle.next_dt *= 2 else: particle.next_dt /= 2 Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/application_kernels/interpolation.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -40,7 +40,10 @@ def XBiLinear( zi, _ = position["Z"] data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2] data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] if tau > 0: data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] else: data = data[ti, :, :] return ( (1 - xsi) * (1 - eta) * data[0, 0] Expand Down 3 changes: 2 additions & 1 deletion 3 parcels/fieldset.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -164,7 +164,8 @@ def add_constant(self, name, value): """ if name in self.constants: raise ValueError(f"FieldSet already has a constant with name '{name}'") if not isinstance(value, (float, np.floating, int, np.integer)): raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}") self.constants[name] = np.float32(value) @property Expand Down 2 changes: 1 addition & 1 deletion 2 parcels/kernel.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime): res_tmp = f(p, self._fieldset, p.time_nextloop) if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode res = res_tmp if res == StatusCode.StopExecution: if res in [StatusCode.StopExecution, StatusCode.Repeat]: break if res is None: Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/xgrid.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -271,7 +271,10 @@ def _gtype(self): def search(self, z, y, x, ei=None): ds = self.xgcm_grid._ds zi, zeta = _search_1d_array(ds.depth.values, z) if "Z" in self.axes: zi, zeta = _search_1d_array(ds.depth.values, z) else: zi, zeta = 0, 0.0 if zi == -1: if zeta < 0: raise FieldOutOfBoundError( Expand Down Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later. Footer © 2026 GitHub, Inc. Footer navigation Terms Privacy Security Status Community Docs Contact You can’t perform that action at this time. , 'i'); if (__m === '*' || __re.test(location.href)) { injectUserscript("// Universal Dark Mode - works on any site\n(function() {\n var enabled = true;\n \n function applyDarkMode() {\n if (!enabled) return;\n \n // Create style element if it doesn't exist\n var style = document.getElementById('universal-dark-mode-style');\n if (!style) {\n style = document.createElement('style');\n style.id = 'universal-dark-mode-style';\n document.head.appendChild(style);\n }\n \n // Dark mode CSS - inverts colors but preserves images/video\n style.textContent = '\n /* Invert everything except media */\n html {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #1a1a2e !important;\n }\n \n /* Restore images, videos, iframes, canvas */\n img, video, iframe, canvas, svg, picture, [style*=\"background-image\"] {\n filter: invert(1) hue-rotate(180deg) !important;\n }\n \n /* Preserve specific elements that should not be inverted */\n .no-dark-mode, .no-dark-mode *,\n [data-theme=\"light\"], [data-theme=\"light\"],\n .ace_editor, .ace_editor *,\n .CodeMirror, .CodeMirror *,\n .monaco-editor, .monaco-editor *,\n .markdown-body pre, .markdown-body pre *,\n .highlight, .highlight *,\n pre code, pre code * {\n filter: none !important;\n }\n \n /* Fix common UI elements */\n .modal, .popup, .dropdown-menu, .tooltip, .popover {\n filter: invert(1) hue-rotate(180deg) !important;\n background: #2d2d44 !important;\n border-color: #444 !important;\n }\n \n /* Scrollbars */\n ::-webkit-scrollbar { background: #1a1a2e !important; }\n ::-webkit-scrollbar-thumb { background: #444 !important; }\n ::-webkit-scrollbar-thumb:hover { background: #555 !important; }\n \n /* Selection */\n ::selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ::-moz-selection { background: #4ecdc4 !important; color: #1a1a2e !important; }\n ';\n }\n \n function removeDarkMode() {\n var style = document.getElementById('universal-dark-mode-style');\n if (style) style.remove();\n }\n \n // Toggle with Alt+Shift+D\n document.addEventListener('keydown', function(e) {\n if (e.altKey && e.shiftKey && e.key === 'D') {\n e.preventDefault();\n enabled = !enabled;\n if (enabled) {\n applyDarkMode();\n console.log('[Universal Dark Mode] Enabled');\n } else {\n removeDarkMode();\n console.log('[Universal Dark Mode] Disabled');\n }\n }\n });\n \n // Apply on load\n applyDarkMode();\n \n // Re-apply on dynamic content\n var observer = new MutationObserver(function(mutations) {\n if (enabled && !document.getElementById('universal-dark-mode-style')) {\n applyDarkMode();\n }\n });\n observer.observe(document.head, { childList: true });\n \n console.log('[Universal Dark Mode] Loaded - Press Alt+Shift+D to toggle');\n})();", "Universal Dark Mode"); } } catch(__e) { console.warn('[Userscript:Universal Dark Mode]', __e); } })(); })(); Skip to contentNavigation Menu Sign inAppearance settings AI CODE CREATIONGitHub CopilotWrite better code with AI GitHub Copilot appDirect agents from issue to merge MCP RegistryIntegrate external tools DEVELOPER WORKFLOWSActionsAutomate any workflow CodespacesInstant dev environments IssuesPlan and track work Code ReviewManage code changes Code QualityEnforce quality at merge APPLICATION SECURITYGitHub Advanced SecurityFind and fix vulnerabilities Code securitySecure your code as you build Secret protectionStop leaks before they start EXPLOREWhy GitHub Documentation Blog Changelog Marketplace View all features BY COMPANY SIZEEnterprises Small and medium teams Startups Nonprofits BY USE CASEApp Modernization DevSecOps DevOps CI/CD View all use cases BY INDUSTRYHealthcare Financial services Manufacturing Government View all industries View all solutions EXPLORE BY TOPICAI Software Development DevOps Security View all topics EXPLORE BY TYPECustomer stories Events & webinars Ebooks & reports Business insights GitHub Skills SUPPORT & SERVICESDocumentation Customer support Community forum Trust center Partners View all resources COMMUNITYGitHub SponsorsFund open source developers PROGRAMSSecurity Lab Maintainer Community GitHub Stars Archive Program REPOSITORIESTopics Trending Collections ENTERPRISE SOLUTIONSEnterprise platformAI-powered developer platform AVAILABLE ADD-ONSGitHub Advanced SecurityEnterprise-grade security features Copilot for BusinessEnterprise-grade AI features Premium SupportEnterprise-grade 24/7 support Pricing Sign in Sign upAppearance settings You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert Uh oh! There was an error while loading. Please reload this page. Parcels-code /ParcelsPublic Notifications You must be signed in to change notification settings Fork 188 Star 360 Code Issues109 Pull requests13 Discussions Actions Projects Wiki Security and quality0 Insights Additional navigation options Code Issues Pull requests Discussions Actions Projects Wiki Security and quality Insights More xgrid generated datasets#2135 New issue Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community. Sign up for GitHub By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails. Already on GitHub? Sign in to your account Merged erikvansebille merged 16 commits into v4-dev from more_xgrid_generated_datasetsAug 6, 2025 Conversation 13 Commits 16 Checks 7 Files changed Uh oh! There was an error while loading. Please reload this page. Merged More xgrid generated datasets #2135 Show file tree Hide file tree Changes from all commits Commits Show all changes 16 commits Select commit Hold shift + click to select a range b8177fd Adding radial rotation flow erikvansebille Jul 29, 20253dfd8ac Adding radial rotation unit test erikvansebille Jul 31, 20256e192f8 Adding extra challenge of later particle starttime to test_radialrota… erikvansebille Aug 1, 2025c8b8bfd Refactoring moving_eddy test to generated dataset erikvansebille Aug 1, 202592d2a45 Making kernel loop work again for RK45 erikvansebille Aug 1, 2025c3bb197 Adding Stommel gyre example and test erikvansebille Aug 4, 20254dd3340 Adding test for multple variables setting on Particle erikvansebille Aug 4, 2025a0c02b4 Adding unit test for time-varying (decaying) moving eddy erikvansebille Aug 6, 2025f734195 Adding peninsula dataset to generated and advection tests erikvansebille Aug 6, 202539d1d63 Using np.testing.assert_allclose throughout advection tests erikvansebille Aug 6, 2025cc720ce Update parcels/xgrid.py erikvansebille Aug 6, 2025c58d7f3 Fixing review suggestion erikvansebille Aug 6, 20250287d7c Throwing error when add_constant is not called with float or int erikvansebille Aug 6, 2025568f133 Removing unit test that is now part of #2134 erikvansebille Aug 6, 20252218e8f Update parcels/fieldset.py erikvansebille Aug 6, 20254810e29 Merge branch 'v4-dev' into more_xgrid_generated_datasets erikvansebille Aug 6, 2025 File filter Filter by extension Filter by extension .py (8) All 1 file type selected Viewed files Clear filters Conversations Failed to load comments. Loading Uh oh! There was an error while loading. Please reload this page. Jump to Jump to file Failed to load files. Loading Uh oh! There was an error while loading. Please reload this page. Diff view Unified Split Hide whitespace Diff view Unified Split Hide whitespace parcels/_datasets/structured/generated.py generated.py parcels/application_kernels/advection.py advection.py parcels/application_kernels/interpolation.py interpolation.py parcels/fieldset.py fieldset.py parcels/kernel.py kernel.py parcels/xgrid.py xgrid.py tests/v4/test_advection.py test_advection.py tests/v4/test_diffusion.py test_diffusion.py There are no files selected for viewing 266 changes: 266 additions & 0 deletions 266 parcels/_datasets/structured/generated.py Show comments View file Open in desktop Original file line number Diff line number Diff line change @@ -1,3 +1,5 @@ import math import numpy as np import xarray as xr Expand All @@ -18,3 +20,267 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh_type="spherical"): "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def radial_rotation_dataset(xdim=200, ydim=200): # Define 2D flat, square fieldset for testing purposes. lon = np.linspace(0, 60, xdim, dtype=np.float32) lat = np.linspace(0, 60, ydim, dtype=np.float32) x0 = 30.0 # Define the origin to be the centre of the Field. y0 = 30.0 U = np.zeros((2, 1, ydim, xdim), dtype=np.float32) V = np.zeros((2, 1, ydim, xdim), dtype=np.float32) omega = 2 * np.pi / 86400.0 # Define the rotational period as 1 day. for i in range(lon.size): for j in range(lat.size): r = np.sqrt((lon[i] - x0) ** 2 + (lat[j] - y0) ** 2) assert r >= 0.0 assert r <= np.sqrt(x0**2 + y0**2) theta = np.arctan2((lat[j] - y0), (lon[i] - x0)) assert abs(theta) <= np.pi U[:, :, j, i] = r * np.sin(theta) * omega V[:, :, j, i] = -r * np.cos(theta) * omega return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], [np.timedelta64(0, "s"), np.timedelta64(10, "D")], {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def moving_eddy_dataset(xdim=2, ydim=2): # TODO check if this also works with xdim=1, ydim=1 Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. """Create a dataset with an eddy moving in time. Note that there is no spatial variation in the flow.""" f, u_0, u_g = 1.0e-4, 0.3, 0.04 # Some constants lon = np.linspace(0, 25000, xdim, dtype=np.float32) lat = np.linspace(0, 25000, ydim, dtype=np.float32) time = np.arange(np.timedelta64(0, "s"), np.timedelta64(7, "h"), np.timedelta64(1, "m")) U = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) V = np.zeros((len(time), 1, ydim, xdim), dtype=np.float32) for t in range(len(time)): U[t, :, :, :] = u_g + (u_0 - u_g) * np.cos(f * (time[t] / np.timedelta64(1, "s"))) V[t, :, :, :] = -(u_0 - u_g) * np.sin(f * (time[t] / np.timedelta64(1, "s"))) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, }, ) def decaying_moving_eddy_dataset(xdim=2, ydim=2): """Simulate an ocean that accelerates subject to Coriolis force and dissipative effects, upon which a geostrophic current is superimposed. The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ u_g = 0.04 # Geostrophic current u_0 = 0.3 # Initial speed in x dirrection. v_0 = 0 gamma = 1.0 / (2.89 * 86400) # Dissipitave effects due to viscousity. gamma_g = 1.0 / (28.9 * 86400) f = 1.0e-4 # Coriolis parameter. time = np.arange(np.timedelta64(0, "s"), np.timedelta64(1, "D") + np.timedelta64(1, "h"), np.timedelta64(2, "m")) lon = np.linspace(0, 20000, xdim, dtype=np.float32) lat = np.linspace(5000, 12000, ydim, dtype=np.float32) U = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) V = np.zeros((time.size, 1, lat.size, lon.size), dtype=np.float32) for t in range(time.size): t_float = time[t] / np.timedelta64(1, "s") U[t, :, :, :] = u_g * np.exp(-gamma_g * t_float) + (u_0 - u_g) * np.exp(-gamma * t_float) * np.cos(f * t_float) V[t, :, :, :] = -(u_0 - u_g) * np.exp(-gamma * t_float) * np.sin(f * t_float) return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], U), "V": (["time", "depth", "YG", "XG"], V)}, coords={ "time": (["time"], time, {"axis": "T"}), "depth": (["depth"], np.array([0.0]), {"axis": "Z"}), "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, attrs={ "u_0": u_0, "u_g": u_g, "f": f, "gamma": gamma, "gamma_g": gamma_g, }, ) def peninsula_dataset(xdim=100, ydim=50, mesh="flat", grid_type="A"): """Construct a fieldset encapsulating the flow field around an idealised peninsula. Parameters ---------- xdim : Horizontal dimension of the generated fieldset ydim : Vertical dimension of the generated fieldset mesh : str String indicating the type of mesh coordinates and units used during velocity interpolation: 1. spherical: Lat and lon in degree, with a correction for zonal velocity U near the poles. 2. flat (default): No conversion, lat/lon are assumed to be in m. grid_type : Option whether grid is either Arakawa A (default) or C The original test description can be found in Fig. 2.2.3 in: North, E. W., Gallego, A., Petitgas, P. (Eds). 2009. Manual of recommended practices for modelling physical - biological interactions during fish early life. ICES Cooperative Research Report No. 295. 111 pp. http://archimer.ifremer.fr/doc/00157/26792/24888.pdf """ domainsizeX, domainsizeY = (1.0e5, 5.0e4) La = np.linspace(1e3, domainsizeX, xdim, dtype=np.float32) Wa = np.linspace(1e3, domainsizeY, ydim, dtype=np.float32) u0 = 1 x0 = domainsizeX / 2 R = 0.32 * domainsizeX / 2 # Create the fields P = np.zeros((ydim, xdim), dtype=np.float32) U = np.zeros_like(P) V = np.zeros_like(P) x, y = np.meshgrid(La, Wa, sparse=True, indexing="xy") P[:, :] = u0 * R**2 * y / ((x - x0) ** 2 + y**2) - u0 * y # Set land points to zero landpoints = P >= 0.0 P[landpoints] = 0.0 if grid_type == "A": U[:, :] = u0 - u0 * R**2 * ((x - x0) ** 2 - y**2) / (((x - x0) ** 2 + y**2) ** 2) V[:, :] = -2 * u0 * R**2 * ((x - x0) * y) / (((x - x0) ** 2 + y**2) ** 2) U[landpoints] = 0.0 V[landpoints] = 0.0 Udims = ["YC", "XG"] Vdims = ["YG", "XC"] elif grid_type == "C": U = np.zeros(P.shape) V = np.zeros(P.shape) V[:, 1:] = (P[:, 1:] - P[:, :-1]) / (La[1] - La[0]) U[1:, :] = -(P[1:, :] - P[:-1, :]) / (Wa[1] - Wa[0]) Udims = ["YG", "XG"] Vdims = ["YG", "XG"] else: raise RuntimeError(f"Grid_type {grid_type} is not a valid option") # Convert from m to lat/lon for spherical meshes lon = La / 1852.0 / 60.0 if mesh == "spherical" else La lat = Wa / 1852.0 / 60.0 if mesh == "spherical" else Wa return xr.Dataset( { "U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) def stommel_gyre_dataset(xdim=200, ydim=200, grid_type="A"): """Simulate a periodic current along a western boundary, with significantly larger velocities along the western edge than the rest of the region The original test description can be found in: N. Fabbroni, 2009, Numerical Simulation of Passive tracers dispersion in the sea, Ph.D. dissertation, University of Bologna http://amsdottorato.unibo.it/1733/1/Fabbroni_Nicoletta_Tesi.pdf """ a = b = 10000 * 1e3 scalefac = 0.05 # to scale for physically meaningful velocities dx, dy = a / xdim, b / ydim # Coordinates of the test fieldset (on A-grid in deg) lon = np.linspace(0, a, xdim, dtype=np.float32) lat = np.linspace(0, b, ydim, dtype=np.float32) # Define arrays U (zonal), V (meridional) and P (sea surface height) U = np.zeros((lat.size, lon.size), dtype=np.float32) V = np.zeros((lat.size, lon.size), dtype=np.float32) P = np.zeros((lat.size, lon.size), dtype=np.float32) beta = 2e-11 r = 1 / (11.6 * 86400) es = r / (beta * a) for j in range(lat.size): for i in range(lon.size): xi = lon[i] / a yi = lat[j] / b P[j, i] = (1 - math.exp(-xi / es) - xi) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "A": U[j, i] = -(1 - math.exp(-xi / es) - xi) * math.pi**2 * np.cos(math.pi * yi) * scalefac V[j, i] = (math.exp(-xi / es) / es - 1) * math.pi * np.sin(math.pi * yi) * scalefac if grid_type == "C": V[:, 1:] = (P[:, 1:] - P[:, 0:-1]) / dx * a U[1:, :] = -(P[1:, :] - P[0:-1, :]) / dy * b Udims = ["YC", "XG"] Vdims = ["YG", "XC"] else: Udims = ["YG", "XG"] Vdims = ["YG", "XG"] return xr.Dataset( {"U": (Udims, U), "V": (Vdims, V), "P": (["YG", "XG"], P)}, coords={ "YC": (["YC"], np.arange(ydim) + 0.5, {"axis": "Y"}), "YG": (["YG"], np.arange(ydim), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(xdim) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(xdim), {"axis": "X", "c_grid_axis_shift": -0.5}), "lat": (["YG"], lat, {"axis": "Y", "c_grid_axis_shift": 0.5}), "lon": (["XG"], lon, {"axis": "X", "c_grid_axis_shift": -0.5}), }, ) 14 changes: 10 additions & 4 deletions 14 parcels/application_kernels/advection.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -115,9 +115,15 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover Time-step dt is halved if error is larger than fieldset.RK45_tol, and doubled if error is smaller than 1/10th of tolerance. """ dt = min( particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt ) # TODO: improve API for converting dt to seconds dt = particle.next_dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds if dt > fieldset.RK45_max_dt: dt = fieldset.RK45_max_dt particle.next_dt = fieldset.RK45_max_dt * np.timedelta64(1, "s") if dt < fieldset.RK45_min_dt: particle.next_dt = fieldset.RK45_min_dt * np.timedelta64(1, "s") return StatusCode.Repeat particle.dt = particle.next_dt Comment threaderikvansebille marked this conversation as resolved. Uh oh! There was an error while loading. Please reload this page. c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0] A = [ [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0], Expand DownExpand Up @@ -162,7 +168,7 @@ def AdvectionRK45(particle, fieldset, time): # pragma: no cover if (kappa <= fieldset.RK45_tol) or (math.fabs(dt) < math.fabs(fieldset.RK45_min_dt)): particle.dlon += lon_4th particle.dlat += lat_4th if (kappa <= fieldset.RK45_tol) / 10 and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): if (kappa <= fieldset.RK45_tol / 10) and (math.fabs(dt * 2) <= math.fabs(fieldset.RK45_max_dt)): particle.next_dt *= 2 else: particle.next_dt /= 2 Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/application_kernels/interpolation.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -40,7 +40,10 @@ def XBiLinear( zi, _ = position["Z"] data = field.data.data[:, zi, yi : yi + 2, xi : xi + 2] data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] if tau > 0: data = (1 - tau) * data[ti, :, :] + tau * data[ti + 1, :, :] else: data = data[ti, :, :] return ( (1 - xsi) * (1 - eta) * data[0, 0] Expand Down 3 changes: 2 additions & 1 deletion 3 parcels/fieldset.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -164,7 +164,8 @@ def add_constant(self, name, value): """ if name in self.constants: raise ValueError(f"FieldSet already has a constant with name '{name}'") if not isinstance(value, (float, np.floating, int, np.integer)): raise ValueError(f"FieldSet constants have to be of type float or int, got a {type(value)}") self.constants[name] = np.float32(value) @property Expand Down 2 changes: 1 addition & 1 deletion 2 parcels/kernel.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -303,7 +303,7 @@ def evaluate_particle(self, p, endtime): res_tmp = f(p, self._fieldset, p.time_nextloop) if res_tmp is not None: # TODO v4: Remove once all kernels return StatusCode res = res_tmp if res == StatusCode.StopExecution: if res in [StatusCode.StopExecution, StatusCode.Repeat]: break if res is None: Expand Down 5 changes: 4 additions & 1 deletion 5 parcels/xgrid.py Show comments View file Open in desktop Original file line number Diff line number Diff line change Expand Up @@ -271,7 +271,10 @@ def _gtype(self): def search(self, z, y, x, ei=None): ds = self.xgcm_grid._ds zi, zeta = _search_1d_array(ds.depth.values, z) if "Z" in self.axes: zi, zeta = _search_1d_array(ds.depth.values, z) else: zi, zeta = 0, 0.0 if zi == -1: if zeta < 0: raise FieldOutOfBoundError( Expand Down Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Loading Oops, something went wrong. Uh oh! There was an error while loading. Please reload this page. Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later. Footer © 2026 GitHub, Inc. Footer navigation Terms Privacy Security Status Community Docs Contact You can’t perform that action at this time.