diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 5864b9288c..701f127bf7 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -24,7 +24,7 @@ from parcels._core.utils.time import is_compatible as datetime_is_compatible from parcels._core.warnings import FieldSetWarning from parcels._python import NOTSET, NotSetType -from parcels._reprs import fieldset_describe +from parcels._repr_utils import fieldset_describe from parcels.interpolators import ( XConstantField, ) diff --git a/src/parcels/_core/particle.py b/src/parcels/_core/particle.py index 6252fba047..2611ccdb92 100644 --- a/src/parcels/_core/particle.py +++ b/src/parcels/_core/particle.py @@ -8,7 +8,7 @@ from parcels._compat import _attrgetter_helper from parcels._core.statuscodes import StatusCode from parcels._core.utils.string import _assert_str_and_python_varname -from parcels._reprs import particleclass_repr, variable_repr +from parcels._repr_utils import particleclass_repr, variable_repr __all__ = ["Particle", "ParticleClass", "Variable"] _TO_WRITE_OPTIONS = [True, False] diff --git a/src/parcels/_core/particlefile.py b/src/parcels/_core/particlefile.py index 50bdf309ce..d12f4eb073 100644 --- a/src/parcels/_core/particlefile.py +++ b/src/parcels/_core/particlefile.py @@ -18,7 +18,7 @@ from parcels._core.particle import ParticleClass from parcels._core.particlesetview import ParticleSetView from parcels._core.utils.time import timedelta_to_float -from parcels._reprs import particlefile_repr +from parcels._repr_utils import particlefile_repr from parcels._typing import PathLike if TYPE_CHECKING: diff --git a/src/parcels/_core/particlesetview.py b/src/parcels/_core/particlesetview.py index 2ad36d4177..030ac10947 100644 --- a/src/parcels/_core/particlesetview.py +++ b/src/parcels/_core/particlesetview.py @@ -1,6 +1,6 @@ import numpy as np -from parcels._reprs import particlesetview_repr +from parcels._repr_utils import particlesetview_repr class ParticleSetView: diff --git a/src/parcels/_core/utils/time.py b/src/parcels/_core/utils/time.py index eee0fe164a..8a20bd1256 100644 --- a/src/parcels/_core/utils/time.py +++ b/src/parcels/_core/utils/time.py @@ -6,7 +6,7 @@ import cftime import numpy as np -from parcels._reprs import timeinterval_repr +from parcels._repr_utils import timeinterval_repr if TYPE_CHECKING: from parcels._typing import TimeLike diff --git a/src/parcels/_reprs.py b/src/parcels/_repr_utils.py similarity index 100% rename from src/parcels/_reprs.py rename to src/parcels/_repr_utils.py diff --git a/tests-v3/test_advection.py b/tests-v3/test_advection.py deleted file mode 100644 index 3d8f06bac3..0000000000 --- a/tests-v3/test_advection.py +++ /dev/null @@ -1,114 +0,0 @@ -import numpy as np -import pytest -import xarray as xr - -from parcels import ( - AdvectionAnalytical, - AdvectionDiffusionEM, - AdvectionDiffusionM1, - AdvectionEE, - AdvectionRK4, - AdvectionRK45, - FieldSet, - Particle, - ParticleSet, -) -from tests.utils import TEST_DATA - -kernel = { - "EE": AdvectionEE, - "RK4": AdvectionRK4, - "RK45": AdvectionRK45, - "AA": AdvectionAnalytical, - "AdvDiffEM": AdvectionDiffusionEM, - "AdvDiffM1": AdvectionDiffusionM1, -} - - -@pytest.fixture -def lon(): - xdim = 200 - return np.linspace(-170, 170, xdim, dtype=np.float32) - - -@pytest.fixture -def lat(): - ydim = 100 - return np.linspace(-80, 80, ydim, dtype=np.float32) - - -@pytest.fixture -def depth(): - zdim = 2 - return np.linspace(0, 30, zdim, dtype=np.float32) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="When refactoring fieldfilebuffer croco support was dropped. This will be fixed in v4.") -def test_advection_2DCROCO(): - fieldset = FieldSet.from_modulefile(TEST_DATA / "fieldset_CROCO2D.py") - - runtime = 1e4 - X = np.array([40e3, 80e3, 120e3]) - Y = np.ones(X.size) * 100e3 - Z = np.zeros(X.size) - pset = ParticleSet(fieldset=fieldset, pclass=Particle, lon=X, lat=Y, depth=Z) - - pset.execute([AdvectionRK4], runtime=runtime, dt=100) - assert np.allclose(pset.depth, Z.flatten(), atol=1e-3) - assert np.allclose(pset.lon_nextloop, [x + runtime for x in X], atol=1e-3) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_analyticalAgrid(): - lon = np.arange(0, 15, dtype=np.float32) - lat = np.arange(0, 15, dtype=np.float32) - U = np.ones((lat.size, lon.size), dtype=np.float32) - V = np.ones((lat.size, lon.size), dtype=np.float32) - fieldset = FieldSet.from_data({"U": U, "V": V}, {"lon": lon, "lat": lat}, mesh="flat") - pset = ParticleSet(fieldset, pclass=Particle, lon=1, lat=1) - - with pytest.raises(NotImplementedError): - pset.execute(AdvectionAnalytical, runtime=1) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1927") -@pytest.mark.parametrize("u", [1, -0.2, -0.3, 0]) -@pytest.mark.parametrize("v", [1, -0.3, 0, -1]) -@pytest.mark.parametrize("w", [None, 1, -0.3, 0, -1]) -@pytest.mark.parametrize("direction", [1, -1]) -def test_uniform_analytical(u, v, w, direction, tmp_zarrfile): - lon = np.arange(0, 15, dtype=np.float32) - lat = np.arange(0, 15, dtype=np.float32) - if w is not None: - depth = np.arange(0, 40, 2, dtype=np.float32) - U = u * np.ones((depth.size, lat.size, lon.size), dtype=np.float32) - V = v * np.ones((depth.size, lat.size, lon.size), dtype=np.float32) - W = w * np.ones((depth.size, lat.size, lon.size), dtype=np.float32) - fieldset = FieldSet.from_data({"U": U, "V": V, "W": W}, {"lon": lon, "lat": lat, "depth": depth}, mesh="flat") - fieldset.W.interp_method = "cgrid_velocity" - else: - U = u * np.ones((lat.size, lon.size), dtype=np.float32) - V = v * np.ones((lat.size, lon.size), dtype=np.float32) - fieldset = FieldSet.from_data({"U": U, "V": V}, {"lon": lon, "lat": lat}, mesh="flat") - fieldset.U.interp_method = "cgrid_velocity" - fieldset.V.interp_method = "cgrid_velocity" - - x0, y0, z0 = 6.1, 6.2, 20 - pset = ParticleSet(fieldset, pclass=Particle, lon=x0, lat=y0, depth=z0) - - outfile = pset.ParticleFile(name=tmp_zarrfile, outputdt=1, chunks=(1, 1)) - pset.execute(AdvectionAnalytical, runtime=4, dt=direction, output_file=outfile) - assert np.abs(pset.lon - x0 - pset.time * u) < 1e-6 - assert np.abs(pset.lat - y0 - pset.time * v) < 1e-6 - if w is not None: - assert np.abs(pset.depth - z0 - pset.time * w) < 1e-4 - - ds = xr.open_zarr(tmp_zarrfile) - times = (direction * ds["time"][:]).values.astype("timedelta64[s]")[0] - timeref = np.arange(1, 5).astype("timedelta64[s]") - assert np.allclose(times, timeref, atol=np.timedelta64(1, "ms")) - lons = ds["lon"][:].values - assert np.allclose(lons, x0 + direction * u * np.arange(1, 5)) diff --git a/tests-v3/test_examples.py b/tests-v3/test_examples.py deleted file mode 100644 index de42fbb86f..0000000000 --- a/tests-v3/test_examples.py +++ /dev/null @@ -1,48 +0,0 @@ -import os -import runpy -import shutil -import sys -import time -from pathlib import Path - -import pytest - -example_folder = (Path(__file__).parent / "../docs/examples").resolve() -example_fnames = [path.name for path in example_folder.glob("*.py")] - - -@pytest.fixture(autouse=True) -def cleanup_generated_data_files(): - """Clean up generated data files from test run. - - Records current folder contents before test, and cleans up any generated `.nc` files - and `.zarr` folders afterwards. For safety this is non-recursive. This function is - only necessary as the scripts being run aren't native pytest tests, so they don't - have access to the `tmpdir` fixture. - - """ - folder_contents = os.listdir() - yield - time.sleep(0.1) # Buffer so that files are closed before we try to delete them. - for fname in os.listdir(): - if fname in folder_contents: - continue - if not (fname.endswith(".nc") or fname.endswith(".zarr")): - continue - - path = Path(fname) - if path.is_dir(): - shutil.rmtree(path) - else: - path.unlink() - print(f"Removed {path}") - - -@pytest.mark.parametrize("example_fname", example_fnames) -def test_example_script(example_fname): - script = str(example_folder / example_fname) - - # Clear sys.argv, otherwise pytest pollutes it with its own arguments. - sys.argv = [sys.argv[0]] - - runpy.run_path(script, run_name="__main__") diff --git a/tests-v3/test_fieldset.py b/tests-v3/test_fieldset.py deleted file mode 100644 index 69295cf871..0000000000 --- a/tests-v3/test_fieldset.py +++ /dev/null @@ -1,417 +0,0 @@ -from datetime import timedelta - -import numpy as np -import pytest -import xarray as xr - -from parcels import ( - AdvectionRK4, - AdvectionRK4_3D, - FieldSet, - Particle, - ParticleSet, - Variable, -) -from parcels.field import VectorField -from parcels.tools.converters import GeographicPolar, Unity -from tests.utils import TEST_DATA - - -def generate_fieldset_data(xdim, ydim, zdim=1, tdim=1): - lon = np.linspace(0.0, 10.0, xdim, dtype=np.float32) - lat = np.linspace(0.0, 10.0, ydim, dtype=np.float32) - depth = np.zeros(zdim, dtype=np.float32) - time = np.zeros(tdim, dtype=np.float64) - if zdim == 1 and tdim == 1: - U, V = np.meshgrid(lon, lat) - dimensions = {"lat": lat, "lon": lon} - else: - U = np.ones((tdim, zdim, ydim, xdim)) - V = np.ones((tdim, zdim, ydim, xdim)) - dimensions = {"lat": lat, "lon": lon, "depth": depth, "time": time} - data = {"U": np.array(U, dtype=np.float32), "V": np.array(V, dtype=np.float32)} - - return (data, dimensions) - - -def to_xarray_dataset(data: dict[str, np.array], dimensions: dict[str, np.array]) -> xr.Dataset: - assert len(dimensions) in [2, 4], "this function only deals with output from generate_fieldset_data()" - - if len(dimensions) == 4: - return xr.Dataset( - { - "U": (["time", "depth", "lat", "lon"], data["U"]), - "V": (["time", "depth", "lat", "lon"], data["V"]), - }, - coords=dimensions, - ) - - return xr.Dataset( - { - "U": (["lat", "lon"], data["U"]), - "V": (["lat", "lon"], data["V"]), - }, - coords=dimensions, - ) - - -@pytest.mark.v4remove -@pytest.mark.xfail(reason="GH1946") -@pytest.fixture -def multifile_fieldset(tmp_path): - stem = "test_subsets" - - timestamps = np.arange(0, 4, 1) * 86400.0 - timestamps = np.expand_dims(timestamps, 1) - - ufiles = [] - vfiles = [] - for index, timestamp in enumerate(timestamps): - data, dimensions = generate_fieldset_data(100, 100) - path = tmp_path / f"{stem}_{index}.nc" - to_xarray_dataset(data, dimensions).pipe(assign_dataset_timestamp_dim, timestamp).to_netcdf(path) - ufiles.append(path) - vfiles.append(path) - - files = {"U": ufiles, "V": vfiles} - variables = {"U": "U", "V": "V"} - dimensions = {"lon": "lon", "lat": "lat", "time": "time"} - return FieldSet.from_netcdf(files, variables, dimensions) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_from_modulefile(): - nemo_fname = str(TEST_DATA / "fieldset_nemo.py") - nemo_error_fname = str(TEST_DATA / "fieldset_nemo_error.py") - - fieldset = FieldSet.from_modulefile(nemo_fname) - - fieldset = FieldSet.from_modulefile(nemo_fname) - assert fieldset.U.grid.lon.shape[1] == 21 - - with pytest.raises(IOError): - FieldSet.from_modulefile(nemo_error_fname) - - FieldSet.from_modulefile(nemo_error_fname, modulename="random_function_name") - - with pytest.raises(IOError): - FieldSet.from_modulefile(nemo_error_fname, modulename="none_returning_function") - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_field_from_netcdf_fieldtypes(): - filenames = { - "varU": { - "lon": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "lat": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "data": str(TEST_DATA / "Uu_eastward_nemo_cross_180lon.nc"), - }, - "varV": { - "lon": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "lat": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "data": str(TEST_DATA / "Vv_eastward_nemo_cross_180lon.nc"), - }, - } - variables = {"varU": "U", "varV": "V"} - dimensions = {"lon": "glamf", "lat": "gphif"} - - # first try without setting fieldtype - fset = FieldSet.from_nemo(filenames, variables, dimensions) - assert isinstance(fset.varU.units, Unity) - - # now try with setting fieldtype - fset = FieldSet.from_nemo(filenames, variables, dimensions, fieldtype={"varU": "U", "varV": "V"}) - assert isinstance(fset.varU.units, GeographicPolar) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_from_agrid_dataset(): - filenames = { - "lon": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "lat": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "data": str(TEST_DATA / "Uu_eastward_nemo_cross_180lon.nc"), - } - variable = {"U": "U"} - dimensions = {"lon": "glamf", "lat": "gphif"} - FieldSet.from_a_grid_dataset(filenames, variable, dimensions) - - -@pytest.mark.v4remove -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_from_cgrid_interpmethod(): - filenames = { - "lon": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "lat": str(TEST_DATA / "mask_nemo_cross_180lon.nc"), - "data": str(TEST_DATA / "Uu_eastward_nemo_cross_180lon.nc"), - } - variable = "U" - dimensions = {"lon": "glamf", "lat": "gphif"} - - with pytest.raises(TypeError): - # should fail because FieldSet.from_c_grid_dataset does not support interp_method - FieldSet.from_c_grid_dataset(filenames, variable, dimensions, interp_method="partialslip") - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("calltype", ["from_nemo"]) -def test_illegal_dimensionsdict(calltype): - with pytest.raises(NameError): - if calltype == "from_data": - data, dimensions = generate_fieldset_data(10, 10) - dimensions["test"] = None - FieldSet.from_data(data, dimensions) - elif calltype == "from_nemo": - fname = str(TEST_DATA / "mask_nemo_cross_180lon.nc") - filenames = {"dx": fname, "mesh_mask": fname} - variables = {"dx": "e1u"} - dimensions = {"lon": "glamu", "lat": "gphiu", "test": "test"} - FieldSet.from_nemo(filenames, variables, dimensions) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("gridtype", ["A", "C"]) -def test_fieldset_dimlength1_cgrid(gridtype): - fieldset = FieldSet.from_data({"U": 0, "V": 0}, {"lon": 0, "lat": 0}) # TODO : Remove from_data - if gridtype == "C": - fieldset.U.interp_method = "cgrid_velocity" - fieldset.V.interp_method = "cgrid_velocity" - try: - fieldset._check_complete() - success = True if gridtype == "A" else False - except NotImplementedError: - success = True if gridtype == "C" else False - assert success - - -def assign_dataset_timestamp_dim(ds, timestamp): - """Expand dim to 'time' and assign timestamp.""" - ds.expand_dims("time") - ds["time"] = timestamp - return ds - - -def addConst(particle, fieldset, time): # pragma: no cover - particle.lon = particle.lon + fieldset.movewest + fieldset.moveeast - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_constant(): - data, dimensions = generate_fieldset_data(100, 100) - fieldset = FieldSet.from_data(data, dimensions) # TODO : Remove from_data - westval = -0.2 - eastval = 0.3 - fieldset.add_constant("movewest", westval) - fieldset.add_constant("moveeast", eastval) - assert fieldset.movewest == westval - - pset = ParticleSet.from_line(fieldset, size=1, pclass=Particle, start=(0.5, 0.5), finish=(0.5, 0.5)) - pset.execute(pset.Kernel(addConst), dt=1, runtime=1) - assert abs(pset.lon[0] - (0.5 + westval + eastval)) < 1e-4 - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("swapUV", [False, True]) -def test_vector_fields(swapUV): - lon = np.linspace(0.0, 10.0, 12, dtype=np.float32) - lat = np.linspace(0.0, 10.0, 10, dtype=np.float32) - U = np.ones((10, 12), dtype=np.float32) - V = np.zeros((10, 12), dtype=np.float32) - data = {"U": U, "V": V} - dimensions = {"U": {"lat": lat, "lon": lon}, "V": {"lat": lat, "lon": lon}} - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") # TODO : Remove from_data - if swapUV: # we test that we can freely edit whatever UV field - UV = VectorField("UV", fieldset.V, fieldset.U) - fieldset.add_vector_field(UV) - - pset = ParticleSet.from_line(fieldset, size=1, pclass=Particle, start=(0.5, 0.5), finish=(0.5, 0.5)) - pset.execute(AdvectionRK4, dt=1, runtime=2) - if swapUV: - assert abs(pset.lon[0] - 0.5) < 1e-9 - assert abs(pset.lat[0] - 1.5) < 1e-9 - else: - assert abs(pset.lon[0] - 1.5) < 1e-9 - assert abs(pset.lat[0] - 0.5) < 1e-9 - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946, originated in GH938") -def test_add_second_vector_field(): - lon = np.linspace(0.0, 10.0, 12, dtype=np.float32) - lat = np.linspace(0.0, 10.0, 10, dtype=np.float32) - U = np.ones((10, 12), dtype=np.float32) - V = np.zeros((10, 12), dtype=np.float32) - data = {"U": U, "V": V} - dimensions = {"U": {"lat": lat, "lon": lon}, "V": {"lat": lat, "lon": lon}} - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") # TODO : Remove from_data - - data2 = {"U2": U, "V2": V} - dimensions2 = {"lon": [ln + 0.1 for ln in lon], "lat": [lt - 0.1 for lt in lat]} - fieldset2 = FieldSet.from_data(data2, dimensions2, mesh="flat") # TODO : Remove from_data - - UV2 = VectorField("UV2", fieldset2.U2, fieldset2.V2) - fieldset.add_vector_field(UV2) - - def SampleUV2(particle, fieldset, time): # pragma: no cover - u, v = fieldset.UV2[time, particle.depth, particle.lat, particle.lon] - particle.dlon += u * particle.dt - particle.dlat += v * particle.dt - - pset = ParticleSet(fieldset, pclass=Particle, lon=0.5, lat=0.5) - pset.execute(AdvectionRK4 + pset.Kernel(SampleUV2), dt=1, runtime=2) - - assert abs(pset.lon[0] - 2.5) < 1e-9 - assert abs(pset.lat[0] - 0.5) < 1e-9 - - -@pytest.mark.v4remove -@pytest.mark.xfail(reason="time_periodic removed in v4") -@pytest.mark.parametrize("use_xarray", [True, False]) -@pytest.mark.parametrize("time_periodic", [86400.0, False]) -@pytest.mark.parametrize("dt_sign", [-1, 1]) -def test_periodic(use_xarray, time_periodic, dt_sign): - lon = np.array([0, 1], dtype=np.float32) - lat = np.array([0, 1], dtype=np.float32) - depth = np.array([0, 1], dtype=np.float32) - tsize = 24 * 60 + 1 - period = 86400 - time = np.linspace(0, period, tsize, dtype=np.float64) - - def temp_func(time): - return 20 + 2 * np.sin(time * 2 * np.pi / period) - - temp_vec = temp_func(time) - - U = np.zeros((tsize, 2, 2, 2), dtype=np.float32) - V = np.zeros((tsize, 2, 2, 2), dtype=np.float32) - V[:, 0, 0, 0] = 1e-5 - W = np.zeros((tsize, 2, 2, 2), dtype=np.float32) - temp = np.zeros((tsize, 2, 2, 2), dtype=np.float32) - temp[:, :, :, :] = temp_vec - D = np.ones((2, 2), dtype=np.float32) # adding non-timevarying field - - full_dims = {"lon": lon, "lat": lat, "depth": depth, "time": time} - dimensions = {"U": full_dims, "V": full_dims, "W": full_dims, "temp": full_dims, "D": {"lon": lon, "lat": lat}} - if use_xarray: - coords = {"lat": lat, "lon": lon, "depth": depth, "time": time} - variables = {"U": "Uxr", "V": "Vxr", "W": "Wxr", "temp": "Txr", "D": "Dxr"} - dimnames = {"lon": "lon", "lat": "lat", "depth": "depth", "time": "time"} - ds = xr.Dataset( - { - "Uxr": xr.DataArray(U, coords=coords, dims=("time", "depth", "lat", "lon")), - "Vxr": xr.DataArray(V, coords=coords, dims=("time", "depth", "lat", "lon")), - "Wxr": xr.DataArray(W, coords=coords, dims=("time", "depth", "lat", "lon")), - "Txr": xr.DataArray(temp, coords=coords, dims=("time", "depth", "lat", "lon")), - "Dxr": xr.DataArray(D, coords={"lat": lat, "lon": lon}, dims=("lat", "lon")), - } - ) - fieldset = FieldSet.from_xarray_dataset( - ds, - variables, - {"U": dimnames, "V": dimnames, "W": dimnames, "temp": dimnames, "D": {"lon": "lon", "lat": "lat"}}, - time_periodic=time_periodic, - allow_time_extrapolation=True, - ) - else: - data = {"U": U, "V": V, "W": W, "temp": temp, "D": D} - fieldset = FieldSet.from_data( - data, dimensions, mesh="flat", time_periodic=time_periodic, allow_time_extrapolation=True - ) # TODO : Remove from_data - - def sampleTemp(particle, fieldset, time): # pragma: no cover - particle.temp = fieldset.temp[time, particle.depth, particle.lat, particle.lon] - # test if we can interpolate UV and UVW together - (particle.u1, particle.v1) = fieldset.UV[time, particle.depth, particle.lat, particle.lon] - (particle.u2, particle.v2, w_) = fieldset.UVW[time, particle.depth, particle.lat, particle.lon] - # test if we can sample a non-timevarying field too - particle.d = fieldset.D[0, 0, particle.lat, particle.lon] - - MyParticle = Particle.add_variables( - [ - Variable("temp", dtype=np.float32, initial=20.0), - Variable("u1", dtype=np.float32, initial=0.0), - Variable("u2", dtype=np.float32, initial=0.0), - Variable("v1", dtype=np.float32, initial=0.0), - Variable("v2", dtype=np.float32, initial=0.0), - Variable("d", dtype=np.float32, initial=0.0), - ] - ) - - pset = ParticleSet(fieldset, pclass=MyParticle, lon=[0.5], lat=[0.5], depth=[0.5]) - pset.execute( - AdvectionRK4_3D + pset.Kernel(sampleTemp), runtime=timedelta(hours=51), dt=timedelta(hours=dt_sign * 1) - ) - - if time_periodic is not False: - t = pset.time[0] - temp_theo = temp_func(t) - elif dt_sign == 1: - temp_theo = temp_vec[-1] - elif dt_sign == -1: - temp_theo = temp_vec[0] - assert np.allclose(temp_theo, pset.temp[0], atol=1e-5) - assert np.allclose(pset.u1[0], pset.u2[0]) - assert np.allclose(pset.v1[0], pset.v2[0]) - assert np.allclose(pset.d[0], 1.0) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_from_data_gridtypes(): - """Simple test for fieldset initialisation from data.""" - xdim, ydim, zdim = 20, 10, 4 - - lon = np.linspace(0.0, 10.0, xdim, dtype=np.float32) - lat = np.linspace(0.0, 10.0, ydim, dtype=np.float32) - depth = np.linspace(0.0, 1.0, zdim, dtype=np.float32) - depth_s = np.ones((zdim, ydim, xdim)) - U = np.ones((zdim, ydim, xdim)) - V = np.ones((zdim, ydim, xdim)) - dimensions = {"lat": lat, "lon": lon, "depth": depth} - data = {"U": np.array(U, dtype=np.float32), "V": np.array(V, dtype=np.float32)} - lonm, latm = np.meshgrid(lon, lat) - for k in range(zdim): - data["U"][k, :, :] = lonm * (depth[k] + 1) + 0.1 - depth_s[k, :, :] = depth[k] - - # Rectilinear Z grid - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") # TODO : Remove from_data - pset = ParticleSet(fieldset, Particle, [0, 0], [0, 0], [0, 0.4]) - pset.execute(AdvectionRK4, runtime=1.5, dt=0.5) - plon = pset.lon - plat = pset.lat - # sol of dx/dt = (init_depth+1)*x+0.1; x(0)=0 - assert np.allclose(plon, [0.17173462592827032, 0.2177736932123214]) - assert np.allclose(plat, [1, 1]) - - # Rectilinear S grid - dimensions["depth"] = depth_s - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") # TODO : Remove from_data - pset = ParticleSet(fieldset, Particle, [0, 0], [0, 0], [0, 0.4]) - pset.execute(AdvectionRK4, runtime=1.5, dt=0.5) - assert np.allclose(plon, pset.lon) - assert np.allclose(plat, pset.lat) - - # Curvilinear Z grid - dimensions["lon"] = lonm - dimensions["lat"] = latm - dimensions["depth"] = depth - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") # TODO : Remove from_data - pset = ParticleSet(fieldset, Particle, [0, 0], [0, 0], [0, 0.4]) - pset.execute(AdvectionRK4, runtime=1.5, dt=0.5) - assert np.allclose(plon, pset.lon) - assert np.allclose(plat, pset.lat) - - # Curvilinear S grid - dimensions["depth"] = depth_s - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") # TODO : Remove from_data - pset = ParticleSet(fieldset, Particle, [0, 0], [0, 0], [0, 0.4]) - pset.execute(AdvectionRK4, runtime=1.5, dt=0.5) - assert np.allclose(plon, pset.lon) - assert np.allclose(plat, pset.lat) diff --git a/tests-v3/test_fieldset_sampling.py b/tests-v3/test_fieldset_sampling.py deleted file mode 100644 index 291c27b885..0000000000 --- a/tests-v3/test_fieldset_sampling.py +++ /dev/null @@ -1,817 +0,0 @@ -import os -from datetime import timedelta -from math import cos, pi - -import numpy as np -import pytest -import xarray as xr - -from parcels import ( - AdvectionRK4, - AdvectionRK4_3D, - Field, - FieldSet, - Geographic, - Particle, - ParticleSet, - Variable, -) -from tests.utils import create_fieldset_global - - -def pclass(): - return Particle.add_variables( - [Variable("u", dtype=np.float32), Variable("v", dtype=np.float32), Variable("p", dtype=np.float32)] - ) - - -def SampleUV(particle, fieldset, time): # pragma: no cover - (particle.u, particle.v) = fieldset.UV[time, particle.depth, particle.lat, particle.lon] - - -def SampleP(particle, fieldset, time): # pragma: no cover - particle.p = fieldset.P[particle] - - -@pytest.fixture -def fieldset(): - return create_fieldset_global() - - -def create_fieldset_geometric(xdim=200, ydim=100): - """Standard earth fieldset with U and V equivalent to lon/lat in m.""" - lon = np.linspace(-180, 180, xdim, dtype=np.float32) - lat = np.linspace(-90, 90, ydim, dtype=np.float32) - V, U = np.meshgrid(lon, lat) - U *= 1000.0 * 1.852 * 60.0 - V *= 1000.0 * 1.852 * 60.0 - data = {"U": U, "V": V} - dimensions = {"lon": lon, "lat": lat} - fieldset = FieldSet.from_data(data, dimensions) - fieldset.U.units = Geographic() - fieldset.V.units = Geographic() - return fieldset - - -@pytest.fixture -def fieldset_geometric(): - return create_fieldset_geometric() - - -def create_fieldset_geometric_polar(xdim=200, ydim=100): - """Standard earth fieldset with U and V equivalent to lon/lat in m - and the inversion of the pole correction applied to U. - """ - lon = np.linspace(-180, 180, xdim, dtype=np.float32) - lat = np.linspace(-90, 90, ydim, dtype=np.float32) - V, U = np.meshgrid(lon, lat) - # Apply inverse of pole correction to U - for i, y in enumerate(lat): - U[i, :] *= cos(y * pi / 180) - U *= 1000.0 * 1.852 * 60.0 - V *= 1000.0 * 1.852 * 60.0 - data = {"U": U, "V": V} - dimensions = {"lon": lon, "lat": lat} - return FieldSet.from_data(data, dimensions, mesh="spherical") - - -@pytest.fixture -def fieldset_geometric_polar(): - return create_fieldset_geometric_polar() - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_sample(fieldset): - """Sample the fieldset using indexing notation.""" - xdim, ydim = 120, 80 - lon = np.linspace(-170, 170, xdim, dtype=np.float32) - lat = np.linspace(-80, 80, ydim, dtype=np.float32) - v_s = np.array([fieldset.UV[0, 0.0, 70.0, x][1] for x in lon]) - u_s = np.array([fieldset.UV[0, 0.0, y, -45.0][0] for y in lat]) - assert np.allclose( - v_s, lon, rtol=1e-5 - ) # Tolerances were rtol=1e-7, increased due to numpy v2 float32 changes (see #1603) - assert np.allclose(u_s, lat, rtol=1e-5) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_sample_eval(fieldset): - """Sample the fieldset using the explicit eval function.""" - xdim, ydim = 60, 60 - lon = np.linspace(-170, 170, xdim, dtype=np.float32) - lat = np.linspace(-80, 80, ydim, dtype=np.float32) - v_s = np.array([fieldset.UV.eval(0, 0.0, 70.0, x)[1] for x in lon]) - u_s = np.array([fieldset.UV.eval(0, 0.0, y, 0.0)[0] for y in lat]) - assert np.allclose( - v_s, lon, rtol=1e-5 - ) # Tolerances were rtol=1e-7, increased due to numpy v2 float32 changes (see #1603) - assert np.allclose(u_s, lat, rtol=1e-5) - - -@pytest.mark.v4remove -@pytest.mark.xfail(reason="Test is directly testing adding the halo. This test should either be adapted or removed.") -def test_fieldset_polar_with_halo(fieldset_geometric_polar): - fieldset_geometric_polar.add_periodic_halo(zonal=5) - pset = ParticleSet(fieldset_geometric_polar, pclass=pclass(), lon=0, lat=0) - pset.execute(runtime=1) - assert pset.lon[0] == 0.0 - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("zdir", [-1, 1]) -def test_verticalsampling(zdir): - dims = (4, 2, 2) - dimensions = { - "lon": np.linspace(0.0, 1.0, dims[2], dtype=np.float32), - "lat": np.linspace(0.0, 1.0, dims[1], dtype=np.float32), - "depth": np.linspace(0.0, 1 * zdir, dims[0], dtype=np.float32), - } - data = {"U": np.zeros(dims, dtype=np.float32), "V": np.zeros(dims, dtype=np.float32)} - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - pset = ParticleSet(fieldset, pclass=Particle, lon=0, lat=0, depth=0.7 * zdir) - pset.execute(AdvectionRK4, dt=1.0, runtime=1.0) - zi, yi, xi = fieldset.U.unravel_index(pset[0].ei) - assert zi == [2] - - -def test_pset_from_field(): - xdim = 10 - ydim = 20 - npart = 10000 - - np.random.seed(123456) - dimensions = { - "lon": np.linspace(0.0, 1.0, xdim, dtype=np.float32), - "lat": np.linspace(0.0, 1.0, ydim, dtype=np.float32), - } - startfield = np.ones((ydim, xdim), dtype=np.float32) - for x in range(xdim): - startfield[:, x] = x - data = { - "U": np.zeros((ydim, xdim), dtype=np.float32), - "V": np.zeros((ydim, xdim), dtype=np.float32), - "start": startfield, - } - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - - densfield = Field( - name="densfield", - data=np.zeros((ydim + 1, xdim + 1), dtype=np.float32), - lon=np.linspace(-1.0 / (xdim * 2), 1.0 + 1.0 / (xdim * 2), xdim + 1, dtype=np.float32), - lat=np.linspace(-1.0 / (ydim * 2), 1.0 + 1.0 / (ydim * 2), ydim + 1, dtype=np.float32), - ) - - fieldset.add_field(densfield) - pset = ParticleSet.from_field(fieldset, size=npart, pclass=Particle, start_field=fieldset.start) - pdens = np.histogram2d(pset.lat, pset.lon, bins=[np.linspace(0.0, 1.0, ydim + 1), np.linspace(0.0, 1.0, xdim + 1)])[ - 0 - ] - assert np.allclose(pdens / sum(pdens.flatten()), startfield / sum(startfield.flatten()), atol=1e-2) - - -@pytest.mark.v4alpha -def test_nearest_neighbor_interpolation2D(): - npart = 81 - dims = (2, 2) - dimensions = { - "lon": np.linspace(0.0, 1.0, dims[1], dtype=np.float32), - "lat": np.linspace(0.0, 1.0, dims[0], dtype=np.float32), - } - data = { - "U": np.zeros(dims, dtype=np.float32), - "V": np.zeros(dims, dtype=np.float32), - "P": np.zeros(dims, dtype=np.float32), - } - data["P"][1, 0] = 1.0 - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - fieldset.P.interp_method = "nearest" - xv, yv = np.meshgrid(np.linspace(0.0, 1.0, int(np.sqrt(npart))), np.linspace(0.0, 1.0, int(np.sqrt(npart)))) - pset = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten()) - pset.execute(SampleP, endtime=1, dt=1) - assert np.allclose(pset.p[(pset.lon < 0.5) & (pset.lat > 0.5)], 1.0, rtol=1e-5) - assert np.allclose(pset.p[(pset.lon > 0.5) | (pset.lat < 0.5)], 0.0, rtol=1e-5) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_nearest_neighbor_interpolation3D(): - npart = 81 - dims = (2, 2, 2) - dimensions = { - "lon": np.linspace(0.0, 1.0, dims[2], dtype=np.float32), - "lat": np.linspace(0.0, 1.0, dims[1], dtype=np.float32), - "depth": np.linspace(0.0, 1.0, dims[0], dtype=np.float32), - } - data = { - "U": np.zeros(dims, dtype=np.float32), - "V": np.zeros(dims, dtype=np.float32), - "P": np.zeros(dims, dtype=np.float32), - } - data["P"][1, 1, 0] = 1.0 - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - fieldset.P.interp_method = "nearest" - xv, yv = np.meshgrid(np.linspace(0, 1.0, int(np.sqrt(npart))), np.linspace(0, 1.0, int(np.sqrt(npart)))) - # combine a pset at 0m with pset at 1m, as meshgrid does not do 3D - pset = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten(), depth=np.zeros(npart)) - pset2 = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten(), depth=np.ones(npart)) - pset.add(pset2) - pset.execute(SampleP, endtime=1, dt=1) - assert np.allclose(pset.p[(pset.lon < 0.5) & (pset.lat > 0.5) & (pset.depth > 0.5)], 1.0, rtol=1e-5) - assert np.allclose(pset.p[(pset.lon > 0.5) | (pset.lat < 0.5) & (pset.depth < 0.5)], 0.0, rtol=1e-5) - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("withDepth", [True, False]) -@pytest.mark.parametrize("arrtype", ["ones", "rand"]) -def test_inversedistance_nearland(withDepth, arrtype): - npart = 81 - dims = (4, 4, 6) if withDepth else (4, 6) - dimensions = { - "lon": np.linspace(0.0, 1.0, dims[-1], dtype=np.float32), - "lat": np.linspace(0.0, 1.0, dims[-2], dtype=np.float32), - } - if withDepth: - dimensions["depth"] = np.linspace(0.0, 1.0, dims[0], dtype=np.float32) - P = np.random.rand(dims[0], dims[1], dims[2]) + 2 if arrtype == "rand" else np.ones(dims, dtype=np.float32) - P[1, 1:2, 1:6] = np.nan # setting some values to land (NaN) - else: - P = np.random.rand(dims[0], dims[1]) + 2 if arrtype == "rand" else np.ones(dims, dtype=np.float32) - P[1:2, 1:6] = np.nan # setting some values to land (NaN) - - data = {"U": np.zeros(dims, dtype=np.float32), "V": np.zeros(dims, dtype=np.float32), "P": P} - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - fieldset.P.interp_method = "linear_invdist_land_tracer" - - xv, yv = np.meshgrid(np.linspace(0.1, 0.9, int(np.sqrt(npart))), np.linspace(0.1, 0.9, int(np.sqrt(npart)))) - pset = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten(), depth=np.zeros(npart)) - if withDepth: - pset2 = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten(), depth=np.ones(npart)) - pset.add(pset2) - pset.execute(SampleP, endtime=1, dt=1) - if arrtype == "rand": - assert np.all((pset.p > 2) & (pset.p < 3)) - else: - assert np.allclose(pset.p, 1.0, rtol=1e-5) - - success = False - try: - fieldset.U.interp_method = "linear_invdist_land_tracer" - fieldset._check_complete() - except NotImplementedError: - success = True - assert success - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("boundaryslip", ["freeslip", "partialslip"]) -@pytest.mark.parametrize("withW", [False, True]) -@pytest.mark.parametrize("withT", [False, True]) -def test_partialslip_nearland_zonal(boundaryslip, withW, withT): - npart = 20 - dims = (3, 9, 3) - U = 0.1 * np.ones(dims, dtype=np.float32) - U[:, 0, :] = np.nan - U[:, -1, :] = np.nan - V = np.zeros(dims, dtype=np.float32) - V[:, 0, :] = np.nan - V[:, -1, :] = np.nan - dimensions = { - "lon": np.linspace(-10, 10, dims[2]), - "lat": np.linspace(0.0, 4.0, dims[1], dtype=np.float32), - "depth": np.linspace(-10, 10, dims[0]), - } - if withT: - dimensions["time"] = [0, 2] - U = np.tile(U, (2, 1, 1, 1)) - V = np.tile(V, (2, 1, 1, 1)) - if withW: - W = 0.1 * np.ones(dims, dtype=np.float32) - W[:, 0, :] = np.nan - W[:, -1, :] = np.nan - if withT: - W = np.tile(W, (2, 1, 1, 1)) - data = {"U": U, "V": V, "W": W} - else: - data = {"U": U, "V": V} - fieldset = FieldSet.from_data(data, dimensions, mesh="flat", interp_method=boundaryslip) - - pset = ParticleSet( - fieldset, pclass=Particle, lon=np.zeros(npart), lat=np.linspace(0.1, 3.9, npart), depth=np.zeros(npart) - ) - kernel = AdvectionRK4_3D if withW else AdvectionRK4 - pset.execute(kernel, endtime=2, dt=1) - if boundaryslip == "partialslip": - assert np.allclose([p.lon for p in pset if p.lat >= 0.5 and p.lat <= 3.5], 0.1) - assert np.allclose([pset[0].lon, pset[-1].lon], 0.06) - assert np.allclose([pset[1].lon, pset[-2].lon], 0.08) - if withW: - assert np.allclose([p.depth for p in pset if p.lat >= 0.5 and p.lat <= 3.5], 0.1) - assert np.allclose([pset[0].depth, pset[-1].depth], 0.06) - assert np.allclose([pset[1].depth, pset[-2].depth], 0.08) - else: - assert np.allclose([p.lon for p in pset], 0.1) - if withW: - assert np.allclose([p.depth for p in pset], 0.1) - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("boundaryslip", ["freeslip", "partialslip"]) -@pytest.mark.parametrize("withW", [False, True]) -def test_partialslip_nearland_meridional(boundaryslip, withW): - npart = 20 - dims = (1, 1, 9) - U = np.zeros(dims, dtype=np.float32) - U[:, :, 0] = np.nan - U[:, :, -1] = np.nan - V = 0.1 * np.ones(dims, dtype=np.float32) - V[:, :, 0] = np.nan - V[:, :, -1] = np.nan - dimensions = {"lon": np.linspace(0.0, 4.0, dims[2], dtype=np.float32), "lat": 0, "depth": 0} - if withW: - W = 0.1 * np.ones(dims, dtype=np.float32) - W[:, :, 0] = np.nan - W[:, :, -1] = np.nan - data = {"U": U, "V": V, "W": W} - interp_method = {"U": boundaryslip, "V": boundaryslip, "W": boundaryslip} - else: - data = {"U": U, "V": V} - interp_method = {"U": boundaryslip, "V": boundaryslip} - fieldset = FieldSet.from_data(data, dimensions, mesh="flat", interp_method=interp_method) - - pset = ParticleSet( - fieldset, pclass=Particle, lat=np.zeros(npart), lon=np.linspace(0.1, 3.9, npart), depth=np.zeros(npart) - ) - kernel = AdvectionRK4_3D if withW else AdvectionRK4 - pset.execute(kernel, endtime=2, dt=1) - if boundaryslip == "partialslip": - assert np.allclose([p.lat for p in pset if p.lon >= 0.5 and p.lon <= 3.5], 0.1) - assert np.allclose([pset[0].lat, pset[-1].lat], 0.06) - assert np.allclose([pset[1].lat, pset[-2].lat], 0.08) - if withW: - assert np.allclose([p.depth for p in pset if p.lon >= 0.5 and p.lon <= 3.5], 0.1) - assert np.allclose([pset[0].depth, pset[-1].depth], 0.06) - assert np.allclose([pset[1].depth, pset[-2].depth], 0.08) - else: - assert np.allclose([p.lat for p in pset], 0.1) - if withW: - assert np.allclose([p.depth for p in pset], 0.1) - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("boundaryslip", ["freeslip", "partialslip"]) -def test_partialslip_nearland_vertical(boundaryslip): - npart = 20 - dims = (9, 1, 1) - U = 0.1 * np.ones(dims, dtype=np.float32) - U[0, :, :] = np.nan - U[-1, :, :] = np.nan - V = 0.1 * np.ones(dims, dtype=np.float32) - V[0, :, :] = np.nan - V[-1, :, :] = np.nan - dimensions = {"lon": 0, "lat": 0, "depth": np.linspace(0.0, 4.0, dims[0], dtype=np.float32)} - data = {"U": U, "V": V} - fieldset = FieldSet.from_data(data, dimensions, mesh="flat", interp_method={"U": boundaryslip, "V": boundaryslip}) - - pset = ParticleSet( - fieldset, pclass=Particle, lon=np.zeros(npart), lat=np.zeros(npart), depth=np.linspace(0.1, 3.9, npart) - ) - pset.execute(AdvectionRK4, endtime=2, dt=1) - if boundaryslip == "partialslip": - assert np.allclose([p.lon for p in pset if p.depth >= 0.5 and p.depth <= 3.5], 0.1) - assert np.allclose([p.lat for p in pset if p.depth >= 0.5 and p.depth <= 3.5], 0.1) - assert np.allclose([pset[0].lon, pset[-1].lon, pset[0].lat, pset[-1].lat], 0.06) - assert np.allclose([pset[1].lon, pset[-2].lon, pset[1].lat, pset[-2].lat], 0.08) - else: - assert np.allclose([p.lon for p in pset], 0.1) - assert np.allclose([p.lat for p in pset], 0.1) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_sample_particle(): - """Sample the fieldset using an array of particles.""" - npart = 120 - lon = np.linspace(-180, 180, 200, dtype=np.float32) - lat = np.linspace(-90, 90, 100, dtype=np.float32) - V, U = np.meshgrid(lon, lat) - data = {"U": U, "V": V} - dimensions = {"lon": lon, "lat": lat} - - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - - 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(SampleUV), endtime=1.0, dt=1.0) - assert np.allclose(pset.v, lon, rtol=1e-6) - - pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0) - pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0) - assert np.allclose(pset.u, lat, rtol=1e-6) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_sample_geographic(fieldset_geometric): - """Sample a fieldset with conversion to geographic units (degrees).""" - 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(SampleUV), endtime=1.0, dt=1.0) - assert np.allclose(pset.v, lon, rtol=1e-6) - - pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0) - pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0) - assert np.allclose(pset.u, lat, rtol=1e-6) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_fieldset_sample_geographic_polar(fieldset_geometric_polar): - """Sample a fieldset with conversion to geographic units and a pole correction.""" - npart = 120 - fieldset = fieldset_geometric_polar - 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(SampleUV), endtime=1.0, dt=1.0) - assert np.allclose(pset.v, lon, rtol=1e-6) - - pset = ParticleSet(fieldset, pclass=pclass(), lat=lat, lon=np.zeros(npart) - 45.0) - pset.execute(pset.Kernel(SampleUV), endtime=1.0, dt=1.0) - assert np.allclose(pset.u, lat, rtol=1e-2) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_meridionalflow_spherical(): - """Create uniform NORTHWARD flow on spherical earth and advect particles. - - As flow is so simple, it can be directly compared to analytical solution. - """ - xdim = 100 - ydim = 200 - - maxvel = 1.0 - dimensions = { - "lon": np.linspace(-180, 180, xdim, dtype=np.float32), - "lat": np.linspace(-90, 90, ydim, dtype=np.float32), - } - data = {"U": np.zeros([ydim, xdim]), "V": maxvel * np.ones([ydim, xdim])} - - fieldset = FieldSet.from_data(data, dimensions, mesh="spherical") - - lonstart = [0, 45] - latstart = [0, 45] - runtime = timedelta(hours=24) - pset = ParticleSet(fieldset, pclass=Particle, lon=lonstart, lat=latstart) - pset.execute(pset.Kernel(AdvectionRK4), runtime=runtime, dt=timedelta(hours=1)) - - assert pset.lat[0] - (latstart[0] + runtime.total_seconds() * maxvel / 1852 / 60) < 1e-4 - assert pset.lon[0] - lonstart[0] < 1e-4 - assert pset.lat[1] - (latstart[1] + runtime.total_seconds() * maxvel / 1852 / 60) < 1e-4 - assert pset.lon[1] - lonstart[1] < 1e-4 - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_zonalflow_spherical(): - """Create uniform EASTWARD flow on spherical earth and advect particles. - - As flow is so simple, it can be directly compared to analytical solution - Note that in this case the cosine conversion is needed - """ - xdim, ydim = 100, 200 - - maxvel = 1.0 - p_fld = 10 - dimensions = { - "lon": np.linspace(-180, 180, xdim, dtype=np.float32), - "lat": np.linspace(-90, 90, ydim, dtype=np.float32), - } - data = {"U": maxvel * np.ones([ydim, xdim]), "V": np.zeros([ydim, xdim]), "P": p_fld * np.ones([ydim, xdim])} - - fieldset = FieldSet.from_data(data, dimensions, mesh="spherical") - - lonstart = [0, 45] - latstart = [0, 45] - runtime = timedelta(hours=24) - pset = ParticleSet(fieldset, pclass=pclass(), lon=lonstart, lat=latstart) - pset.execute(pset.Kernel(AdvectionRK4) + SampleP, runtime=runtime, dt=timedelta(hours=1)) - - assert pset.lat[0] - latstart[0] < 1e-4 - assert ( - pset.lon[0] - (lonstart[0] + runtime.total_seconds() * maxvel / 1852 / 60 / cos(latstart[0] * pi / 180)) < 1e-4 - ) - assert abs(pset.p[0] - p_fld) < 1e-4 - assert pset.lat[1] - latstart[1] < 1e-4 - assert ( - pset.lon[1] - (lonstart[1] + runtime.total_seconds() * maxvel / 1852 / 60 / cos(latstart[1] * pi / 180)) < 1e-4 - ) - assert abs(pset.p[1] - p_fld) < 1e-4 - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -def test_random_field(): - """Sampling test that tests for overshoots by sampling a field of random numbers between 0 and 1.""" - xdim, ydim = 20, 20 - npart = 100 - - np.random.seed(123456) - dimensions = { - "lon": np.linspace(0.0, 1.0, xdim, dtype=np.float32), - "lat": np.linspace(0.0, 1.0, ydim, dtype=np.float32), - } - data = { - "U": np.zeros((ydim, xdim), dtype=np.float32), - "V": np.zeros((ydim, xdim), dtype=np.float32), - "P": np.random.uniform(0, 1.0, size=(ydim, xdim)), - "start": np.ones((ydim, xdim), dtype=np.float32), - } - - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - pset = ParticleSet.from_field(fieldset, size=npart, pclass=pclass(), start_field=fieldset.start) - pset.execute(SampleP, endtime=1.0, dt=1.0) - sampled = pset.p - assert (sampled >= 0.0).all() - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("allow_time_extrapolation", [True, False]) -def test_sampling_out_of_bounds_time(allow_time_extrapolation): - xdim, ydim, tdim = 10, 10, 10 - - dimensions = { - "lon": np.linspace(0.0, 1.0, xdim, dtype=np.float32), - "lat": np.linspace(0.0, 1.0, ydim, dtype=np.float32), - "time": np.linspace(0.0, 1.0, tdim, dtype=np.float64), - } - data = { - "U": np.zeros((tdim, ydim, xdim), dtype=np.float32), - "V": np.zeros((tdim, ydim, xdim), dtype=np.float32), - "P": np.transpose(np.ones((xdim, ydim, 1), dtype=np.float32) * dimensions["time"]), - } - - fieldset = FieldSet.from_data(data, dimensions, mesh="flat", allow_time_extrapolation=allow_time_extrapolation) - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=-1.0) - if allow_time_extrapolation: - pset.execute(SampleP, endtime=-0.9, dt=0.1) - assert np.allclose(pset.p, 0.0, rtol=1e-5) - else: - with pytest.raises(RuntimeError): - pset.execute(SampleP, endtime=-0.9, dt=0.1) - - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=0) - pset.execute(SampleP, runtime=0.1, dt=0.1) - assert np.allclose(pset.p, 0.0, rtol=1e-5) - - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=0.5) - pset.execute(SampleP, runtime=0.1, dt=0.1) - assert np.allclose(pset.p, 0.5, rtol=1e-5) - - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=1.0) - pset.execute(SampleP, runtime=0.1, dt=0.1) - assert np.allclose(pset.p, 1.0, rtol=1e-5) - - pset = ParticleSet(fieldset, pclass=pclass(), lon=[0.5], lat=[0.5], time=2.0) - if allow_time_extrapolation: - pset.execute(SampleP, runtime=0.1, dt=0.1) - assert np.allclose(pset.p, 1.0, rtol=1e-5) - else: - with pytest.raises(RuntimeError): - pset.execute(SampleP, runtime=0.1, dt=0.1) - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="When refactoring fieldfilebuffer croco support was dropped. This will be fixed in v4.") -def test_sampling_3DCROCO(): - data_path = os.path.join(os.path.dirname(__file__), "test_data/") - fieldset = FieldSet.from_modulefile(data_path + "fieldset_CROCO3D.py") - - SampleP = Particle.add_variable("p", initial=0.0) - - def SampleU(particle, fieldset, time): # pragma: no cover - particle.p = fieldset.U[time, particle.depth, particle.lat, particle.lon, particle] - - pset = ParticleSet(fieldset, pclass=SampleP, lon=120e3, lat=50e3, depth=-0.4) - pset.execute(SampleU, endtime=1, dt=1) - assert np.isclose(pset.p, 1.0) - - -@pytest.mark.v4alpha -@pytest.mark.xfail( - reason="Now timestamps has been removed, and filebuffer expects different files. This needs to be rewritten/removed." -) -@pytest.mark.parametrize("npart", [1, 10]) -def test_sampling_multigrids_non_vectorfield_from_file(npart, tmpdir): - xdim, ydim = 100, 200 - filepath = tmpdir.join("test_subsets") - U = Field( - "U", - np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - V = Field( - "V", - np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - B = Field( - "B", - np.ones((3 * ydim, 4 * xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, 4 * xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, 3 * ydim, dtype=np.float32), - ) - fieldset = FieldSet(U, V) - fieldset.add_field(B, "B") - fieldset.write(filepath) - fieldset = None - - ufiles = [filepath + "U.nc"] * 4 - vfiles = [filepath + "V.nc"] * 4 - bfiles = [filepath + "B.nc"] * 4 - timestamps = np.arange(0, 4, 1) * 86400.0 - timestamps = np.expand_dims(timestamps, 1) - files = {"U": ufiles, "V": vfiles, "B": bfiles} - variables = {"U": "vozocrtx", "V": "vomecrty", "B": "B"} - dimensions = {"lon": "nav_lon", "lat": "nav_lat"} - fieldset = FieldSet.from_netcdf(files, variables, dimensions, timestamps=timestamps, allow_time_extrapolation=True) - - fieldset.add_constant("sample_depth", 2.5) - assert fieldset.U.grid is fieldset.V.grid - assert fieldset.U.grid is not fieldset.B.grid - - TestParticle = Particle.add_variable("sample_var", initial=0.0) - - pset = ParticleSet.from_line(fieldset, pclass=TestParticle, start=[0.3, 0.3], finish=[0.7, 0.7], size=npart) - - def test_sample(particle, fieldset, time): # pragma: no cover - particle.sample_var += fieldset.B[time, fieldset.sample_depth, particle.lat, particle.lon] - - kernels = pset.Kernel(AdvectionRK4) + pset.Kernel(test_sample) - pset.execute(kernels, runtime=10, dt=1) - assert np.allclose(pset.sample_var, 10.0) - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("npart", [1, 10]) -def test_sampling_multigrids_non_vectorfield(npart): - xdim, ydim = 100, 200 - U = Field( - "U", - np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - V = Field( - "V", - np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - B = Field( - "B", - np.ones((3 * ydim, 4 * xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, 4 * xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, 3 * ydim, dtype=np.float32), - ) - fieldset = FieldSet(U, V) - fieldset.add_field(B, "B") - fieldset.add_constant("sample_depth", 2.5) - assert fieldset.U.grid is fieldset.V.grid - assert fieldset.U.grid is not fieldset.B.grid - - TestParticle = Particle.add_variable("sample_var", initial=0.0) - - pset = ParticleSet.from_line(fieldset, pclass=TestParticle, start=[0.3, 0.3], finish=[0.7, 0.7], size=npart) - - def test_sample(particle, fieldset, time): # pragma: no cover - particle.sample_var += fieldset.B[time, fieldset.sample_depth, particle.lat, particle.lon] - - kernels = pset.Kernel(AdvectionRK4) + pset.Kernel(test_sample) - pset.execute(kernels, runtime=10, dt=1) - assert np.allclose(pset.sample_var, 10.0) - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.parametrize("ugridfactor", [1, 10]) -def test_sampling_multiple_grid_sizes(ugridfactor): - xdim, ydim = 10, 20 - U = Field( - "U", - np.zeros((ydim * ugridfactor, xdim * ugridfactor), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim * ugridfactor, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim * ugridfactor, dtype=np.float32), - ) - V = Field( - "V", - np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - fieldset = FieldSet(U, V) - pset = ParticleSet(fieldset, pclass=Particle, lon=[0.8], lat=[0.9]) - - if ugridfactor > 1: - assert fieldset.U.grid is not fieldset.V.grid - else: - assert fieldset.U.grid is fieldset.V.grid - pset.execute(AdvectionRK4, runtime=10, dt=1) - assert np.isclose(pset.lon[0], 0.8) - assert np.all((0 <= pset.xi) & (pset.xi < xdim * ugridfactor)) - - -@pytest.mark.v4future -@pytest.mark.xfail(reason="GH1946") -def test_multiple_grid_addlater_error(): - xdim, ydim = 10, 20 - U = Field( - "U", - np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - V = Field( - "V", - np.zeros((ydim, xdim), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim, dtype=np.float32), - ) - fieldset = FieldSet(U, V) - - pset = ParticleSet(fieldset, pclass=Particle, lon=[0.8], lat=[0.9]) # noqa ; to trigger fieldset._check_complete - - P = Field( - "P", - np.zeros((ydim * 10, xdim * 10), dtype=np.float32), - lon=np.linspace(0.0, 1.0, xdim * 10, dtype=np.float32), - lat=np.linspace(0.0, 1.0, ydim * 10, dtype=np.float32), - ) - - fail = False - try: - fieldset.add_field(P) - except RuntimeError: - fail = True - assert fail - - -def test_fieldset_sampling_updating_order(tmp_zarrfile): - def calc_p(t, y, x): - return 10 * t + x + 0.2 * y - - dims = [2, 4, 5] - dimensions = { - "lon": np.linspace(0.0, 1.0, dims[2], dtype=np.float32), - "lat": np.linspace(0.0, 1.0, dims[1], dtype=np.float32), - "time": np.arange(dims[0], dtype=np.float32), - } - - p = np.zeros(dims, dtype=np.float32) - for i, x in enumerate(dimensions["lon"]): - for j, y in enumerate(dimensions["lat"]): - for n, t in enumerate(dimensions["time"]): - p[n, j, i] = calc_p(t, y, x) - - data = { - "U": 0.5 * np.ones(dims, dtype=np.float32), - "V": np.zeros(dims, dtype=np.float32), - "P": p, - } - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - - xv, yv = np.meshgrid(np.arange(0, 1, 0.5), np.arange(0, 1, 0.5)) - pset = ParticleSet(fieldset, pclass=pclass(), lon=xv.flatten(), lat=yv.flatten()) - - def SampleP(particle, fieldset, time): # pragma: no cover - particle.p = fieldset.P[time, particle.depth, particle.lat, particle.lon] - - kernels = [AdvectionRK4, SampleP] - - pfile = pset.ParticleFile(tmp_zarrfile, outputdt=1) - pset.execute(kernels, endtime=1, dt=1, output_file=pfile) - - ds = xr.open_zarr(tmp_zarrfile) - for t in range(len(ds["obs"])): - for i in range(len(ds["trajectory"])): - assert np.isclose( - ds["p"].values[i, t], - calc_p(float(ds["time"].values[i, t]) / 1e9, ds["lat"].values[i, t], ds["lon"].values[i, t]), - ) diff --git a/tests-v3/test_interpolation.py b/tests-v3/test_interpolation.py deleted file mode 100644 index 67eab64204..0000000000 --- a/tests-v3/test_interpolation.py +++ /dev/null @@ -1,64 +0,0 @@ -import pytest - -import parcels._interpolation as interpolation -from tests.utils import create_fieldset_zeros_3d - - -@pytest.fixture -def tmp_interpolator_registry(): - """Resets the interpolator registry after the test. Vital when testing manipulating the registry.""" - old_2d = interpolation._interpolator_registry_2d.copy() - old_3d = interpolation._interpolator_registry_3d.copy() - yield - interpolation._interpolator_registry_2d = old_2d - interpolation._interpolator_registry_3d = old_3d - - -@pytest.mark.usefixtures("tmp_interpolator_registry") -def test_interpolation_registry(): - @interpolation.register_3d_interpolator("test") - @interpolation.register_2d_interpolator("test") - def some_function(): - return "test" - - assert "test" in interpolation.get_2d_interpolator_registry() - assert "test" in interpolation.get_3d_interpolator_registry() - - f = interpolation.get_2d_interpolator_registry()["test"] - g = interpolation.get_3d_interpolator_registry()["test"] - assert f() == g() == "test" - - -@pytest.mark.v4remove -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.usefixtures("tmp_interpolator_registry") -def test_interpolator_override(): - fieldset = create_fieldset_zeros_3d() - - @interpolation.register_3d_interpolator("linear") - def test_interpolator(ctx: interpolation.InterpolationContext3D): - raise NotImplementedError - - with pytest.raises(NotImplementedError): - fieldset.U[0, 0.5, 0.5, 0.5] - - -@pytest.mark.v4remove -@pytest.mark.xfail(reason="GH1946") -@pytest.mark.usefixtures("tmp_interpolator_registry") -def test_full_depth_provided_to_interpolators(): - """The full depth needs to be provided to the interpolation schemes as some interpolators - need to know whether they are at the surface or bottom of the water column. - - https://github.com/OceanParcels/Parcels/pull/1816#discussion_r1908840408 - """ - xdim, ydim, zdim = 10, 11, 12 - fieldset = create_fieldset_zeros_3d(xdim=xdim, ydim=ydim, zdim=zdim) - - @interpolation.register_3d_interpolator("linear") - def test_interpolator2(ctx: interpolation.InterpolationContext3D): - assert ctx.data.shape[1] == zdim - # The array z dimension is the same as the fieldset z dimension - return 0 - - fieldset.U[0.5, 0.5, 0.5, 0.5] diff --git a/tests-v3/test_kernel_execution.py b/tests-v3/test_kernel_execution.py deleted file mode 100644 index 6a61f29860..0000000000 --- a/tests-v3/test_kernel_execution.py +++ /dev/null @@ -1,84 +0,0 @@ -import numpy as np -import pytest - -from parcels import ( - FieldSet, - Particle, - ParticleSet, -) -from tests.utils import create_fieldset_unit_mesh - - -@pytest.fixture -def fieldset_unit_mesh(): - return create_fieldset_unit_mesh() - - -@pytest.mark.parametrize("kernel_type", ["update_lon", "update_dlon"]) -def test_execution_order(kernel_type): - fieldset = FieldSet.from_data( - {"U": [[0, 1], [2, 3]], "V": np.ones((2, 2))}, {"lon": [0, 2], "lat": [0, 2]}, mesh="flat" - ) - - def MoveLon_Update_Lon(particle, fieldset, time): # pragma: no cover - particle.lon += 0.2 - - def MoveLon_Update_dlon(particle, fieldset, time): # pragma: no cover - particle.dlon += 0.2 - - def SampleP(particle, fieldset, time): # pragma: no cover - particle.p = fieldset.U[time, particle.depth, particle.lat, particle.lon] - - SampleParticle = Particle.add_variable("p", dtype=np.float32, initial=0.0) - - MoveLon = MoveLon_Update_dlon if kernel_type == "update_dlon" else MoveLon_Update_Lon - - kernels = [MoveLon, SampleP] - lons = [] - ps = [] - for dir in [1, -1]: - pset = ParticleSet(fieldset, pclass=SampleParticle, lon=0, lat=0) - pset.execute(kernels[::dir], endtime=1, dt=1) - lons.append(pset.lon) - ps.append(pset.p) - - if kernel_type == "update_dlon": - assert np.isclose(lons[0], lons[1]) - assert np.isclose(ps[0], ps[1]) - assert np.allclose(lons[0], 0) - else: - assert np.isclose(ps[0] - ps[1], 0.1) - assert np.allclose(lons[0], 0.2) - - -def test_multi_kernel_duplicate_varnames(fieldset_unit_mesh): - # Testing for merging of two Kernels with the same variable declared - # Should throw a warning, but go ahead regardless - def Kernel1(particle, fieldset, time): # pragma: no cover - add_lon = 0.1 - particle.dlon += add_lon - - def Kernel2(particle, fieldset, time): # pragma: no cover - add_lon = -0.3 - particle.dlon += add_lon - - pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=[0.5], lat=[0.5]) - pset.execute([Kernel1, Kernel2], endtime=2.0, dt=1.0) - assert np.allclose(pset.lon, 0.3, rtol=1e-5) - - -def test_update_kernel_in_script(fieldset_unit_mesh): - # Testing what happens when kernels are updated during runtime of a script - # Should throw a warning, but go ahead regardless - def MoveEast(particle, fieldset, time): # pragma: no cover - add_lon = 0.1 - particle.dlon += add_lon - - def MoveWest(particle, fieldset, time): # pragma: no cover - add_lon = -0.3 - particle.dlon += add_lon - - pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=[0.5], lat=[0.5]) - pset.execute(pset.Kernel(MoveEast), endtime=1.0, dt=1.0) - pset.execute(pset.Kernel(MoveWest), endtime=3.0, dt=1.0) - assert np.allclose(pset.lon, 0.3, rtol=1e-5) # should be 0.5 + 0.1 - 0.3 = 0.3 diff --git a/tests-v3/test_kernel_language.py b/tests-v3/test_kernel_language.py deleted file mode 100644 index eaae53f211..0000000000 --- a/tests-v3/test_kernel_language.py +++ /dev/null @@ -1,301 +0,0 @@ -import random -from contextlib import nullcontext as does_not_raise - -import numpy as np -import pytest - -from parcels import ( - Field, - Kernel, - Particle, - ParticleSet, - Variable, -) -from tests.common_kernels import DoNothing -from tests.utils import create_fieldset_unit_mesh - - -def expr_kernel(name, pset, expr): - pycode = (f"def {name}(particle, fieldset, time):\n" - f" particle.p = {expr}") # fmt: skip - return Kernel(kernels=None, fieldset=pset.fieldset, pclass=pset._pclass, funccode=pycode, funcname=name) - - -@pytest.fixture -def fieldset_unit_mesh(): - return create_fieldset_unit_mesh() - - -@pytest.mark.parametrize( - "name, expr, result", - [ - ("Add", "2 + 5", 7), - ("Sub", "6 - 2", 4), - ("Mul", "3 * 5", 15), - ("Div", "24 / 4", 6), - ], -) -def test_expression_int(name, expr, result): - """Test basic arithmetic expressions.""" - npart = 10 - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet( - create_fieldset_unit_mesh(mesh="spherical"), - pclass=TestParticle, - lon=np.linspace(0.0, 1.0, npart), - lat=np.zeros(npart) + 0.5, - ) - pset.execute(expr_kernel(f"Test{name}", pset, expr), endtime=1.0, dt=1.0) - assert np.all([p.p == result for p in pset]) - - -@pytest.mark.parametrize( - "name, expr, result", - [ - ("Add", "2. + 5.", 7), - ("Sub", "6. - 2.", 4), - ("Mul", "3. * 5.", 15), - ("Div", "24. / 4.", 6), - ("Pow", "2 ** 3", 8), - ], -) -def test_expression_float(name, expr, result): - """Test basic arithmetic expressions.""" - npart = 10 - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet( - create_fieldset_unit_mesh(mesh="spherical"), - pclass=TestParticle, - lon=np.linspace(0.0, 1.0, npart), - lat=np.zeros(npart) + 0.5, - ) - pset.execute(expr_kernel(f"Test{name}", pset, expr), endtime=1.0, dt=1.0) - assert np.all([p.p == result for p in pset]) - - -@pytest.mark.parametrize( - "name, expr, result", - [ - ("True", "True", True), - ("False", "False", False), - ("And", "True and False", False), - ("Or", "True or False", True), - ("Equal", "5 == 5", True), - ("NotEqual", "3 != 4", True), - ("Lesser", "5 < 3", False), - ("LesserEq", "3 <= 5", True), - ("Greater", "4 > 2", True), - ("GreaterEq", "2 >= 4", False), - ("CheckNaN", "math.nan != math.nan", True), - ], -) -def test_expression_bool(name, expr, result): - """Test basic arithmetic expressions.""" - npart = 10 - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet( - create_fieldset_unit_mesh(mesh="spherical"), - pclass=TestParticle, - lon=np.linspace(0.0, 1.0, npart), - lat=np.zeros(npart) + 0.5, - ) - pset.execute(expr_kernel(f"Test{name}", pset, expr), endtime=1.0, dt=1.0) - assert np.all(result == pset.p) - - -def test_while_if_break(): - """Test while, if and break commands.""" - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=[0], lat=[0]) - - def kernel(particle, fieldset, time): # pragma: no cover - while particle.p < 30: - if particle.p > 9: - break - particle.p += 1 - if particle.p > 5: - particle.p *= 2.0 - - pset.execute(kernel, endtime=1.0, dt=1.0) - assert np.allclose(pset.p, 20.0, rtol=1e-12) - - -def test_nested_if(): - """Test nested if commands.""" - TestParticle = Particle.add_variables( - [Variable("p0", dtype=np.int32, initial=0), Variable("p1", dtype=np.int32, initial=1)] - ) - pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=0, lat=0) - - def kernel(particle, fieldset, time): # pragma: no cover - if particle.p1 >= particle.p0: - var = particle.p0 - if var + 1 < particle.p1: - particle.p1 = -1 - - pset.execute(kernel, endtime=10, dt=1.0) - assert np.allclose([pset.p0[0], pset.p1[0]], [0, 1]) - - -def test_pass(): - """Test pass commands.""" - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=0, lat=0) - - def kernel(particle, fieldset, time): # pragma: no cover - particle.p = -1 - pass - - pset.execute(kernel, endtime=10, dt=1.0) - assert np.allclose(pset[0].p, -1) - - -def test_dt_as_variable_in_kernel(): - pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=Particle, lon=0, lat=0) - - def kernel(particle, fieldset, time): # pragma: no cover - dt = 1.0 # noqa - - pset.execute(kernel, endtime=10, dt=1.0) - - -def test_varname_as_fieldname(): - """Tests for error thrown if variable has same name as Field.""" - fset = create_fieldset_unit_mesh(mesh="spherical") - fset.add_field(Field("speed", 10, lon=0, lat=0)) - fset.add_constant("vertical_speed", 0.1) - particle = Particle.add_variable("speed") - pset = ParticleSet(fset, pclass=particle, lon=0, lat=0) - - def kernel_particlename(particle, fieldset, time): # pragma: no cover - particle.speed = fieldset.speed[particle] - - pset.execute(kernel_particlename, endtime=1, dt=1.0) - assert pset[0].speed == 10 - - def kernel_varname(particle, fieldset, time): # pragma: no cover - vertical_speed = fieldset.vertical_speed # noqa - - pset.execute(kernel_varname, endtime=1, dt=1.0) - - -def test_if_withfield(fieldset_unit_mesh): - """Test combination of if and Field sampling commands.""" - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet(fieldset_unit_mesh, pclass=TestParticle, lon=[0], lat=[0]) - - def kernel(particle, fieldset, time): # pragma: no cover - u, v = fieldset.UV[time, 0, 0, 1.0] - particle.p = 0 - if fieldset.U[time, 0, 0, 1.0] == u: - particle.p += 1 - if fieldset.U[time, 0, 0, 1.0] == fieldset.U[time, 0, 0, 1.0]: - particle.p += 1 - if True: - particle.p += 1 - if fieldset.U[time, 0, 0, 1.0] == u and 1 == 1: - particle.p += 1 - if ( - fieldset.U[time, 0, 0, 1.0] == fieldset.U[time, 0, 0, 1.0] - and fieldset.U[time, 0, 0, 1.0] == fieldset.U[time, 0, 0, 1.0] - ): - particle.p += 1 - if fieldset.U[time, 0, 0, 1.0] == u: - particle.p += 1 - else: - particle.p += 1000 - if fieldset.U[time, 0, 0, 1.0] == 3: - particle.p += 1000 - else: - particle.p += 1 - - pset.execute(kernel, endtime=1.0, dt=1.0) - assert np.allclose(pset.p, 7.0, rtol=1e-12) - return - - -def test_print(fieldset_unit_mesh, capfd): - """Test print statements.""" - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet(fieldset_unit_mesh, pclass=TestParticle, lon=[0.5], lat=[0.5]) - - def kernel(particle, fieldset, time): # pragma: no cover - particle.p = 1e-3 - tmp = 5 - print(f"{particle.trajectory} {particle.p:f} {tmp:f}") - - pset.execute(kernel, endtime=1.0, dt=1.0, verbose_progress=False) - out, err = capfd.readouterr() - lst = out.split(" ") - tol = 1e-8 - assert ( - abs(float(lst[0]) - pset.trajectory[0]) < tol - and abs(float(lst[1]) - pset.p[0]) < tol - and abs(float(lst[2]) - 5) < tol - ) - - def kernel2(particle, fieldset, time): # pragma: no cover - tmp = 3 - print(f"{tmp:f}") - - pset.execute(kernel2, endtime=2.0, dt=1.0, verbose_progress=False) - out, err = capfd.readouterr() - lst = out.split(" ") - assert abs(float(lst[0]) - 3) < tol - - -def test_fieldset_access(fieldset_unit_mesh): - pset = ParticleSet(fieldset_unit_mesh, pclass=Particle, lon=0, lat=0) - - def kernel(particle, fieldset, time): # pragma: no cover - particle.lon = fieldset.U.grid.lon[2] - - pset.execute(kernel, endtime=1, dt=1.0) - assert pset.lon[0] == fieldset_unit_mesh.U.grid.lon[2] - - -@pytest.mark.parametrize("concat", [False, True]) -def test_random_kernel_concat(fieldset_unit_mesh, concat): - TestParticle = Particle.add_variable("p", dtype=np.float32, initial=0) - pset = ParticleSet(fieldset_unit_mesh, pclass=TestParticle, lon=0, lat=0) - - def RandomKernel(particle, fieldset, time): # pragma: no cover - particle.p += random.uniform(0, 1) - - def AddOne(particle, fieldset, time): # pragma: no cover - particle.p += 1.0 - - kernels = [RandomKernel, AddOne] if concat else RandomKernel - pset.execute(kernels, runtime=1) - assert pset.p > 1 if concat else pset.p < 1 - - -def test_dt_modif_by_kernel(): - TestParticle = Particle.add_variable("age", dtype=np.float32, initial=0) - pset = ParticleSet(create_fieldset_unit_mesh(mesh="spherical"), pclass=TestParticle, lon=[0.5], lat=[0]) - - def modif_dt(particle, fieldset, time): # pragma: no cover - particle.age += particle.dt - particle.dt = 2 - - endtime = 4 - pset.execute(modif_dt, endtime=endtime + 1, dt=1.0) - assert np.isclose(pset.time[0], endtime) - - -@pytest.mark.parametrize( - ("dt", "expectation"), [(1e-2, does_not_raise()), (1e-5, does_not_raise()), (1e-6, pytest.raises(ValueError))] -) -def test_small_dt(dt, expectation): - npart = 10 - pset = ParticleSet( - create_fieldset_unit_mesh(mesh="spherical"), - pclass=Particle, - lon=np.zeros(npart), - lat=np.zeros(npart), - time=np.arange(0, npart) * dt * 10, - ) - - with expectation: - pset.execute(DoNothing, dt=dt, runtime=dt * 101) - assert np.allclose([p.time for p in pset], dt * 100) diff --git a/tests-v3/test_particles.py b/tests-v3/test_particles.py deleted file mode 100644 index 172d2a4291..0000000000 --- a/tests-v3/test_particles.py +++ /dev/null @@ -1,88 +0,0 @@ -from operator import attrgetter - -import numpy as np -import pytest - -from parcels import ( - AdvectionRK4, - Particle, - ParticleSet, - Variable, -) -from tests.utils import create_fieldset_zeros_unit_mesh - - -@pytest.fixture -def fieldset(): - return create_fieldset_zeros_unit_mesh() - - -def test_print(fieldset): - TestParticle = Particle.add_variable("p", to_write=True) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=[0, 1], lat=[0, 1]) - print(pset) - - -def test_variable_init(fieldset): - """Test that checks correct initialisation of custom variables.""" - npart = 10 - extra_vars = [ - Variable("p_float", dtype=np.float32, initial=10.0), - Variable("p_double", dtype=np.float64, initial=11.0), - ] - TestParticle = Particle.add_variables(extra_vars) - TestParticle = TestParticle.add_variable("p_int", np.int32, initial=12.0) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart)) - - def addOne(particle, fieldset, time): # pragma: no cover - particle.p_float += 1.0 - particle.p_double += 1.0 - particle.p_int += 1 - - pset.execute(pset.Kernel(AdvectionRK4) + addOne, runtime=1.0, dt=1.0) - assert np.allclose([p.p_float for p in pset], 11.0, rtol=1e-12) - assert np.allclose([p.p_double for p in pset], 12.0, rtol=1e-12) - assert np.allclose([p.p_int for p in pset], 13, rtol=1e-12) - - -@pytest.mark.parametrize("type", ["np.int8", "mp.float", "np.int16"]) -def test_variable_unsupported_dtypes(fieldset, type): - """Test that checks errors thrown for unsupported dtypes.""" - TestParticle = Particle.add_variable("p", dtype=type, initial=10.0) - with pytest.raises((RuntimeError, TypeError)): - ParticleSet(fieldset, pclass=TestParticle, lon=[0], lat=[0]) - - -def test_variable_special_names(fieldset): - """Test that checks errors thrown for special names.""" - for vars in ["z", "lon"]: - TestParticle = Particle.add_variable(vars, dtype=np.float32, initial=10.0) - with pytest.raises(AttributeError): - ParticleSet(fieldset, pclass=TestParticle, lon=[0], lat=[0]) - - -@pytest.mark.parametrize("coord_type", [np.float32, np.float64]) -def test_variable_init_relative(fieldset, coord_type): - """Test that checks relative initialisation of custom variables.""" - npart = 10 - lonlat_type = np.float64 if coord_type == "double" else np.float32 - - TestParticle = Particle.add_variables( - [ - Variable("p_base", dtype=lonlat_type, initial=10.0), - Variable("p_relative", dtype=lonlat_type, initial=attrgetter("p_base")), - Variable("p_lon", dtype=lonlat_type, initial=attrgetter("lon")), - Variable("p_lat", dtype=lonlat_type, initial=attrgetter("lat")), - ] - ) - - lon = np.linspace(0, 1, npart, dtype=lonlat_type) - lat = np.linspace(1, 0, npart, dtype=lonlat_type) - pset = ParticleSet(fieldset, pclass=TestParticle, lon=lon, lat=lat, lonlatdepth_dtype=coord_type) - # Adjust base variable to test for aliasing effects - for p in pset: - p.p_base += 3.0 - assert np.allclose([p.p_base for p in pset], 13.0, rtol=1e-12) - assert np.allclose([p.p_relative for p in pset], 10.0, rtol=1e-12) - assert np.allclose([p.p_lon for p in pset], lon, rtol=1e-12) - assert np.allclose([p.p_lat for p in pset], lat, rtol=1e-12) diff --git a/tests-v3/test_particlesets.py b/tests-v3/test_particlesets.py deleted file mode 100644 index 73bbdb204a..0000000000 --- a/tests-v3/test_particlesets.py +++ /dev/null @@ -1,258 +0,0 @@ -import numpy as np -import pytest - -from parcels import ( - CurvilinearZGrid, - Field, - FieldSet, - Particle, - ParticleSet, - Variable, -) -from tests.utils import create_fieldset_zeros_simple - - -@pytest.fixture -def fieldset(): - return create_fieldset_zeros_simple() - - -@pytest.fixture -def pset(fieldset): - npart = 10 - pset = ParticleSet(fieldset, pclass=Particle, lon=np.linspace(0, 1, npart), lat=np.zeros(npart)) - return pset - - -def test_pset_create_list_with_customvariable(fieldset): - npart = 100 - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) - - MyParticle = Particle.add_variable("v") - - v_vals = np.arange(npart) - pset = ParticleSet(fieldset, lon=lon, lat=lat, v=v_vals, pclass=MyParticle) - assert np.allclose([p.lon for p in pset], lon, rtol=1e-12) - assert np.allclose([p.lat for p in pset], lat, rtol=1e-12) - assert np.allclose([p.v for p in pset], v_vals, rtol=1e-12) - - -@pytest.mark.parametrize("restart", [True, False]) -def test_pset_create_fromparticlefile(fieldset, restart, tmp_zarrfile): - lon = np.linspace(0, 1, 10, dtype=np.float32) - lat = np.linspace(1, 0, 10, dtype=np.float32) - - TestParticle = Particle.add_variable("p", np.float32, initial=0.33) - TestParticle = TestParticle.add_variable("p2", np.float32, initial=1, to_write=False) - TestParticle = TestParticle.add_variable("p3", np.float64, to_write="once") - - pset = ParticleSet(fieldset, lon=lon, lat=lat, depth=[4] * len(lon), pclass=TestParticle, p3=np.arange(len(lon))) - pfile = pset.ParticleFile(tmp_zarrfile, outputdt=1) - - def Kernel(particle, fieldset, time): # pragma: no cover - particle.p = 2.0 - if particle.lon == 1.0: - particle.delete() - - pset.execute(Kernel, runtime=2, dt=1, output_file=pfile) - - pset_new = ParticleSet.from_particlefile( - fieldset, pclass=TestParticle, filename=tmp_zarrfile, restart=restart, repeatdt=1 - ) - - for var in ["lon", "lat", "depth", "time", "p", "p2", "p3"]: - assert np.allclose([getattr(p, var) for p in pset], [getattr(p, var) for p in pset_new]) - - if restart: - assert np.allclose([p.trajectory for p in pset], [p.trajectory for p in pset_new]) - pset_new.execute(Kernel, runtime=2, dt=1) - assert len(pset_new) == 3 * len(pset) - assert pset[0].p3.dtype == np.float64 - - -@pytest.mark.parametrize("lonlatdepth_dtype", [np.float64, np.float32]) -def test_pset_create_field(fieldset, lonlatdepth_dtype): - npart = 100 - np.random.seed(123456) - shape = (fieldset.U.lat.size, fieldset.U.lon.size) - K = Field("K", lon=fieldset.U.lon, lat=fieldset.U.lat, data=np.ones(shape, dtype=np.float32)) - pset = ParticleSet.from_field( - fieldset, size=npart, pclass=Particle, start_field=K, lonlatdepth_dtype=lonlatdepth_dtype - ) - assert (np.array([p.lon for p in pset]) <= K.lon[-1]).all() - assert (np.array([p.lon for p in pset]) >= K.lon[0]).all() - assert (np.array([p.lat for p in pset]) <= K.lat[-1]).all() - assert (np.array([p.lat for p in pset]) >= K.lat[0]).all() - assert isinstance(pset[0].lat, lonlatdepth_dtype) - - -def test_pset_create_field_curvi(): - npart = 100 - np.random.seed(123456) - r_v = np.linspace(0.25, 2, 20) - theta_v = np.linspace(0, np.pi / 2, 200) - dtheta = theta_v[1] - theta_v[0] - dr = r_v[1] - r_v[0] - (r, theta) = np.meshgrid(r_v, theta_v) - - x = -1 + r * np.cos(theta) - y = -1 + r * np.sin(theta) - grid = CurvilinearZGrid(x, y) - - u = np.ones(x.shape) - v = np.where(np.logical_and(theta > np.pi / 4, theta < np.pi / 3), 1, 0) - - ufield = Field("U", u, grid=grid) - vfield = Field("V", v, grid=grid) - fieldset = FieldSet(ufield, vfield) - pset = ParticleSet.from_field(fieldset, size=npart, pclass=Particle, start_field=fieldset.V) - - lons = np.array([p.lon + 1 for p in pset]) - lats = np.array([p.lat + 1 for p in pset]) - thetas = np.arctan2(lats, lons) - rs = np.sqrt(lons * lons + lats * lats) - - test = np.pi / 4 - dtheta < thetas - test *= thetas < np.pi / 3 + dtheta - test *= rs > 0.25 - dr - test *= rs < 2 + dr - assert np.all(test) - - -def test_pset_create_with_time(fieldset): - npart = 100 - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) - time = 5.0 - pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle, time=time) - assert np.allclose([p.time for p in pset], time, rtol=1e-12) - pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle, time=[time] * npart) - assert np.allclose([p.time for p in pset], time, rtol=1e-12) - pset = ParticleSet.from_line(fieldset, size=npart, start=(0, 1), finish=(1, 0), pclass=Particle, time=time) - assert np.allclose([p.time for p in pset], time, rtol=1e-12) - - -def test_pset_repeated_release(fieldset): - npart = 10 - time = np.arange(0, npart, 1) # release 1 particle every second - pset = ParticleSet(fieldset, lon=np.zeros(npart), lat=np.zeros(npart), pclass=Particle, time=time) - assert np.allclose([p.time for p in pset], time) - - def IncrLon(particle, fieldset, time): # pragma: no cover - particle.dlon += 1.0 - - pset.execute(IncrLon, dt=1.0, runtime=npart + 1) - assert np.allclose([p.lon for p in pset], np.arange(npart, 0, -1)) - - -def test_pset_repeatdt_check_dt(fieldset): - pset = ParticleSet(fieldset, lon=[0], lat=[0], pclass=Particle, repeatdt=5) - - def IncrLon(particle, fieldset, time): # pragma: no cover - particle.lon = 1.0 - - pset.execute(IncrLon, dt=2, runtime=21) - assert np.allclose([p.lon for p in pset], 1) # if p.dt is nan, it won't be executed so p.lon will be 0 - - -def test_pset_access(fieldset): - npart = 100 - lon = np.linspace(0, 1, npart, dtype=np.float32) - lat = np.linspace(1, 0, npart, dtype=np.float32) - pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle) - assert pset.size == 100 - assert np.allclose([pset[i].lon for i in range(pset.size)], lon, rtol=1e-12) - assert np.allclose([pset[i].lat for i in range(pset.size)], lat, rtol=1e-12) - - -def test_pset_custom_pclass(fieldset): - npart = 100 - TestParticle = Particle.add_variable([Variable("p", np.float32, initial=0.33), Variable("n", np.int32, initial=2)]) - - pset = ParticleSet(fieldset, pclass=TestParticle, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart)) - assert pset.size == npart - assert np.allclose([p.p - 0.33 for p in pset], np.zeros(npart), atol=1e-5) - assert np.allclose([p.n - 2 for p in pset], np.zeros(npart), rtol=1e-12) - - -def test_pset_add_execute(fieldset): - npart = 10 - - def AddLat(particle, fieldset, time): # pragma: no cover - particle.dlat += 0.1 - - pset = ParticleSet(fieldset, lon=[], lat=[], pclass=Particle) - for _ in range(npart): - pset += ParticleSet(pclass=Particle, lon=0.1, lat=0.1, fieldset=fieldset) - for _ in range(4): - pset.execute(pset.Kernel(AddLat), runtime=1.0, dt=1.0) - assert np.allclose(np.array([p.lat for p in pset]), 0.4, rtol=1e-12) - - -@pytest.mark.xfail(reason="Particle removal has not been implemented yet") -def test_pset_remove_particle(fieldset): - npart = 100 - lon = np.linspace(0, 1, npart) - lat = np.linspace(1, 0, npart) - pset = ParticleSet(fieldset, lon=lon, lat=lat, pclass=Particle) - for ilon, ilat in zip(lon[::-1], lat[::-1], strict=True): - assert pset.lon[-1] == ilon - assert pset.lat[-1] == ilat - pset.remove_indices(pset[-1]) - assert pset.size == 0 - - -@pytest.mark.parametrize("staggered_grid", ["Agrid", "Cgrid"]) -def test_from_field_exact_val(staggered_grid): - """ - Tests the creation of a ParticleSet from a field with exact values - on both A-grid and C-grid staggered grids. Verifies that particles - are initialized correctly within the masked region and that their - properties match the expected field values. - """ - xdim = 4 - ydim = 3 - - lon = np.linspace(-1, 2, xdim, dtype=np.float32) - lat = np.linspace(50, 52, ydim, dtype=np.float32) - - dimensions = {"lat": lat, "lon": lon} - if staggered_grid == "Agrid": - U = np.zeros((ydim, xdim), dtype=np.float32) - V = np.zeros((ydim, xdim), dtype=np.float32) - data = {"U": np.array(U, dtype=np.float32), "V": np.array(V, dtype=np.float32)} - mask = np.array([[1, 1, 0, 0], - [1, 1, 1, 0], - [1, 1, 1, 1]]) # fmt: skip - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - - FMask = Field("mask", mask, lon, lat) - fieldset.add_field(FMask) - elif staggered_grid == "Cgrid": - U = np.array([[0, 0, 0, 0], - [1, 0, 0, 0], - [1, 1, 0, 0]]) # fmt: skip - V = np.array([[0, 1, 0, 0], - [0, 1, 0, 0], - [0, 1, 1, 0]]) # fmt: skip - data = {"U": np.array(U, dtype=np.float32), "V": np.array(V, dtype=np.float32)} - mask = np.array([[-1, -1, -1, -1], [-1, 1, 0, 0], [-1, 1, 1, 0]]) - fieldset = FieldSet.from_data(data, dimensions, mesh="flat") - fieldset.U.interp_method = "cgrid_velocity" - fieldset.V.interp_method = "cgrid_velocity" - - FMask = Field("mask", mask, lon, lat, interp_method="cgrid_tracer") - fieldset.add_field(FMask) - - SampleParticle = Particle.add_variable("mask", initial=0) - - def SampleMask(particle, fieldset, time): # pragma: no cover - particle.mask = fieldset.mask[particle] - - pset = ParticleSet.from_field(fieldset, size=400, pclass=SampleParticle, start_field=FMask, time=0) - pset.execute(SampleMask, dt=1, runtime=1) - assert np.allclose([p.mask for p in pset], 1) - assert (np.array([p.lon for p in pset]) <= 1).all() - test = np.logical_or(np.array([p.lon for p in pset]) <= 0, np.array([p.lat for p in pset]) >= 51) - assert test.all() diff --git a/tests-v3/tools/test_warnings.py b/tests-v3/tools/test_warnings.py deleted file mode 100644 index 7dc03c66ad..0000000000 --- a/tests-v3/tools/test_warnings.py +++ /dev/null @@ -1,43 +0,0 @@ -import numpy as np -import pytest - -from parcels import ( - AdvectionRK45, - FieldSet, - FieldSetWarning, - KernelWarning, - Particle, - ParticleSet, -) -from tests.utils import TEST_DATA - - -@pytest.mark.v4alpha -@pytest.mark.xfail(reason="From_pop is not supported during v4-alpha development. This will be reconsidered in v4.") -def test_fieldset_warning_pop(): - filenames = str(TEST_DATA / "POPtestdata_time.nc") - variables = {"U": "U", "V": "V", "W": "W", "T": "T"} - dimensions = {"lon": "lon", "lat": "lat", "depth": "w_deps", "time": "time"} - with pytest.warns(FieldSetWarning, match="General s-levels are not supported in B-grid.*"): - # b-grid with s-levels and POP output in meters warning - FieldSet.from_pop(filenames, variables, dimensions, mesh="flat") - - -def test_kernel_warnings(): - # RK45 warnings - lat = [0, 1, 5, 10] - lon = [0, 1, 5, 10] - u = [[1, 1, 1, 1] for _ in range(4)] - v = [[1, 1, 1, 1] for _ in range(4)] - fieldset = FieldSet.from_data(data={"U": u, "V": v}, dimensions={"lon": lon, "lat": lat}) - pset = ParticleSet( - fieldset=fieldset, - pclass=Particle.add_variable("next_dt", dtype=np.float32, initial=1), - lon=[0], - lat=[0], - depth=[0], - time=[0], - next_dt=1, - ) - with pytest.warns(KernelWarning): - pset.execute(AdvectionRK45, runtime=1, dt=1) diff --git a/tests-v3/tools/test_helpers.py b/tests/test_decorators.py similarity index 79% rename from tests-v3/tools/test_helpers.py rename to tests/test_decorators.py index c3499b55a1..060bde6d7c 100644 --- a/tests-v3/tools/test_helpers.py +++ b/tests/test_decorators.py @@ -1,17 +1,6 @@ import pytest -import parcels.tools._helpers as helpers -from parcels.tools._helpers import deprecated, deprecated_made_private - - -def test_format_list_items_multiline(): - expected = """[ - item1, - item2, - item3 -]""" - assert helpers._format_list_items_multiline(["item1", "item2", "item3"], 1) == expected - assert helpers._format_list_items_multiline([], 1) == "[]" +from parcels._decorators import deprecated, deprecated_made_private def test_deprecated(): diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index a97e2cd879..49101320c1 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -11,6 +11,7 @@ StatusCode, Variable, VectorField, + read_particlefile, ) from parcels._core.index_search import _search_time_index from parcels._datasets.structured.generated import simple_UV_dataset @@ -193,22 +194,8 @@ def test_interpolation_mesh_type(mesh, npart=10): assert v == 0.0 -interp_methods = { - "linear": XLinear, -} - - -@pytest.mark.xfail(reason="ParticleFile not implemented yet") -@pytest.mark.parametrize( - "interp_name", - [ - "linear", - # "freeslip", - # "nearest", - # "cgrid_velocity", - ], -) -def test_interp_regression_v3(interp_name): +@pytest.mark.parametrize(("interp_name", "interp_method"), [("linear", XLinear)]) +def test_interp_regression_v3(interp_name, interp_method): """Test that the v4 versions of the interpolation are the same as the v3 versions.""" ds_input = xr.open_dataset(str(TEST_DATA / f"test_interpolation_data_random_{interp_name}.nc")) ydim = ds_input["U"].shape[2] @@ -247,20 +234,20 @@ def test_interp_regression_v3(interp_name): ) fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") - assert fieldset.U.interp_method == interp_methods[interp_name] - assert fieldset.V.interp_method == interp_methods[interp_name] - assert fieldset.W.interp_method == interp_methods[interp_name] + assert isinstance(fieldset.U.interp_method, interp_method) + assert isinstance(fieldset.V.interp_method, interp_method) + assert isinstance(fieldset.W.interp_method, interp_method) x, y, z = np.meshgrid(np.linspace(0, 1, 7), np.linspace(0, 1, 13), np.linspace(0, 1, 5)) TestP = Particle.add_variable(Variable("pid", dtype=np.int32, initial=0)) pset = ParticleSet(fieldset, pclass=TestP, x=x, y=y, z=z, pid=np.arange(x.size)) - def DeleteParticle(particle, fieldset, time): - if particle.state >= 50: - particle.state = StatusCode.Delete + def DeleteParticle(particles, fieldset): + any_error = particles.state >= 50 # This captures all Errors + particles[any_error].state = StatusCode.Delete - outfile = ParticleFile(f"test_interpolation_v4_{interp_name}", outputdt=np.timedelta64(1, "s")) + outfile = ParticleFile(f"test_interpolation_v4_{interp_name}.parquet", outputdt=np.timedelta64(1, "s"), mode="w") pset.execute( [AdvectionRK4_3D, DeleteParticle], runtime=np.timedelta64(4, "s"), @@ -270,9 +257,43 @@ def DeleteParticle(particle, fieldset, time): print(str(TEST_DATA / f"test_interpolation_jit_{interp_name}.zarr")) ds_v3 = xr.open_zarr(str(TEST_DATA / f"test_interpolation_jit_{interp_name}.zarr")) - ds_v4 = xr.open_zarr(f"test_interpolation_v4_{interp_name}.zarr") + ds_v4 = read_particlefile(f"test_interpolation_v4_{interp_name}.parquet") + + v3_starts = np.column_stack([ds_v3.lon[:, 0].values, ds_v3.lat[:, 0].values, ds_v3.z[:, 0].values]) + unique_starts_v3, inverse_indices = np.unique(v3_starts, axis=0, return_inverse=True) + + v4_pid_to_data = {pid: ds_v4.filter(ds_v4["particle_id"] == pid) for pid in ds_v4["particle_id"].unique()} + + for start_lon, start_lat, start_z in unique_starts_v3: + # Find particles in v3 with this starting position + ind_v3 = np.where( + inverse_indices + == np.where( + (unique_starts_v3[:, 0] == start_lon) + & (unique_starts_v3[:, 1] == start_lat) + & (unique_starts_v3[:, 2] == start_z) + )[0][0] + )[0][0] + + # Find particles in v4 with this starting position using vectorized filter + v4_mask = (ds_v4["x"] == start_lon) & (ds_v4["y"] == start_lat) & (ds_v4["z"] == start_z) + ind_v4 = ds_v4.filter(v4_mask)["particle_id"].unique().to_numpy() + + v3_lon = ds_v3.lon[ind_v3, :].values + v3_lat = ds_v3.lat[ind_v3, :].values + v3_z = ds_v3.z[ind_v3, :].values + + # Use cached v4 data + v4_data = v4_pid_to_data[ind_v4[0]] + v4_lon = v4_data["x"].to_numpy()[:-1] + v4_lat = v4_data["y"].to_numpy()[:-1] + v4_z = v4_data["z"].to_numpy()[:-1] + + # Skip if all NaN + if np.all(np.isnan(v3_lon)) or np.all(np.isnan(v4_lon)): + continue - tol = 1e-6 - np.testing.assert_allclose(ds_v3.lon, ds_v4.lon, atol=tol) - np.testing.assert_allclose(ds_v3.lat, ds_v4.lat, atol=tol) - np.testing.assert_allclose(ds_v3.z, ds_v4.z, atol=tol) + tol = 1e-6 + np.testing.assert_allclose(v3_lon, v4_lon, atol=tol) + np.testing.assert_allclose(v3_lat, v4_lat, atol=tol) + np.testing.assert_allclose(v3_z, v4_z, atol=tol) diff --git a/tests/test_kernel.py b/tests/test_kernel.py index 7f9b139337..dd6db2d00d 100644 --- a/tests/test_kernel.py +++ b/tests/test_kernel.py @@ -2,11 +2,16 @@ import pytest from parcels import ( + FieldSet, + KernelWarning, + Particle, ParticleSet, + Variable, ) from parcels._core.kernel import Kernel +from parcels._datasets.structured.generated import simple_UV_dataset from parcels.kernels import AdvectionRK4, AdvectionRK45 -from tests.common_kernels import MoveEast, MoveNorth +from tests.common_kernels import DoNothing, MoveEast, MoveNorth def test_unknown_var_in_kernel(fieldset): @@ -107,6 +112,18 @@ def test_RK45Kernel_error_no_next_dt(fieldset): Kernel(kernels=[AdvectionRK45], pset=pset) +def test_rk45_kernel_warnings(fieldset): + pset = ParticleSet( + fieldset=fieldset, + pclass=Particle.add_variable(Variable("next_dt", dtype=np.float32, initial=1)), + x=[0], + y=[0], + next_dt=1, + ) + with pytest.warns(KernelWarning): + pset.execute(AdvectionRK45, runtime=1, dt=1) + + def test_kernel_signature(fieldset): pset = ParticleSet(fieldset, x=[0.5], y=[0.5]) @@ -145,3 +162,66 @@ def kernel_with_forced_kwarg(particles, *, fieldset=0): match="Parameter 'fieldset' has incorrect parameter kind. Expected POSITIONAL_OR_KEYWORD, got KEYWORD_ONLY", ): Kernel(kernels=[kernel_with_forced_kwarg], pset=pset) + + +@pytest.mark.parametrize("kernel_type", ["update_lon", "update_dlon"]) +def test_execution_order(kernel_type): + ds = simple_UV_dataset(dims=(1, 1, 2, 2), mesh="flat") + ds["U"].data[:, :] = [[0, 1], [2, 3]] + ds["lon"].data = [0, 2] + fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat") + + def MoveLon_Update_X(particles, fieldset): # pragma: no cover + particles.x += 0.2 + + def MoveLon_Update_DX(particles, fieldset): # pragma: no cover + particles.dx += 0.2 + + def SampleP(particles, fieldset): # pragma: no cover + particles.p = fieldset.U[particles] + print(particles.x, particles.p, fieldset.U[particles]) + + SampleParticle = Particle.add_variable(Variable("p", dtype=np.float32, initial=0.0)) + + MoveLon = MoveLon_Update_DX if kernel_type == "update_dlon" else MoveLon_Update_X + + kernels = [MoveLon, SampleP] + lons = [] + ps = [] + for dir in [1, -1]: + pset = ParticleSet(fieldset, pclass=SampleParticle, x=0, y=0) + pset.execute(kernels[::dir], runtime=1, dt=1) + lons.append(pset.x) + ps.append(pset.p) + + if kernel_type == "update_dlon": + assert np.isclose(lons[0], lons[1]) + assert np.isclose(ps[0], ps[1]) + assert np.allclose(lons[0], 0.2) + else: + assert np.isclose(ps[0] - ps[1], 0.1) + assert np.allclose(lons[0], 0.2) + + +@pytest.mark.xfail(reason="Modifying dt in a kernel doesn't work GH2765") +def test_dt_modify_in_kernel(fieldset): + TestParticle = Particle.add_variable(Variable("age", dtype=np.float32, initial=0)) + pset = ParticleSet(fieldset, pclass=TestParticle, x=[0.5], y=[0]) + + def ModifyDt(particles, fieldset): # pragma: no cover + particles.age += particles.dt + particles.dt = 2 + + runtime = 10 + expected_age = 1 + 2 * (runtime - 2) # 1 for the first step; 2 for the remaining steps (except last) + pset.execute(ModifyDt, runtime=runtime, dt=1.0) + np.testing.assert_allclose(pset.t[0], runtime) + np.testing.assert_allclose(pset.age[0], expected_age) + + +@pytest.mark.parametrize("dt", [1e-2, 1e-5, 1e-6, 1e-9]) +def test_small_dt(fieldset, dt): + pset = ParticleSet(fieldset, x=[0], y=[0]) + + pset.execute(DoNothing, dt=dt, runtime=dt * 100) + assert np.allclose([p.t for p in pset], dt * 100) diff --git a/tests/test_particle.py b/tests/test_particle.py index 62eb65cffc..fe17130e1e 100644 --- a/tests/test_particle.py +++ b/tests/test_particle.py @@ -117,6 +117,13 @@ def test_particleclass_add_variable_in_loop(): assert var1.to_write == var2.to_write +def test_variable_special_names(fieldset): + """Test that checks errors thrown for special names.""" + for vars in ["z", "x"]: + with pytest.raises(ValueError): + Particle.add_variable(Variable(vars, dtype=np.float32, initial=10.0)) + + def test_particleclass_add_variable_collision(): p_initial = ParticleClass(variables=[Variable("vara", dtype=np.float32)]) diff --git a/tests/test_repr_utils.py b/tests/test_repr_utils.py new file mode 100644 index 0000000000..37d8342684 --- /dev/null +++ b/tests/test_repr_utils.py @@ -0,0 +1,11 @@ +from parcels._repr_utils import _format_list_items_multiline + + +def test_format_list_items_multiline(): + expected = """[ + item1, + item2, + item3 +]""" + assert _format_list_items_multiline(["item1", "item2", "item3"], 1) == expected + assert _format_list_items_multiline([], 1) == "[]" diff --git a/tests/test_sigmagrids.py b/tests/test_sigmagrids.py index 313dc14a4d..8337091eb0 100644 --- a/tests/test_sigmagrids.py +++ b/tests/test_sigmagrids.py @@ -3,7 +3,7 @@ import parcels import parcels.tutorial from parcels import Particle, ParticleSet, Variable -from parcels.kernels import AdvectionRK2_3D_CROCO, SampleOmegaCroco, convert_z_to_sigma_croco +from parcels.kernels import AdvectionRK2, AdvectionRK2_3D_CROCO, SampleOmegaCroco, convert_z_to_sigma_croco def test_conversion_3DCROCO(): @@ -44,6 +44,28 @@ def test_conversion_3DCROCO(): np.testing.assert_allclose(sigma, s_xroms, atol=1e-3) +def test_advection_2DCROCO(): + ds_fields = parcels.tutorial.open_dataset("CROCOidealized_data/data") + + fields = { + "U": ds_fields["u"], + "V": ds_fields["v"], + } + ds_fset = parcels.convert.croco_to_sgrid(fields=fields, coords=ds_fields) + fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset) + fieldset = fieldset.to_windowed_arrays() + + runtime = 10_000 + X = np.array([40e3, 80e3, 120e3]) + Y = np.ones(X.size) * 100e3 + Z = np.zeros(X.size) + pset = ParticleSet(fieldset=fieldset, x=X, y=Y, z=Z) + + pset.execute([AdvectionRK2], runtime=runtime, dt=100) + assert np.allclose(pset.z, Z.flatten(), atol=1e-3) + assert np.allclose(pset.x, [x + runtime for x in X], atol=1e-3) + + def test_advection_3DCROCO(): ds_fields = parcels.tutorial.open_dataset("CROCOidealized_data/data")