Implement XConstantField interpolation

Author: erikvansebille

src/parcels/_core/fieldset.py

@@ -17,6 +17,7 @@
 from parcels._core.xgrid import _DEFAULT_XGCM_KWARGS, XGrid
 from parcels._logger import logger
 from parcels._typing import Mesh
+from parcels.interpolators import XConstantField
 
 if TYPE_CHECKING:
     from parcels._core.basegrid import BaseGrid
@@ -116,7 +117,7 @@ def add_field(self, field: Field, name: str | None = None):
 
         self.fields[name] = field
 
-    def add_constant_field(self, name: str, value, mesh: Mesh = "flat"):
+    def add_constant_field(self, name: str, value, mesh: Mesh = "spherical"):
         """Wrapper function to add a Field that is constant in space,
            useful e.g. when using constant horizontal diffusivity
 
@@ -134,14 +135,24 @@ def add_constant_field(self, name: str, value, mesh: Mesh = "flat"):
                correction for zonal velocity U near the poles.
             2. flat: No conversion, lat/lon are assumed to be in m.
         """
-        ds = xr.Dataset({name: (["time", "lat", "lon", "depth"], np.full((1, 1, 1, 1), value))})
-        grid = XGrid(xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS))
+        ds = xr.Dataset(
+            {name: (["time", "depth", "lat", "lon"], np.full((1, 1, 1, 1), value))},
+            coords={
+                "time": (["time"], [np.timedelta64(0, "s")], {"axis": "T"}),
+                "depth": (["depth"], [0], {"axis": "Z"}),
+                "lat": (["lat"], [0], {"axis": "Y", "c_grid_axis_shift": -0.5}),
+                "lon": (["lon"], [0], {"axis": "X", "c_grid_axis_shift": -0.5}),
+                "lat_C": (["lat_C"], [0.5], {"axis": "Y"}),  # TODO why is this needed?
+                "lon_C": (["lon_C"], [0.5], {"axis": "X"}),  # TODO why is this needed?
+            },
+        )
+        grid = XGrid(xgcm.Grid(ds, **_DEFAULT_XGCM_KWARGS), mesh=mesh)
         self.add_field(
             Field(
                 name,
                 ds[name],
                 grid,
-                interp_method=None,  # TODO : Need to define an interpolation method for constants
+                interp_method=XConstantField,
             )
         )
 

src/parcels/interpolators.py

@@ -20,6 +20,7 @@
     "CGrid_Velocity",
     "UXPiecewiseConstantFace",
     "UXPiecewiseLinearNode",
+    "XConstantField",
     "XFreeslip",
     "XLinear",
     "XLinearInvdistLandTracer",
@@ -135,6 +136,15 @@ def XLinear(
     return value.compute() if is_dask_collection(value) else value
 
 
+def XConstantField(
+    particle_positions: dict[str, float | np.ndarray],
+    grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]],
+    field: Field,
+):
+    """Returning the single value of a Constant Field (with a size=(1,1,1,1) array)"""
+    return field.data[0, 0, 0, 0].values
+
+
 def CGrid_Velocity(
     particle_positions: dict[str, float | np.ndarray],
     grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]],

tests/test_diffusion.py

@@ -23,12 +23,10 @@ def test_fieldKh_Brownian(mesh):
     grid = XGrid.from_dataset(ds, mesh=mesh)
     U = Field("U", ds["U"], grid, interp_method=XLinear)
     V = Field("V", ds["V"], grid, interp_method=XLinear)
-    ds["Kh_zonal"] = (["time", "depth", "YG", "XG"], np.full((2, 1, 2, 2), kh_zonal))
-    ds["Kh_meridional"] = (["time", "depth", "YG", "XG"], np.full((2, 1, 2, 2), kh_meridional))
-    Kh_zonal = Field("Kh_zonal", ds["Kh_zonal"], grid=grid, interp_method=XLinear)
-    Kh_meridional = Field("Kh_meridional", ds["Kh_meridional"], grid=grid, interp_method=XLinear)
     UV = VectorField("UV", U, V)
-    fieldset = FieldSet([U, V, UV, Kh_zonal, Kh_meridional])
+    fieldset = FieldSet([U, V, UV])
+    fieldset.add_constant_field("Kh_zonal", kh_zonal, mesh=mesh)
+    fieldset.add_constant_field("Kh_meridional", kh_meridional, mesh=mesh)
 
     npart = 100
     runtime = np.timedelta64(2, "h")

tests/test_fieldset.py

@@ -59,7 +59,6 @@ def test_fieldset_add_constant_field(fieldset):
     lat = ds["lat"].mean()
     lon = ds["lon"].mean()
 
-    pytest.xfail(reason="Not yet implemented interpolation.")
     assert fieldset.test_constant_field[time, z, lat, lon] == 1.0