[DO NOT MERGE] Feature/fieldset as uxarray obj #1934

fluidnumerics-joe · 2025-03-11T17:01:30Z

Chose the correct base branch (main for v3 changes, v4-dev for v4 changes)
Fixes #
Added tests
Added documentation

This is needed to provide (at least) a `_check_complete` call when passing a uxarray dataset to a particleset during initialization

Add draft of an eval method for the UxFieldSet class. Basic interpolation schemes are provided for node registered and face registered fields.

This commit contains a number of "to do" items marked in comments, with more incoming soon. A draft PR will be opened, FOR COMPARISON PURPOSES ONLY, so that we can assess changes required for supporting uxarray

for more information, see https://pre-commit.ci

…cels/Parcels into feature/fieldset-as-uxarray-obj

fluidnumerics-joe · 2025-03-11T17:28:29Z

parcels/application_kernels/advection.py

@@ -14,6 +16,15 @@
 ]


+def UxAdvectionEuler(particle, fieldset: UXFieldSet, time):


We'll need to have advection routines for UXArray that are distinct from XArray, unless we are able to get uniformity in the velocity field naming conventions between structured and unstructured data sets.

fluidnumerics-joe · 2025-03-11T17:29:42Z

parcels/kernel.py

-                    g._load_chunk = np.where(
-                        g._load_chunk == g._chunk_loaded_touched, g._chunk_deprecated, g._load_chunk
-                    )
+        # if pset.fieldset is not None:


This didn't seem to be needing to execute for the UXArray data, so I've just commented this out. Not sure what the plans are here, given that I recall us discussing removing chunking and deferred loading.

This is gone in some of the recent updates

fluidnumerics-joe · 2025-03-11T17:31:33Z

parcels/particleset.py

@@ -148,7 +149,10 @@ def ArrayClass_init(self, *args, **kwargs):
            pid_orig = np.arange(lon.size)

        if depth is None:
-            mindepth = self.fieldset.gridset.dimrange("depth")[0]
+            if type(self.fieldset) == UXFieldSet:
+                mindepth = 0  # TO DO :  get the min depth from the fieldset.uxgrid


I suppose we'll need a new implementation for fieldset.gridset.dimrange . Since the new proposed fieldset is meant to wrap around either an XArray or UXArray dataset, which each have a grid baked into the structure, the gridset subclass here is irrelevant.

fluidnumerics-joe · 2025-03-11T17:34:54Z

parcels/particleset.py

+            lonlatdepth_dtype = np.float32  # To do : get precision from fieldset
+        else:
+            if isinstance(fieldset.U, Field) and (not fieldset.U.allow_time_extrapolation):
+                _warn_particle_times_outside_fieldset_time_bounds(time, fieldset.U.grid.time_full)


I suppose we'll need an updated implementation of _warn_particle_times_outside_fieldset_time_bounds .. Alternatively, we'll need to replace the fieldset.U.grid.time_full with a fieldset.time_full call.

Note that I'm assuming that a UXarray or XArray dataset have a single time dimension for all fields in the dataset. Getting the time extents for the dataset wouldn't depend on the individual fields (e.g. U or V); rather it'd be associated with the dimensions of the dataset.

fluidnumerics-joe · 2025-03-11T17:35:37Z

parcels/particleset.py

-        if lonlatdepth_dtype is None:
-            lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U)
+        if type(fieldset) == UXFieldSet:
+            lonlatdepth_dtype = np.float32  # To do : get precision from fieldset


An obvious to-do item here; we need to be able to get the precision from the XArray or UXArray dataset underneath the thin fieldset wrapper.

fluidnumerics-joe · 2025-03-11T17:36:34Z

parcels/particleset.py

@@ -191,7 +200,10 @@ def ArrayClass_init(self, *args, **kwargs):
            self._repeatkwargs = kwargs
            self._repeatkwargs.pop("partition_function", None)

