diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb index ffe6f6a78..795743b86 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb @@ -147,11 +147,12 @@ "metadata": {}, "outputs": [], "source": [ - "u, v = fieldset.UV.eval(\n", - " np.array([0]), np.array([0]), np.array([20]), np.array([50]), apply_conversion=False\n", - ") # do not convert m/s to deg/s\n", + "u, v = fieldset.UV.eval(np.array([0]), np.array([0]), np.array([60]), np.array([50]))\n", + "u *= 1852 * 60 * np.cos(np.deg2rad(60)) # convert from degrees s^-1 to m s^-1\n", + "v *= 1852 * 60 # convert from degrees s^-1 to m s\n", "print(f\"(u, v) = ({u[0]:.3f}, {v[0]:.3f})\")\n", - "assert np.isclose(u, 1.0, atol=1e-3)" + "assert np.isclose(u, 1.0, atol=1e-3)\n", + "assert np.isclose(v, 0.0, atol=1e-3)" ] }, { diff --git a/docs/user_guide/examples/tutorial_unitconverters.ipynb b/docs/user_guide/examples/tutorial_unitconverters.ipynb index 4bcf19ede..e75a68642 100644 --- a/docs/user_guide/examples/tutorial_unitconverters.ipynb +++ b/docs/user_guide/examples/tutorial_unitconverters.ipynb @@ -15,7 +15,9 @@ "source": [ "In most applications, Parcels works with `spherical` meshes, where longitude and latitude are given in degrees, while depth is given in meters. But it is also possible to use `flat` meshes, where longitude and latitude are given in meters (note that the dimensions are then still called `longitude` and `latitude` for consistency reasons).\n", "\n", - "In all cases, velocities are given in m/s. Parcels seamlessly converts between meters and degrees, under the hood. For transparency, this guide explains how this works.\n" + "In all cases, velocities are given in m s-1. But for advection (such as `particles.dlat = v * particles.dt`), the velocity needs to be in degrees s-1. \n", + "\n", + "Parcels seamlessly converts between meters s-1 and degrees s-1, under the hood. For transparency, this guide explains how this works.\n" ] }, { @@ -130,7 +132,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Indeed, if we multiply the value of the U field with 1852 \\* 60 \\* cos(lat) (the number of meters in 1 degree of longitude) and the value of the V field with 1852 \\* 60 (the number of meters in 1 degree of latitude), we get the expected 1 m s-1 for `u` and `v`.\n" + "Indeed, if we multiply the value of `u` with 1852 \\* 60 \\* cos(lat) (the number of meters in 1 degree of longitude) and the value of `v` with 1852 \\* 60 (the number of meters in 1 degree of latitude), we get the expected 1 m s-1 for `u` and `v`.\n" ] }, { @@ -149,31 +151,6 @@ "assert np.isclose(v, 1.0)" ] }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can also interpolate the Field to these values directly by using the `eval()` method and setting `applyConversion=False`, as below\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(\n", - " fieldset.UV.eval(\n", - " time,\n", - " z,\n", - " lat,\n", - " lon,\n", - " apply_conversion=False,\n", - " )\n", - ")" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -185,7 +162,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Sampling `U` and `V` separately will _not_ convert to degrees s-1, so these velocities cannot be used directly for advection on spherical coordinates. This is one of the main reasons to always use the `UV` VectorField for velocity sampling in Parcels.\n" + "Sampling `U` and `V` separately will _not_ convert to degrees s-1, so these velocities _should not_ be used directly for advection on spherical coordinates. This is one of the main reasons to always use the `UV` VectorField for velocity sampling in Parcels.\n" ] }, { @@ -204,7 +181,7 @@ "source": [ "## Unit conversion for other fields such as diffusivity\n", "\n", - "For other fields such as diffusivity, Parcels does not apply similar unit conversions when using a `spherical` mesh. For example, if we define a diffusivity field with value 10 m2 s-1, Parcels will not convert this to degrees2 s-1 under the hood. \n", + "For other fields such as diffusivity, Parcels does not apply similar unit conversions when using a `spherical` mesh. For example, if we define a diffusivity field in m2 s-1, Parcels will _not_ convert this to degrees2 s-1 under the hood. \n", "\n", "If you want to work with diffusivity in degrees2 s-1 (for example to move particles using a random walk), you will have to convert this yourself in your kernel. \n", "\n", diff --git a/docs/user_guide/v4-migration.md b/docs/user_guide/v4-migration.md index af0056d2d..6e2f4743d 100644 --- a/docs/user_guide/v4-migration.md +++ b/docs/user_guide/v4-migration.md @@ -40,7 +40,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat - `Field.eval()` returns an array of floats instead of a single float (related to the vectorization) - `Field.eval()` does not throw OutOfBounds or other errors - The `NestedField` class has been removed. See the Nested Grids how-to guide for how to set up Nested Grids in v4. -- `applyConversion` has been renamed to `apply_conversion` and only works for VectorFields. Conversion of units should be handled in Kernels. +- `applyConversion` has been removed. Interpolation on VectorFields automatically converts from m/s to degrees/s for spherical meshes. Other conversion of units should be handled in Interpolators or Kernels. ## GridSet diff --git a/src/parcels/_core/field.py b/src/parcels/_core/field.py index 831946005..e404ad362 100644 --- a/src/parcels/_core/field.py +++ b/src/parcels/_core/field.py @@ -268,12 +268,10 @@ def __init__( else: self.vector_type = "2D" - # Setting the interpolation method dynamically if vector_interp_method is None: - self._vector_interp_method = None - else: - assert_same_function_signature(vector_interp_method, ref=ZeroInterpolator_Vector, context="Interpolation") - self._vector_interp_method = vector_interp_method + raise ValueError("vector_interp_method must be provided for VectorField initialization.") + assert_same_function_signature(vector_interp_method, ref=ZeroInterpolator_Vector, context="Interpolation") + self._vector_interp_method = vector_interp_method def __repr__(self): return vectorfield_repr(self) @@ -287,7 +285,7 @@ def vector_interp_method(self, method: Callable): assert_same_function_signature(method, ref=ZeroInterpolator_Vector, context="Interpolation") self._vector_interp_method = method - def eval(self, time: datetime, z, y, x, particles=None, apply_conversion=True): + def eval(self, time: datetime, z, y, x, particles=None): """Interpolate vectorfield values in space and time. Parameters @@ -300,10 +298,6 @@ def eval(self, time: datetime, z, y, x, particles=None, apply_conversion=True): particles : ParticleSet, optional If provided, used to associate results with particle indices and to update particle state and element indices. Defaults to None. - apply_conversion : bool, default True - If True and the underlying grid is spherical, the horizontal velocity - components (U, V) are converted from m s^-1 to degrees s^-1. - If False, raw interpolated values are returned. Returns ------- @@ -327,26 +321,7 @@ def eval(self, time: datetime, z, y, x, particles=None, apply_conversion=True): particle_positions, grid_positions = _get_positions(self.U, time, z, y, x, particles, _ei) - if self._vector_interp_method is None: - u = self.U._interp_method(particle_positions, grid_positions, self.U) - v = self.V._interp_method(particle_positions, grid_positions, self.V) - - if apply_conversion and self.U.grid._mesh == "spherical": - u /= 1852 * 60.0 * np.cos(np.deg2rad(y)) - v /= 1852 * 60.0 - - if "3D" in self.vector_type: - w = self.W._interp_method(particle_positions, grid_positions, self.W) - else: - w = 0.0 - else: - (u, v, w) = self._vector_interp_method(particle_positions, grid_positions, self) - - if apply_conversion: - if self.U.grid._mesh == "spherical": - meshJac = 1852 * 60.0 * np.cos(np.deg2rad(y)) - u = u / meshJac - v = v / meshJac + (u, v, w) = self._vector_interp_method(particle_positions, grid_positions, self) for vel in (u, v, w): _update_particle_states_interp_value(particles, vel) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 26350e182..f11241294 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -20,7 +20,14 @@ from parcels._logger import logger from parcels._reprs import fieldset_repr from parcels._typing import Mesh -from parcels.interpolators import UxPiecewiseConstantFace, UxPiecewiseLinearNode, XConstantField, XLinear +from parcels.interpolators import ( + Ux_Velocity, + UxPiecewiseConstantFace, + UxPiecewiseLinearNode, + XConstantField, + XLinear, + XLinear_Velocity, +) if TYPE_CHECKING: from parcels._core.basegrid import BaseGrid @@ -266,9 +273,11 @@ def from_fesom2(cls, ds: ux.UxDataset): if "W" in ds.data_vars: fields["W"] = Field("W", ds["W"], grid, _select_uxinterpolator(ds["U"])) - fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"]) + fields["UVW"] = VectorField( + "UVW", fields["U"], fields["V"], fields["W"], vector_interp_method=Ux_Velocity + ) else: - fields["UV"] = VectorField("UV", fields["U"], fields["V"]) + fields["UV"] = VectorField("UV", fields["U"], fields["V"], vector_interp_method=Ux_Velocity) for varname in set(ds.data_vars) - set(fields.keys()): fields[varname] = Field(varname, ds[varname], grid, _select_uxinterpolator(ds[varname])) @@ -350,11 +359,13 @@ def from_sgrid_conventions( if "U" in ds.data_vars and "V" in ds.data_vars: fields["U"] = Field("U", ds["U"], grid, XLinear) fields["V"] = Field("V", ds["V"], grid, XLinear) - fields["UV"] = VectorField("UV", fields["U"], fields["V"]) + fields["UV"] = VectorField("UV", fields["U"], fields["V"], vector_interp_method=XLinear_Velocity) if "W" in ds.data_vars: fields["W"] = Field("W", ds["W"], grid, XLinear) - fields["UVW"] = VectorField("UVW", fields["U"], fields["V"], fields["W"]) + fields["UVW"] = VectorField( + "UVW", fields["U"], fields["V"], fields["W"], vector_interp_method=XLinear_Velocity + ) for varname in set(ds.data_vars) - set(fields.keys()) - skip_vars: fields[varname] = Field(varname, ds[varname], grid, XLinear) diff --git a/src/parcels/_datasets/structured/generated.py b/src/parcels/_datasets/structured/generated.py index 5bcf00c9c..8faeec829 100644 --- a/src/parcels/_datasets/structured/generated.py +++ b/src/parcels/_datasets/structured/generated.py @@ -15,6 +15,7 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh="spherical"): max_lon = 180.0 if mesh == "spherical" else 1e6 + max_lat = 90.0 if mesh == "spherical" else 1e6 return xr.Dataset( {"U": (["time", "depth", "YG", "XG"], np.zeros(dims)), "V": (["time", "depth", "YG", "XG"], np.zeros(dims))}, @@ -25,7 +26,7 @@ def simple_UV_dataset(dims=(360, 2, 30, 4), maxdepth=1, mesh="spherical"): "YG": (["YG"], np.arange(dims[2]), {"axis": "Y", "c_grid_axis_shift": -0.5}), "XC": (["XC"], np.arange(dims[3]) + 0.5, {"axis": "X"}), "XG": (["XG"], np.arange(dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), - "lat": (["YG"], np.linspace(-90, 90, dims[2]), {"axis": "Y", "c_grid_axis_shift": -0.5}), + "lat": (["YG"], np.linspace(-max_lat, max_lat, dims[2]), {"axis": "Y", "c_grid_axis_shift": -0.5}), "lon": (["XG"], np.linspace(-max_lon, max_lon, dims[3]), {"axis": "X", "c_grid_axis_shift": -0.5}), }, ).pipe( diff --git a/src/parcels/interpolators.py b/src/parcels/interpolators.py index 4d4382abe..61b58650b 100644 --- a/src/parcels/interpolators.py +++ b/src/parcels/interpolators.py @@ -20,10 +20,12 @@ "CGrid_Velocity", "UxPiecewiseConstantFace", "UxPiecewiseLinearNode", + "Ux_Velocity", "XConstantField", "XFreeslip", "XLinear", "XLinearInvdistLandTracer", + "XLinear_Velocity", "XNearest", "XPartialslip", "ZeroInterpolator", @@ -146,6 +148,25 @@ def XConstantField( return field.data[0, 0, 0, 0].values +def XLinear_Velocity( + particle_positions: dict[str, float | np.ndarray], + grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]], + vectorfield: VectorField, +): + """Trilinear interpolation on a regular grid for VectorFields of velocity.""" + u = XLinear(particle_positions, grid_positions, vectorfield.U) + v = XLinear(particle_positions, grid_positions, vectorfield.V) + if vectorfield.grid._mesh == "spherical": + u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"])) + v /= 1852 * 60 + + if vectorfield.W: + w = XLinear(particle_positions, grid_positions, vectorfield.W) + else: + w = 0.0 + return u, v, w + + def CGrid_Velocity( particle_positions: dict[str, float | np.ndarray], grid_positions: dict[_XGRID_AXES, dict[str, int | float | np.ndarray]], @@ -268,9 +289,10 @@ def CGrid_Velocity( U = (1 - xsi) * U0 + xsi * U1 V = (1 - eta) * V0 + eta * V1 - meshJac = 1852 * 60.0 if grid._mesh == "spherical" else 1 - - jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac + if grid._mesh == "spherical": + jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * 1852 * 60.0 + else: + jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) u = ( (-(1 - eta) * U - (1 - xsi) * V) * px[0] @@ -288,6 +310,11 @@ def CGrid_Velocity( u = u.compute() v = v.compute() + if grid._mesh == "spherical": + conversion = 1852 * 60.0 * np.cos(np.deg2rad(particle_positions["lat"])) + u /= conversion + v /= conversion + # check whether the grid conversion has been applied correctly xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3] dlon = xx - lon @@ -682,3 +709,22 @@ def UxPiecewiseLinearNode( zk = field.grid.z.values[zi] zkp1 = field.grid.z.values[zi + 1] return (fzk * (zkp1 - z) + fzkp1 * (z - zk)) / (zkp1 - zk) # Linear interpolation in the vertical direction + + +def Ux_Velocity( + particle_positions: dict[str, float | np.ndarray], + grid_positions: dict[_UXGRID_AXES, dict[str, int | float | np.ndarray]], + vectorfield: VectorField, +): + """Interpolation kernel for Vectorfields of velocity on a UxGrid.""" + u = vectorfield.U._interp_method(particle_positions, grid_positions, vectorfield.U) + v = vectorfield.V._interp_method(particle_positions, grid_positions, vectorfield.V) + if vectorfield.grid._mesh == "spherical": + u /= 1852 * 60 * np.cos(np.deg2rad(particle_positions["lat"])) + v /= 1852 * 60 + + if "3D" in vectorfield.vector_type: + w = vectorfield.W._interp_method(particle_positions, grid_positions, vectorfield.W) + else: + w = 0.0 + return u, v, w diff --git a/tests-v3/test_fieldset_sampling.py b/tests-v3/test_fieldset_sampling.py index 1a3c68623..291c27b88 100644 --- a/tests-v3/test_fieldset_sampling.py +++ b/tests-v3/test_fieldset_sampling.py @@ -29,10 +29,6 @@ def SampleUV(particle, fieldset, time): # pragma: no cover (particle.u, particle.v) = fieldset.UV[time, particle.depth, particle.lat, particle.lon] -def SampleUVNoConvert(particle, fieldset, time): # pragma: no cover - (particle.u, particle.v) = fieldset.UV.eval(time, particle.depth, particle.lat, particle.lon, applyConversion=False) - - def SampleP(particle, fieldset, time): # pragma: no cover particle.p = fieldset.P[particle] @@ -439,24 +435,6 @@ def test_fieldset_sample_geographic(fieldset_geometric): assert np.allclose(pset.u, lat, rtol=1e-6) -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_sample_geographic_noconvert(fieldset_geometric): - """Sample a fieldset without conversion to geographic units.""" - npart = 120 - fieldset = fieldset_geometric - lon = np.linspace(-170, 170, npart) - lat = np.linspace(-80, 80, npart) - - pset = ParticleSet(fieldset, pclass=pclass(), lon=lon, lat=np.zeros(npart) + 70.0) - pset.execute(pset.Kernel(SampleUVNoConvert), endtime=1.0, dt=1.0) - assert np.allclose(pset.v, lon * 1000 * 1.852 * 60, rtol=1e-6) - - pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0) - pset.execute(pset.Kernel(SampleUVNoConvert), endtime=1.0, dt=1.0) - assert np.allclose(pset.u, lat * 1000 * 1.852 * 60, rtol=1e-6) - - @pytest.mark.v4alpha @pytest.mark.xfail(reason="GH1946") def test_fieldset_sample_geographic_polar(fieldset_geometric_polar): diff --git a/tests/test_advection.py b/tests/test_advection.py index ab882b15a..cc94ccd28 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -14,7 +14,7 @@ simple_UV_dataset, stommel_gyre_dataset, ) -from parcels.interpolators import CGrid_Velocity, XLinear +from parcels.interpolators import CGrid_Velocity, XLinear, XLinear_Velocity from parcels.kernels import ( AdvectionDiffusionEM, AdvectionDiffusionM1, @@ -35,13 +35,18 @@ def test_advection_zonal(mesh, npart=10): ds["U"].data[:] = 1.0 fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) - pset = ParticleSet(fieldset, lon=np.zeros(npart) + 20.0, lat=np.linspace(0, 80, npart)) - pset.execute(AdvectionRK4, runtime=np.timedelta64(2, "h"), dt=np.timedelta64(15, "m")) + runtime = 7200 + startlat = np.linspace(0, 80, npart) + startlon = 20.0 + np.zeros(npart) + pset = ParticleSet(fieldset, lon=startlon, lat=startlat) + pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) + expected_dlon = runtime if mesh == "spherical": - assert (np.diff(pset.lon) > 1.0e-4).all() - else: - assert (np.diff(pset.lon) < 1.0e-4).all() + expected_dlon /= 1852 * 60 * np.cos(np.deg2rad(pset.lat)) + + np.testing.assert_allclose(pset.lon - startlon, expected_dlon, atol=1e-5) + np.testing.assert_allclose(pset.lat, startlat, atol=1e-5) def test_advection_zonal_with_particlefile(tmp_store): @@ -88,17 +93,41 @@ def test_advection_zonal_periodic(): np.testing.assert_allclose(pset.lat, 0.5, atol=1e-5) -def test_horizontal_advection_in_3D_flow(npart=10): - """Flat 2D zonal flow that increases linearly with z from 0 m/s to 1 m/s.""" - ds = simple_UV_dataset(mesh="flat") +@pytest.mark.parametrize("mesh", ["spherical", "flat"]) +def test_advection_meridional(mesh, npart=10): + """All particles move the same in meridional direction, regardless of latitude.""" + ds = simple_UV_dataset(mesh=mesh) + ds["V"].data[:] = 1.0 + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) + + runtime = 7200 + startlat = np.linspace(0, 80, npart) + startlon = 20.0 + np.zeros(npart) + pset = ParticleSet(fieldset, lon=startlon, lat=startlat) + pset.execute(AdvectionRK4, runtime=runtime, dt=np.timedelta64(15, "m")) + + expected_dlat = runtime + if mesh == "spherical": + expected_dlat /= 1852 * 60 + + np.testing.assert_allclose(pset.lon, startlon, atol=1e-5) + np.testing.assert_allclose(pset.lat - startlat, expected_dlat, atol=1e-4) + + +@pytest.mark.parametrize("mesh", ["spherical", "flat"]) +def test_horizontal_advection_in_3D_flow(mesh, npart=10): + """2D zonal flow that increases linearly with z from 0 m/s to 1 m/s.""" + ds = simple_UV_dataset(mesh=mesh) ds["U"].data[:] = 1.0 ds["U"].data[:, 0, :, :] = 0.0 # Set U to 0 at the surface - fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") + fieldset = FieldSet.from_sgrid_conventions(ds, mesh=mesh) pset = ParticleSet(fieldset, lon=np.zeros(npart), lat=np.zeros(npart), z=np.linspace(0.1, 0.9, npart)) pset.execute(AdvectionRK4, runtime=np.timedelta64(2, "h"), dt=np.timedelta64(15, "m")) expected_lon = pset.z * pset.time + if mesh == "spherical": + expected_lon /= 1852 * 60 * np.cos(np.deg2rad(pset.lat)) np.testing.assert_allclose(pset.lon, expected_lon, atol=1.0e-1) @@ -184,10 +213,10 @@ def test_length1dimensions(u, v, w): # TODO: Refactor this test to be more read grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) - fields = [U, V, VectorField("UV", U, V)] + fields = [U, V, VectorField("UV", U, V, vector_interp_method=XLinear_Velocity)] if w: W = Field("W", ds["W"], grid, interp_method=XLinear) - fields.append(VectorField("UVW", U, V, W)) + fields.append(VectorField("UVW", U, V, W, vector_interp_method=XLinear_Velocity)) fieldset = FieldSet(fields) x0, y0, z0 = 2, 8, -4 @@ -207,7 +236,7 @@ def test_radialrotation(npart=10): grid = XGrid.from_dataset(ds, mesh="flat") U = parcels.Field("U", ds["U"], grid, interp_method=XLinear) V = parcels.Field("V", ds["V"], grid, interp_method=XLinear) - UV = parcels.VectorField("UV", U, V) + UV = parcels.VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = parcels.FieldSet([U, V, UV]) dt = np.timedelta64(30, "s") @@ -248,10 +277,10 @@ def test_moving_eddy(kernel, rtol): if kernel in [AdvectionRK2_3D, AdvectionRK4_3D]: # Using W to test 3D advection (assuming same velocity as V) W = Field("W", ds["V"], grid, interp_method=XLinear) - UVW = VectorField("UVW", U, V, W) + UVW = VectorField("UVW", U, V, W, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, W, UVW]) else: - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, UV]) if kernel in [AdvectionDiffusionEM, AdvectionDiffusionM1]: # Add zero diffusivity field for diffusion kernels @@ -299,7 +328,7 @@ def test_decaying_moving_eddy(kernel, rtol): grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, UV]) start_lon, start_lat = 10000, 10000 @@ -347,7 +376,7 @@ def test_stommelgyre_fieldset(kernel, rtol, grid_type): npart = 2 ds = stommel_gyre_dataset(grid_type=grid_type) grid = XGrid.from_dataset(ds, mesh="flat") - vector_interp_method = None if grid_type == "A" else CGrid_Velocity + vector_interp_method = XLinear_Velocity if grid_type == "A" else CGrid_Velocity U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) P = Field("P", ds["P"], grid, interp_method=XLinear) @@ -391,7 +420,7 @@ def test_peninsula_fieldset(kernel, rtol, grid_type): U = Field("U", ds["U"], grid, interp_method=XLinear) V = Field("V", ds["V"], grid, interp_method=XLinear) P = Field("P", ds["P"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) fieldset = FieldSet([U, V, P, UV]) dt = np.timedelta64(30, "m") diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 1a2515bcb..23b205cc9 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -13,7 +13,7 @@ from parcels._datasets.structured.generic import datasets as datasets_structured from parcels._datasets.structured.generic import datasets_sgrid from parcels._datasets.unstructured.generic import datasets as datasets_unstructured -from parcels.interpolators import XLinear +from parcels.interpolators import XLinear, XLinear_Velocity from tests import utils ds = datasets_structured["ds_2d_left"] @@ -25,7 +25,7 @@ def fieldset() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet( [U, V, UV], @@ -165,7 +165,6 @@ def test_fieldset_init_incompatible_calendars(): grid = XGrid.from_dataset(ds1, mesh="flat") U = Field("U", ds1["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds1["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) ds2 = ds.copy() ds2["time"] = ( @@ -177,7 +176,7 @@ def test_fieldset_init_incompatible_calendars(): incompatible_calendar = Field("test", ds2["data_g"], grid2, interp_method=XLinear) with pytest.raises(CalendarError, match="Expected field '.*' to have calendar compatible with datetime object"): - FieldSet([U, V, UV, incompatible_calendar]) + FieldSet([U, V, incompatible_calendar]) def test_fieldset_add_field_incompatible_calendars(fieldset): diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index 2f6df98a5..5b433c3e0 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -25,7 +25,7 @@ from parcels._core.utils.time import TimeInterval, timedelta_to_float from parcels._datasets.structured.generated import peninsula_dataset from parcels._datasets.structured.generic import datasets -from parcels.interpolators import XLinear +from parcels.interpolators import XLinear, XLinear_Velocity from parcels.kernels import AdvectionRK4 from tests.common_kernels import DoNothing @@ -37,7 +37,7 @@ def fieldset() -> FieldSet: # TODO v4: Move into a `conftest.py` file and remov grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, XLinear) V = Field("V", ds["V_A_grid"], grid, XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet( [U, V, UV], diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index b989356dc..5cab91ee3 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -23,7 +23,7 @@ from parcels._datasets.structured.generated import simple_UV_dataset from parcels._datasets.structured.generic import datasets as datasets_structured from parcels._datasets.unstructured.generic import datasets as datasets_unstructured -from parcels.interpolators import UxPiecewiseConstantFace, UxPiecewiseLinearNode, XLinear +from parcels.interpolators import Ux_Velocity, UxPiecewiseConstantFace, UxPiecewiseLinearNode, XLinear, XLinear_Velocity from parcels.kernels import AdvectionEE, AdvectionRK2, AdvectionRK4, AdvectionRK4_3D, AdvectionRK45 from tests.common_kernels import DoNothing from tests.utils import DEFAULT_PARTICLES @@ -35,7 +35,7 @@ def fieldset() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet([U, V, UV]) @@ -47,7 +47,7 @@ def fieldset_no_time_interval() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet([U, V, UV]) @@ -494,7 +494,7 @@ def test_uxstommelgyre_pset_execute(): grid=grid, interp_method=UxPiecewiseConstantFace, ) - UV = VectorField(name="UV", U=U, V=V) + UV = VectorField(name="UV", U=U, V=V, vector_interp_method=Ux_Velocity) fieldset = FieldSet([UV, UV.U, UV.V, P]) pset = ParticleSet( fieldset, @@ -540,7 +540,7 @@ def test_uxstommelgyre_multiparticle_pset_execute(): grid=grid, interp_method=UxPiecewiseConstantFace, ) - UVW = VectorField(name="UVW", U=U, V=V, W=W) + UVW = VectorField(name="UVW", U=U, V=V, W=W, vector_interp_method=Ux_Velocity) fieldset = FieldSet([UVW, UVW.U, UVW.V, UVW.W, P]) pset = ParticleSet( fieldset, @@ -579,7 +579,7 @@ def test_uxstommelgyre_pset_execute_output(): grid=grid, interp_method=UxPiecewiseConstantFace, ) - UV = VectorField(name="UV", U=U, V=V) + UV = VectorField(name="UV", U=U, V=V, vector_interp_method=Ux_Velocity) fieldset = FieldSet([UV, UV.U, UV.V, P]) pset = ParticleSet( fieldset, diff --git a/tests/test_particlesetview.py b/tests/test_particlesetview.py index 847878f54..d5d60161a 100644 --- a/tests/test_particlesetview.py +++ b/tests/test_particlesetview.py @@ -4,7 +4,7 @@ from parcels import Field, FieldSet, Particle, ParticleSet, Variable, VectorField, XGrid from parcels._core.statuscodes import StatusCode from parcels._datasets.structured.generic import datasets as datasets_structured -from parcels.interpolators import XLinear +from parcels.interpolators import XLinear, XLinear_Velocity @pytest.fixture @@ -13,7 +13,7 @@ def fieldset() -> FieldSet: grid = XGrid.from_dataset(ds, mesh="flat") U = Field("U", ds["U_A_grid"], grid, interp_method=XLinear) V = Field("V", ds["V_A_grid"], grid, interp_method=XLinear) - UV = VectorField("UV", U, V) + UV = VectorField("UV", U, V, vector_interp_method=XLinear_Velocity) return FieldSet([U, V, UV]) diff --git a/tests/test_uxarray_fieldset.py b/tests/test_uxarray_fieldset.py index cf3c56dd9..b6895b52d 100644 --- a/tests/test_uxarray_fieldset.py +++ b/tests/test_uxarray_fieldset.py @@ -13,6 +13,7 @@ ) from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels.interpolators import ( + Ux_Velocity, UxPiecewiseConstantFace, UxPiecewiseLinearNode, ) @@ -47,6 +48,7 @@ def uv_fesom_channel(ds_fesom_channel) -> VectorField: grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"], mesh="flat"), interp_method=UxPiecewiseConstantFace, ), + vector_interp_method=Ux_Velocity, ) return UV @@ -73,6 +75,7 @@ def uvw_fesom_channel(ds_fesom_channel) -> VectorField: grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"], mesh="flat"), interp_method=UxPiecewiseLinearNode, ), + vector_interp_method=Ux_Velocity, ) return UVW @@ -118,11 +121,12 @@ def test_fesom2_square_delaunay_uniform_z_coordinate_eval(): U=Field(name="U", data=ds.U, grid=grid, interp_method=UxPiecewiseConstantFace), V=Field(name="V", data=ds.V, grid=grid, interp_method=UxPiecewiseConstantFace), W=Field(name="W", data=ds.W, grid=grid, interp_method=UxPiecewiseLinearNode), + vector_interp_method=Ux_Velocity, ) P = Field(name="p", data=ds.p, grid=grid, interp_method=UxPiecewiseLinearNode) fieldset = FieldSet([UVW, P, UVW.U, UVW.V, UVW.W]) - (u, v, w) = fieldset.UVW.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0], apply_conversion=False) + (u, v, w) = fieldset.UVW.eval(time=[0.0], z=[1.0], y=[30.0], x=[30.0]) assert np.allclose([u.item(), v.item(), w.item()], [1.0, 1.0, 0.0], rtol=1e-3, atol=1e-6) assert np.isclose(