-        ngrids = fieldset.gridset.size
+        if type(fieldset) == UXFieldSet:
+            ngrids = 1


Here, we explicitly set ngrids to 1. However, with the idea that the fieldset will store a list of either XArray or UXArray datasets, we will want to set ngrids to the length of the list of datasets.

fluidnumerics-joe · 2025-03-11T17:36:57Z

parcels/particleset.py

+        if type(fieldset) == UXFieldSet:
+            ngrids = 1
+        else:
+            ngrids = fieldset.gridset.size


gridset attribute will be going away with the updated fieldset.

fluidnumerics-joe · 2025-03-11T17:37:52Z

parcels/particleset.py

@@ -967,7 +979,10 @@ def execute(
        if runtime is not None and endtime is not None:
            raise RuntimeError("Only one of (endtime, runtime) can be specified")

-        mintime, maxtime = self.fieldset.gridset.dimrange("time_full")
+        if type(self.fieldset) == UXFieldSet:
+            mintime, maxtime = self.fieldset.get_time_range()


fieldset.get_time_range is meant to replace fieldset.gridset.dimrange("time_full")

fluidnumerics-joe · 2025-03-11T17:39:06Z

parcels/particleset.py

@@ -1037,18 +1052,19 @@ def execute(

            time_at_startofloop = time

-            next_input = self.fieldset.computeTimeChunk(time, dt)
+            #next_input = self.fieldset.computeTimeChunk(time, dt)


Some work is going to be needed to get this implemented. Though I suspect some work is being done in this area currently in the v4-dev. Would be best to discuss together.

fluidnumerics-joe · 2025-03-11T17:39:54Z

parcels/particleset.py


+            next_time = endtime


This was hard-coded, just to get things working for the simple case with one time level in the UXArray dataset.

fluidnumerics-joe · 2025-03-11T17:41:50Z

parcels/uxfieldset.py

+    """A FieldSet class that holds hydrodynamic data needed to execute particles
+    in a UXArray.Dataset"""
+    # Change uxds to ds_list - which is a list of either uxDataset or xarray dataset
+    def __init__(self, uxds: ux.UxDataset, time_origin: float | np.datetime64 | np.timedelta64 | cftime.datetime = 0):


The uxds argument will be changed to ds (standing for "DataSets"), which will be a List(uxarray.UxDataset | xarray.Dataset)

fluidnumerics-joe · 2025-03-11T17:43:21Z

parcels/uxfieldset.py

+        assert self.uxds.v is not None, "UXFieldSet does not provide v velocity data"
+        assert self.uxds.uxgrid is not None, "UXFieldSet does not provide a grid"
+
+    def _face_interp(self, field, time, z, y, x, ei):


For now, hard-coded interpolation schemes are provided for face registered data and node registered data.

We may want to discuss with the Geomar folks if this dichotomy is intuitive for them or if they have a different way of viewing the interpolation problem.

fluidnumerics-joe · 2025-03-11T17:45:28Z

parcels/uxfieldset.py

+            if not self._point_is_in_face(y,x,ei):
+                # If the point is not in the previously defined face, then
+                # search for the face again.
+                # To do : Update the search here to do nearest neighbors search, rather than spatial hash - [email protected]


For UXArray data, we ought to add a method that does a nearest neighbors search, rather than a hash query when the particle is not found in the ei face.

I suppose this code here, and the proposed code, ought to placed in a _search_uxindices method

fluidnumerics-joe · 2025-03-11T17:46:21Z

parcels/uxfieldset.py

+                r = self._face_interp(field, time, z, y, x, fi)
+            else:
+                r = self._node_interp(field, time, z, y, x, fi)
+


For 3-D interpolation, we need to add vertical interpolation. Currently, all interpolation uses the fields at the upper most vertical layer

fluidnumerics-joe · 2025-03-11T17:47:49Z

parcels/uxfieldset.py

+                # If the point is not in the previously defined face, then
+                # search for the face again.
+                # To do : Update the search here to do nearest neighbors search, rather than spatial hash - [email protected]
+                print(f"Position : {x}, {y}")


print statements can probable come out now.. these were needed for debugging..

fluidnumerics-joe · 2025-03-11T17:48:36Z

parcels/uxfieldset.py

+        vim1 = nodes[i - 1]
+        vi = nodes[i]
+        vi1 = nodes[(i + 1) % n]
+        a0 = _triangle_area(vim1, vi, vi1)


These areas need to be clipped with a minimum value. As we found in uxarray, when a particle lies on a face edge or vertex, these areas can be zero, leading to division by zero.

fluidnumerics-joe · 2025-03-12T16:15:12Z

parcels/application_kernels/advection.py

+def UxAdvectionEuler(particle, fieldset: UXFieldSet, time):
+    """Advection of particles using Explicit Euler (aka Euler Forward) integration.
+    on an unstructured grid."""
+    vel, ei = fieldset.eval(["u","v"],time,particle.depth,particle.lat,particle.lon, particle.ei[0])


Push particle.ei to eval and select which igrid.

Also, push the xarray/uxarray "wrapper" down to the field class

fluidnumerics-joe added 4 commits February 27, 2025 16:48

Add simple wrapper for uxarray dataset in a uxfieldset

0dfa7e1

This is needed to provide (at least) a `_check_complete` call when passing a uxarray dataset to a particleset during initialization

Fix particleset checks to work with uxfieldset

6821e20

Add uxadvectionEuler method for euler forward integration

f1802ce

Add draft of an eval method for the UxFieldSet class. Basic interpolation schemes are provided for node registered and face registered fields.

Adds fixes to force this to work for FESOM2 example

e1cbc09

This commit contains a number of "to do" items marked in comments, with more incoming soon. A draft PR will be opened, FOR COMPARISON PURPOSES ONLY, so that we can assess changes required for supporting uxarray

fluidnumerics-joe marked this pull request as draft March 11, 2025 17:01

pre-commit-ci bot and others added 3 commits March 11, 2025 17:02

[pre-commit.ci] auto fixes from pre-commit.com hooks

4caf74c

for more information, see https://pre-commit.ci

Merge branch 'feature/fieldset-as-uxarray-obj' of github.com:OceanPar…

43f93c1

…cels/Parcels into feature/fieldset-as-uxarray-obj

Fix missed merge conflicts

d683511

fluidnumerics-joe commented Mar 11, 2025

View reviewed changes

fluidnumerics-joe commented Mar 12, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[DO NOT MERGE] Feature/fieldset as uxarray obj #1934

[DO NOT MERGE] Feature/fieldset as uxarray obj #1934

fluidnumerics-joe commented Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 12, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025 •

edited

Loading

fluidnumerics-joe Mar 11, 2025 •

edited

Loading

fluidnumerics-joe Mar 11, 2025 •

edited

Loading

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 11, 2025

fluidnumerics-joe Mar 12, 2025

		@@ -14,6 +16,15 @@
		]


		def UxAdvectionEuler(particle, fieldset: UXFieldSet, time):

[DO NOT MERGE] Feature/fieldset as uxarray obj #1934

Are you sure you want to change the base?

[DO NOT MERGE] Feature/fieldset as uxarray obj #1934

Conversation

fluidnumerics-joe commented Mar 11, 2025

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fluidnumerics-joe Mar 11, 2025 • edited Loading

Choose a reason for hiding this comment

fluidnumerics-joe Mar 11, 2025 • edited Loading

Choose a reason for hiding this comment

fluidnumerics-joe Mar 11, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fluidnumerics-joe Mar 11, 2025 •

edited

Loading

fluidnumerics-joe Mar 11, 2025 •

edited

Loading

fluidnumerics-joe Mar 11, 2025 •

edited

Loading